1. Introduction
Internal solitons in the ocean (ISWs) are, arguably, the most commonly observed kind of solitary waves (or structures close to them) in nature. Despite a relatively long history of such observations, many related theoretical problems remain unsolved. It refers, among others, to strongly nonlinear solitons that cannot be described by the classical Korteweg-de Vries equation and its numerous modifications. Strongly nonlinear solitons are rather common in the ocean, particularly on oceanic shelves where they are generated by the barotropic tide. There are many works on the observation and description of soliton-like waves in the shelf regions [Reference Apel, Ostrovsky, Stepanyants and Lynch1–Reference Orr and Mignerey7]. An even more complex factor is that in the oceanic conditions, internal solitons always interact with other motions of different scales, including currents, turbulence and surface waves. All that contributes to the ‘internal weather’ and, in particular, to the dynamics of biological components and impurities responsible for the ecological environment. Here we limit ourselves to one, albeit broad, class of processes: the interactions between internal waves (IWs) and small-scale turbulence. This, in turn, can be conditionally broken into several problems. The one is the attenuation of internal waves due to turbulent viscosity and diffusion [Reference Monin and Yaglom8–Reference Ostrovsky11]. Another is how internal waves generate turbulence, which is observed almost everywhere in the ocean, including its deep areas. It is often associated with shear instability that, according to the classical Miles–Howard theory, is possible if Ri < ¼, where Ri is the Richardson number [Reference Miles12, Reference Howard13], and the eventual generation of turbulence. At the same time, turbulence is regularly observed for larger values of Ri. For a developed turbulence when the semi-empirical description is valid, the turbulence can be maintained under a softer condition, Ri < 1. However, in many cases, turbulence exists at Ri > 1, i.e., for a much stronger stratification [Reference Moum, Farmer, Smyth, Armi and Vagle14–Reference Palmer, Stephenson, Inall, Balfour, Düsterhus and Green19]. The explanation of that is based on a modification of the classical semi-empirical theory of turbulence suggested in Ostrovsky and Troitskaya [Reference Ostrovsky and Troitskaya20], where the kinetic equation for the distribution function of velocity and density fluctuations was used. This approach reduces the number of arbitrary constants in the semi-empirical equations and, more importantly, it adds an equation for the potential energy of turbulence, which is proportional to the variance of density fluctuations. As a result, the restriction on the value of Ri can be completely removed so that turbulence energy can remain finite at any finite velocity shear. The application of this theory to the turbulence dynamics under the action of shear flow in different areas of the ocean was studied in refs [Reference Ostrovsky and Troitskaya20–Reference Ostrovsky, Gladskikh and Soustova24]. The case of internal waves is similar but more complicated because the results depend not only on the IW amplitude but also on the ratio between the time scales of the wave and those of generation and relaxation of turbulent energy. There exists extensive literature regarding turbulence generation by instability and breaking of internal waves (e.g. [Reference D'Asaro and Lien25–Reference Rodda27], and references therein). Here we are interested in the effect of non-breaking IW supporting turbulence long after a breaking event. This effect was observed in the laboratory experiment [Reference Matusov, Ostrovsky and Tsimring28], where a standing internal wave amplified turbulence penetrating downwards from a vibrating perforated grid. For the oceanic observations, we discuss below the effect of strong internal solitons on turbulence at large Richardson numbers.
Theoretical descriptions of stratified turbulence and the corresponding numerical modeling methods, besides the classical semi-empirical (gradient) schemes (RANS), use Large Eddy Simulation (LES) and direct numerical simulation (DNS). The difference between these schemes is based on the degree of resolution of turbulent pulsations and their energy spectrum. Over the past 30 years, new approaches have been developed to study the characteristics of stratified turbulence based on the spectral approach (e.g [Reference Galperin and Sukoriansky29]). These models potentially allow for a more consistent description of the unsteady processes of interaction between turbulence and internal waves than RANS models. More references can be found in Gladskikh et al. [Reference Gladskikh, Ostrovsky, Troitskaya, Soustova and Mortikov23]. For the problems considered here, we refer to the direct numerical simulation of the IW action on turbulence [Reference Druzhinin and Ostrovsky30].
Returning to the theme of this paper, we consider the model mentioned above, where the given mean fluid velocity field u(r, t) depends on time due to an internal wave. In these cases, the general equations derived in [Reference Ostrovsky and Troitskaya20] (see also [Reference Ostrovsky, Soustova, Troitskaya and Gladskikh21, Reference Gladskikh, Ostrovsky, Troitskaya, Soustova and Mortikov23]) are reduced to the equations for the turbulent kinetic energy (TKE, K) and turbulent potential energy (P) in the form [Reference Monin and Yaglom8, Reference Ostrovsky and Troitskaya20]:
\begin{align}
\frac{{\partial K}}{{\partial t}} + \left( {{\mathbf{u}}\nabla } \right)K &= L\sqrt K {\left( {\frac{{\partial {u_i}}}{{\partial {x_k}}} + \frac{{\partial {u_k}}}{{\partial {x_i}}}} \right)^2} - {N^2}L\sqrt K \left[ {1 - \frac{{3P}}{K}(1 - G)} \right] - \frac{{C{K^{3/2}}}}{L} + \frac{5}{3}\frac{\partial }{{\partial {x_i}}}\left( {L\sqrt K \frac{{\partial K}}{{\partial {x_i}}}} \right), \nonumber\\
\frac{{\partial P}}{{\partial t}} + \left( {{\mathbf{u}}\nabla } \right)P &= {N^2}L\sqrt K \left[ {1 - \frac{{3P}}{K}(1 - G)} \right] - \frac{{D{K^{1/2}}}}{L}P + \frac{\partial }{{\partial {x_i}}}\left( {\frac{{L\sqrt K }}{{\rho {N^2}}}\frac{{\partial (\rho {N^2}P)}}{{\partial {x_i}}}} \right).
\end{align} Here, 
${N^2}(z) = (g/\rho )d\rho /dz$ is the squared Brunt–Vaisala (buoyancy) frequency, g is the gravity acceleration, and ρ(z) is the mean density. C and D are empirical constants (in what follows, they are taken C = D = 0.09), L is the outer scale of turbulence, and G is the anisotropy parameter, which depends on the ratio s = 
$\frac{L_z}{L_x}$ of the vertical to horizontal correlation scales of turbulence. As shown in Ostrovsky and Troitskaya [Reference Ostrovsky and Troitskaya20], if s ≪1 (strongly anisotropic “pancake” turbulence), then G is close to 1, and the common gradient theory, neglecting potential energy, is valid. In the quasi-isotropic case, typical of small-scale turbulence, s is close to 1 and G significantly differs from 1, so that equations (1) must be solved together. Here, G = 0.5 is taken.
2. Turbulence affected by a periodic IW
Some characteristic features of the process can be elucidated from considering the effect of a sinusoidal (linear or weakly nonlinear) internal wave mode; note that equations (1) for turbulence remain nonlinear. A single IW mode with the velocity components uz = W, ux = U for such waves has the form
\begin{equation}W = A\sin ({k_z}z)\cos (\omega t - {k_x}x),{\text{ }}U = A({k_z}/{k_x})\cos ({k_z}z)\sin (\omega t - {k_x}x).\end{equation}Figure 1 shows two cases, both with the last terms responsible for diffusion neglected in (1). Here, equations (2) are local, and the depth z enters as a parameter.
Effect of the internal waves of different frequencies on a turbulent layer in a given vertical cross-section (x = 0). Parameters in (2) are x = 0, and z = 5 m, L = 1. (a) A = 0.02 m/s, ω = 0.04 rad/s, N = 0.049 rad/s; (b) A = 0.01 m/s, ω = 0.0019 rad/s, N = 0.01645 rad/s. Blue: kinetic, orange: potential energy. The dashed line in the bottom plot is for the classic model without potential energy (only the first equation (1) with G = 1).
Figure 1 illustrates two qualitatively different kinds of turbulence evolution. For a high-frequency internal wave, the turbulent energy reaches a constant level with relatively small oscillations at a frequency 
$2\omega $ around it. For a longer wave, turbulence energy follows the wave energy. Note that within the classical theory, such a process, among others, was briefly considered by Ivanov et al. [Reference Ivanov, Ostrovsky, Soustova and Tsimring31]. In that case, the turbulence level can drop to practically zero (dashed line in the bottom plot) when the semi-empirical approach becomes inapplicable. In the system (1), where the turbulent potential energy is significant, such an effect is less pronounced, as there is no critical level of Ri.
The above results are valid if the vertical scale of the wave mode is large (small kz). Otherwise, the diffusion terms in (1) are important. Figure 2 illustrates their effect.
In these cases, the saturation of the energy level occurs more slowly, and in the high-frequency case (a), an additional modulation with a period of about 4000 remains.
3. Effect of solitons
3. 1. The post-soliton stage
The action of localized internal waves, including solitary waves, on turbulence depends on the relationship between the duration of the soliton and the characteristic time of turbulence variation, so that it may either follow the wave profile in a quasi-stationary way or react with a delay. Rather universal is the turbulence decay after the soliton impact, when u = 0. Due to the Kolmogorov-type dissipation, the values of K and P become small. At this stage, the last two terms in both equations (1) are of the order of K 3/2, whereas the terms with N 2 are of the order of K 1/2. Hence, asymptotically, the terms with N 2 would dominate. Since these terms have opposite signs for positive K and P, the two equations (1) can be compatible only if these terms go to zero. As a result, asymptotically we have
\begin{equation}3P = K/(1 - G),{\text{ }}G \ne 1.\end{equation}If (3) is met, equations (1) are separated, and the solution for K at K(t 0) = K 0 is
\begin{equation}K(t \geq {t_0}) = \frac{{{K_0}}}{{{{\left[ {1 + C\sqrt {{K_0}} (t - {t_0})/2L} \right]}^2}}}.\end{equation}As will be seen below, these simple relations can be true even shortly after the soliton.
3.2. The effect of Gardner solitons
As the first example, consider the action of solitons, which are the solutions of the classic Gardner equation (e.g. [Reference Apel, Ostrovsky, Stepanyants and Lynch1]). The vertical shear of the horizontal fluid velocity in a wave mode corresponding to (1) has the form uz = U (t, x ) fz (z). In the dimensionless variables, the factor U satisfies the equation
\begin{equation}{U_t} + {U_x} + 6(U - {U^2}){U_x} + {U_{xxx}} = 0.\end{equation}Here (5) is written in the original variables x and t, adding the velocity of a long linear wave, c 0 = 1. In physics, the corresponding equation is still weakly nonlinear, but unlike the KdV solitons, its solitary solutions can broaden with amplitude, up to its limiting value when the soliton becomes flat-top. This is also characteristic of strongly nonlinear solitons [Reference Apel, Ostrovsky, Stepanyants and Lynch1, Reference Ostrovsky and Grue2]. The solitary solutions of (5) are [Reference Apel, Ostrovsky, Stepanyants and Lynch1]:
\begin{align}
U &= \frac{k}{2}\left[ {\tanh \frac{k}{2}\left( {x - t - {k^2}t + \Delta } \right) - \tanh \frac{k}{2}\left( {x - t - {k^2}t - \Delta } \right)} \right], \nonumber\\
\Delta &= {k^{ - 1}}\text{arctanh} (k),{\text{ }}0 \leq k \leq 1.
\end{align}At small k, this solution is close to the KdV soliton, whereas near k = 1, it defines a flat-top soliton. The examples are plotted in Figure 3.
The solitons (6) for k = 0.6, 0.95, and 1–10−8 (from the smaller to the larger).
Consider first the local, non-diffusive model. In eqs. (1), we let N 2 = 0.25 and (df/dz)2 = 0.2 at a given level of z. The local Richardson number is Ri = N 2/uz 2 = 1.25/U 2. Figure 4 shows the time dependence of Ri for three values of the parameter k in (6). It exceeds unity for all times, so that the classic gradient model would prevent any support of turbulence by such a wave.
Variation of local Richardson number for the solitons shown in Figure 3, with the same values of k.
The corresponding solutions of equations (1) are shown in Figure 5. As in the previous section, the local, non-diffusive model is valid if the vertical scale of the velocity field is sufficiently large. The diffusive case was calculated for a mode proportional to 
$\cos ({k_z}z)$ with kz = 0.1π, as in the case of a sinusoidal wave.
Evolution of kinetic and potential energies at z = 5. From lower to upper pairs of curves: k = 0.6, 0.95, and 1–10–8. Blue-kinetic, maroon-potential energy of turbulence. (a) Without diffusion, (b) with diffusion.
Again, because of the effect of potential energy, turbulence exists at all Richardson numbers. Note that, in agreement with (3), here K = 1.5 P for, roughly, 
$t \geq 10$. The inertial character of the process is evident, and the solutions (3) and (4) are valid almost immediately after the soliton. Figure 6 shows a 3D plot showing the space-time behavior of the kinetic energy of turbulence under the action of a flat-top soliton.
3D plot of kinetic turbulent energy for k = 1–10−8.
Note that, in this case, the diffusive terms in eqs. (1) significantly diminish the level of turbulent energy.
3.3. Effect of strong internal solitons on turbulence in the ocean
Applications to real field data in the ocean and the atmosphere are more complicated due to the simultaneous action of different uncontrolled factors. As a result, the available publications rarely present sufficient quantitative data to be used in the theory. Here the model (1) is applied to the oceanic data described by Moum et al. [Reference Moum, Farmer, Smyth, Armi and Vagle14], where the effect of strongly nonlinear internal solitons on small-scale turbulence was observed off the North-Western Pacific coast of the USA. Even though rather scarce quantitative data are given, the figures illustrating the effect allow us to make reasonable estimates. Some preliminary results have been illustrated in our short presentation [Reference Ostrovsky, Gladskikh and Soustova24]. Note also that some topics of the present paper have been briefly outlined in the review paper of one of the authors published earlier [Reference Ostrovsky11] with a reference to this paper as being prepared. Here we further study this problem.
For the time dependence of soliton parameters, we use the echosounder images of the leading soliton in the observed group shown in several figures of [Reference Moum, Farmer, Smyth, Armi and Vagle14]. In its Figure 7, the plots are supplemented by isolines of density (isopycnals), color images of fluid velocity, and turbulent kinetic energy dissipation rate. The latter is a commonly measured turbulence characteristic in the ocean, and where possible, we shall verify theoretical results by comparison with the data of the turbulent dissipation rate.
Time profile of the leading soliton, digitized (dashed) and interpolated (solid) from the echosounder image shown in [Reference Moum, Farmer, Smyth, Armi and Vagle14]. The black dot marked in Figure 9 of that paper shows the approximate position of the contact device, measuring the wave vertical profile as shown in their Figure 15.
Figure 7 here shows the time profile of the leading solitary wave in the group.
Note that the corresponding isopycnal is depressed from about 10 m to about 36 m, which testifies to a very strong nonlinearity. Moreover, according to Figure 2 of [Reference Moum, Farmer, Smyth, Armi and Vagle14], the total water depth in the observation area is 100 m or slightly more. Therefore, the wave moves down a significant part of the total water layer.
3.4. Quasistatic approximation
Consider first the quasi-static approximation, supposing that the soliton duration is longer than the time needed for saturation of the turbulence parameters. For the field experiment (Figure 15 of [Reference Moum, Farmer, Smyth, Armi and Vagle14]), vertical profiles of horizontal velocity and density are shown for a vertical cross-section marked by the dot in Figure 7. We digitized and interpolated these profiles, as shown in Figure 8. The motion of the density jump (here at about 30 m depth) seen in Figure 8b will be considered further in the paper.
Depth dependencies of horizontal velocity u (left) and normalized excess density 
$\sigma = (\rho - {\rho _0})/{\rho _0}$ (right). Here it is taken 
${\rho _0}$ =1000 kg/m3. Dashed lines: digitized plots of those in [Reference Moum, Farmer, Smyth, Armi and Vagle14]. Solid lines: their polynomial interpolations.
As mentioned in Moum et al. [Reference Moum, Farmer, Smyth, Armi and Vagle14], the data for fluid velocity obtained from ADCP is rough, so further averaging is justified. Then, we found the depth dependence of the functions (du/dz)2 and N 2 entering the system (1), where u is the horizontal velocity. Their variation in depth is shown in Figure 9.
Interpolated depth dependencies of functions (du/dz)2 (a) and N 2 (b). Here, u is the horizontal component of fluid velocity (for a long wave considered here, the vertical velocity variation can be neglected in this context).
Figure 10 shows the resulting profile of the Richardson number Ri = N 2/(du/dz)2. For almost all depths considered (from 15 to 42 m), Ri > 1, and again, the turbulence can be supported by the wave only because of the effect of its potential energy.
Depth dependence of Richardson number Ri = N 2/(du/dz)2.
Substituting (du/dz)2 and N 2 as functions of z into (1), we obtain the depth dependence of kinetic and potential energies as shown in Figure 11. Figure 12 shows the corresponding distribution of the TKE dissipation rate 
$\varepsilon = C{K^{3/2}}/L$, which is a rather common value measured in oceanic experiments as a characteristic of turbulence level [Reference Moum, Farmer, Smyth, Armi and Vagle14, Reference Forryan, Martin, Srokosz, Popova, Painter and Renner15].
Depth dependence of kinetic (solid line) and potential (dashed line) energies in the quasistatic approximation, L = 1 m.
Depth dependence of TKE dissipation rate in the quasistatic approximation.
3.5. Non-stationary processes
The above results qualitatively correspond to the data of [Reference Moum, Farmer, Smyth, Armi and Vagle14]. Indeed, the turbulence level has two maxima, one near the water surface (supposedly because of wind wave breaking), another at a depth of about 38 m, not far from the experimental data. However, there are some discrepancies. In particular, the TKE dissipation rate (about 5.10−5 m2/s3) exceeds by about an order the data of the cited paper indicated by a color bar in its Figure 6. To evaluate the applicability of the quasistatic approach, we considered a transient process by solving non-stationary equations (1) for the same data as those used above. One result is shown in Figure 13.
Growth of turbulent kinetic (left) and potential (right) energies from small initial values to saturation at three depths, 10, 30, and 40 m.
Figure 14 shows the corresponding variation of dissipation rate.
Growth of the turbulence dissipation rate at different depths.
As seen from these figures, the transient process takes about 10–15 min, which is comparable to the duration of the soliton shown in Figure 7. Therefore, to determine the dynamics of turbulence in the field of a strongly nonlinear soliton of internal waves, it is necessary to take into account the time dependence of the squared horizontal velocity shear 
${\left( {\frac{{\partial u}}{{\partial z}}} \right)^2}$and the squared Brunt–Vaisala frequency N2, which are present in system (1).
3.6. Dynamics of turbulence in the pycnocline
To describe the variation of the parameters 
${\left( {\frac{{\partial u}}{{\partial z}}} \right)^2}{\text{and }}$N2 at the pycnocline, we begin with the data shown in Figure 8. As mentioned above, there is a density jump at a depth of 30 m in the density profile (Figure 8, dashed line), where the turbulence is concentrated, as seen from the color scattering images shown in [Reference Moum, Farmer, Smyth, Armi and Vagle14]. According to Figures 6 and 15 of [Reference Moum, Farmer, Smyth, Armi and Vagle14], the velocity u 1 over the pycnocline is almost independent of depth, except for the near-surface area. To find it, we use a relation between u1 and the pycnocline displacement in a strong soliton given in [Reference Ostrovsky and Grue2]:
\begin{equation}{u_1} = \frac{{c(\eta - {h_1})}}{\eta }.\end{equation}Variation of squared vertical shear of the horizontal fluid velocity.
Here, c is the wave velocity (0.6 m/s in this case), h 1 is the thickness of the upper layer (10–12 m in the experiment [Reference Moum, Farmer, Smyth, Armi and Vagle14]) over the pycnocline, and η is the local depth of the pycnocline (this notation is slightly different from that in [Reference Ostrovsky and Grue2], which is due to the choice of the starting point for depth measurements). The temporal profile of the leading soliton η(t) is shown in Figure 7 above. To determine the dependence of the horizontal velocity on time, we assume that the stratification remains unchanged in the pycnocline at the soliton length (about 300 m), and equal to N 2 = 0.0008 s−2. For this purpose, the soliton-like displacement of the pycnocline shown in Figure 7 was approximated as:
\begin{equation} \eta = - 12 + \left( {14{\text{Tanh}}\left[ {0.007\left( {t - 350} \right)} \right] - 14{\text{Tanh}}\left[ {0.0035\left( {t + 350} \right)} \right]} \right)\,\,{\text{m}}{\text{.}}\end{equation} Using (7), we obtain the time dependence of 
${\left( {\frac{{\partial u}}{{\partial z}}} \right)^2}$at the pycnocline in the form
\begin{equation}{\left( {\frac{{\partial u}}{{\partial z}}} \right)^2} = \frac{{({{1.1.10}^{ - 3}}{{(14\operatorname{Tanh} [0.007(t - 350)] - 14\operatorname{Tanh} [0.0035(t + 350)])}^2}}}{{({{1.1.10}^{ - 3}}{{(14\operatorname{Tanh} [0.007(t - 350)] - 14\operatorname{Tanh} [0.0035(t + 350)])}^2}}}{\text{ }}{{\text{s}}^{ - 2}}.\end{equation}It is shown in Figure 15.
The obtained values of N2 and
${\text{ }}{\left( {\frac{{\partial u}}{{\partial z}}} \right)^2}$ are substituted into the local equations for K and P since the diffusion terms are more than three orders of magnitude smaller than the other terms in (1).
Then, using equations (1), the kinetic and potential energies (Figure 16a) and the kinetic energy dissipation rate (Figure 16b) are calculated.
(a) The kinetic (blue) and potential (yellow) turbulent energy densities at the thermocline. (b) TKE dissipation rate variation along the soliton.
The obtained results have good agreement with the data of [Reference Moum, Farmer, Smyth, Armi and Vagle14], particularly with its Figure 7. First, the order of calculated maximum dissipation rate is about 5.10−6 m2/s2, which corresponds to the maximum values shown in the color bar in Figure 7 of [Reference Moum, Farmer, Smyth, Armi and Vagle14]. Note that the level of dissipation rate before the soliton in the same figure, which can roughly be taken for background, is two orders smaller than that. Second, due to the delay of turbulence development, the maximums of energy and dissipation rate are shifted towards the rear part of the soliton, which also agrees with the sound scattering intensity distribution shown in [Reference Moum, Farmer, Smyth, Armi and Vagle14].
Note that the authors of [Reference Moum, Farmer, Smyth, Armi and Vagle14] explain the existence of turbulence by the presence of microstructure with a presumably small Richardson number, Ri < 1/4. However, as shown here, the observed values and distributions of turbulent energy can be supported by the IW even at Richardson numbers significantly exceeding unity.
4. Conclusions
The modification of the Reynolds-type equations for turbulence described above was suggested as early as 1987. However, its applications to the observed processes in the stratified flows began only recently [Reference Ostrovsky, Soustova, Troitskaya and Gladskikh21–Reference Ostrovsky, Gladskikh and Soustova24]. Here, the effects of the internal waves, particularly internal solitons, have been considered. It is confirmed that the finite-energy turbulence can exist at large Reynolds numbers (strong stratification), and for a strong internal soliton, the theoretical estimates agree with the observational data, both qualitatively and by order of magnitude. This also helps to explain the ubiquitous presence of turbulence in the ocean. Indeed, even when turbulence is generated by wave breaking events, its long-time existence can be due to the support by non-breaking waves which are more common in most areas of the ocean.
Among the promising future developments, there is the description of the mutual action of internal solitons and turbulence when a soliton dissipates during interaction. Attenuation of internal waves on turbulence was studied in detail for sinusoidal internal waves, see a brief review of that in [Reference Ostrovsky11]. For these processes, the variation of turbulence scales in time can be important [Reference Zilitinkevich, Elperin, Kleeorin, Rogachevskii and Esau32]. This work is in progress.
Acknowledgements
The work was supported by the RSF project No. 23-27-00002. I.A.S. and A. M. K. acknowledge the state assignment of the IAP RAS on the topic FFUF-2025-0026.
The authors are grateful to D.S. Gladskikh for valuable discussions.
Open access
