The study of the effect of small-scale turbulence on internal gravity waves propagation in a pycnocline

This paper presents the results of modeling the interaction between internal waves (IWs) and turbulence using direct numerical simulation (DNS). Turbulence is excited and supported by a random forcing localized in a vertical layer separated from the pycnocline. The main attention is paid to the internal wave damping due to turbulence and comparison of the results with those obtained theoretically by using the semi-empirical approach. It is shown that the IW damping rate predicted by the theory agrees well with the DNS results when turbulence is sufficiently strong to be only weakly perturbed by the internal wave; however, the theory overestimates the damping rate of IWs for a weaker turbulence. The DNS parameters are matched to the parameters of the laboratory experiment, and an extrapolation to the oceanic scales is also provided.


Introduction
The interaction between internal gravity waves (IWs) and small-scale turbulence is a significant factor governing the dynamics of the upper ocean.One of the important aspects of this interaction is the phenomenon of IW damping by turbulence (Phillips, 1977).
Early theoretical studies of this phenomenon used sa emi-empirical approach (LeBlond, 1966;Ostrovsky and Soustova, 1979;Ivanov et al., 1983).An explicit expression for the IW damping rate by turbulence based on the semi-empirical closure was derived by Ostrovsky and Zaborskikh (1996) in application to the typical oceanic conditions (vertical profiles of the mean density, turbulent kinetic energy and mean current velocity).
The phenomena of damping of IWs by turbulence was observed in laboratory experiments by Phillips (1977), Kantha (1980), and Barenblatt (1978).A similar phenomena of the damping of surface waves by turbulence generated by a submerged oscillating grid was also experimentally observed by Olmez and Milgram (1992).Quantitative measurements of the IW damping by turbulence were first performed in a laboratory experiment by Ostrovsky et al. (1996).In this experiment, IWs were generated in the pycnocline by a wavemaker, and small-scale turbulence was induced by an oscillating grid at some level above the pycnocline.Measurements and comparison of the IW amplitudes with and without turbulence showed an effective enhancement of the decay rate of IW under the effect of turbulence which was found to be in good agreement with the theoretical prediction by Ostrovsky and Zaborskikh (1996).
However, the known theoretical results have significant limitations.In the theory developed by Ostrovsky and Zaborskikh (1996) it was assumed that the IW is sufficiently weak to allow linearization of the semi-empirical equations, and the experiment by Ostrovsky et al. (1996) was designed correspondingly.A stronger IW can affect the average turbulent energy; moreover, as experimentally observed by Ostrovsky et al. (1989), a strong internal wave can increase the average turbulent energy, which can be a source of the ubiquitous presence of turbulent spots in the ocean.Also, O. A. Druzhinin et al.: The effect of turbulence on IW propagation the applicability of semi-empirical, Reynolds-type equations with a simple closure hypothesis can be limited as well.
In this paper we present the results of direct numerical simulation (DNS) to study of the effect of small-scale turbulence on internal gravity waves propagation in a densitystratified fluid with a pycnocline.As far as we are aware this is the first such study.The DNS results are compared with the semi-empirical model prediction, and explore the applicability limits of the semi-empirical approach concerning the IW damping by turbulence.Thus, for comparison with the theoretical results of Ostrovsky and Zaborskikh (1996), the parameters of the laboratory experiment by Ostrovsky et al. (1996) are employed.Further the geophysical estimates for oceanic scales are also provided.

Basic equations and numerical method
We consider a stably stratified fluid with a pycnocline (Fig. 1).Stationary turbulence is maintained by a forcing located at some distance above the pycnocline.The first mode of the internal wave propagating in the pycnocline from left to right is prescribed as an initial condition.Periodic boundary conditions in the x and y directions and Neumann (zero normal gradient) boundary condition in the z direction are considered.The thickness of the pycnocline, L 0 , and the buoyancy frequency in the middle of the pycnocline, (where g is the gravity acceleration and ρ 0 (z) is the fluid density), are chosen to define the characteristic length and timescales, L 0 and T 0 = 1/N 0 , which are used further to write the governing equations in the dimensionless form.
The Navier-Stokes equations for the fluid velocity are written under the Boussinesq approximation in the dimensionless form as (Phillips, 1977) The equation for the fluid density is written as (3) In Eqs.
(1)-(3), U i (i = x, y, z) is the instantaneous fluid velocity, and ρ and P are the instantaneous deviations of the fluid density and pressure from the respective hydrostatic profiles, x i = x, y, z are the Cartesian coordinates; Re and F r as the Reynolds and Froude numbers are defined as   3) are normalized by the length, time and velocity scales, L 0 , T 0 and U 0 = L 0 /T 0 .Note that since the timescale is defined as T 0 = 1/N 0 and the velocity scale U 0 = L 0 /T 0 = L 0 N 0 , the Froude number in DNS equals (is identical to) unity, F r = 1.The density deviation ρ is normalized by the density jump across the pycnocline, ρ 0 (Fig. 1).
The dimensionless reference profile of the buoyancy frequency, N ref (z), in the left-hand side of Eq. ( 3) is prescribed in the form: where z p defines the pycnocline location.The corresponding dimensionless reference density profile, ρ ref (z), is then obtained from Eq. ( 5) as where constant ρ ref (−∞) can be arbitrary since its value does not influence the integration of Eqs. ( 1)-(3).Thus, for convenience, ρ ref (−∞) is set equal to 1.5, and that the reference density profile is rewritten in the following form: The dimensionless instantaneous density is obtained as a sum of ρ ref (z) and the instantaneous density deviation, ρ.Note that a corresponding dimensional density can be obtained as a sum ρ 0 + ρ 0 0.5[1 − tanh 2(z − z p )] + ρ where ρ 0 is the dimensional reference (undisturbed) density above the pycnocline.
In DNS we prescribe the Reynolds number to be Re = 2000.This is sufficiently large to allow one to consider the viscous damping of IWs as sufficiently weak as compared to the damping caused by turbulence.We also neglect the diffusion effect on N ref (z) and ρ ref (z) profiles, i.e. the profile N ref (z) in Eq. (3) does not evolve with time.The Prandtl number Pr is set equal to unity.
Equations ( 1)-( 3) are discretized in a cubic domain with sizes 0 ≤ x ≤ 40, −10 ≤ y ≤ 10 and 0 ≤ z ≤ 20 by employing a finite difference method of the second-order accuracy on a uniform rectangular staggered grid consisting of 400 × 200 × 200 nodes in the x, y and z directions, respectively.The integration is advanced in time using the Adams-Bashforth method with time step t = 0.01.The Poisson equation for the pressure is solved by FFT transform over x and y coordinates, and Gauss elimination method over z coordinate (Fletcher, 1991;Druzhinin, 2003).The Neumann (zero normal gradient) boundary condition is prescribed for all fields in the vertical (z) direction in the horizontal (x, y) planes at z = 0 and z = 20, and periodic boundary conditions are prescribed in the horizontal (x) and spanwise (y) directions.
In order to induce turbulence in DNS there should be a mechanism to compensate for its decay by the viscous dissipation.In the presence of a mean shear, the nature of the evolving flow may depend critically on the representation and location of turbulence (cf.e.g., Pham et al., 2009).In the present study, there is no mean shear flow, and a random, divergence-free forcing, f i , is employed in the r.h.s. of Eq. (1) to support turbulence.This forcing is considered with the following form: (8) where i = x, y, z, and z f defines the location of the turbulence layer.U f i (x, y, z) is a homogeneous isotropic field with a given power spectrum as follows: where wavenumber k f defines the spectral location of the energy peak.Factor E 0 is chosen so that the amplitude (i.e. an absolute maximum value) of the field U f i (x, y, z) in each direction equals unity.Then parameter F 0 defines the turbulence intensity and hereafter is called the forcing amplitude, and ω f is the forcing frequency.

Internal waves
The initial condition for the velocity and density fields is prescribed as a first mode of internal wave field with wavelength, λ (and wavenumber k = 2π/λ), and frequency, ω.The solution of the linearized Eqs. ( 1)-(3) for the progressive internal wave propagating from left to right in the x direction can be defined as (Phillips, 1977) The initial conditions for IWs field are taken from Eqs. (10)-( 12) at t = 0. Function W (z) in Eqs. ( 10)-( 12) is obtained as an eigenfunction of the well-known boundary problem (Phillips, 1977): , where W 0 is the IW velocity amplitude.The problem in Eq. ( 13) was solved by the shooting method with matching at the pycnocline center (Hazel, 1972).The distribution of the first mode IW vertical velocity and the dispersion relation, ω(k), for wavenumbers in the range 0.3 < k < 6 obtained numerically for four cases of wavelength λ = 4, 5, 8, and 10 are presented in Fig. 2a.The figure shows that, as expected, the energy of the first mode is concentrated around the pycnocline.DNS was performed with initial conditions in Eqs. ( 10  0.531, period T ≈ 12), and λ = 10 (frequency ω = 0.489, period T ≈ 13).The IW amplitude was prescribed as W 0 = 0.05. Figure 2c shows isopycnal displacements obtained in DNS at different times with initial conditions prescribed for IW with wavelength λ = 10.In this case, the amplitude of the isopycnal displacement is about a ≈ 0.1, and the wave slope is about ka = 2π a/λ ≈ 0.06 which may be regarded small enough to ensure that non-linear effects during the IW propagation in the pycnocline remain negligible.Similar results were obtained for IW with wavelengths λ = 4, 5 and 8.

Forced turbulence and its effect on the pycnocline
In order to investigate how turbulence affects the internal wave DNS was performed with only the forcing on and no IW fields imposed initially.The mid-pycnocline level was prescribed at z p = 8 and the forcing level was set at z f = 10.The values of z p and z f were chosen to ensure that the effect of turbulent mixing on the pycnocline remained sufficiently small for the considered forcing amplitudes in the range of 1 < F 0 < 4. On the other hand, the effect of turbulence on the IW propagation should be sufficiently large to be detected in DNS.The frequency of the forcing was prescribed as ω f = 10, which is about 20 times larger than the typical IW frequency, and k f was set equal to unity.This choice is made to ensure that the forced turbulence does not generate effectively internal waves in the pycnocline.At that frequency, the maximum amplitude of density fluctuations induced by turbulence in the pycnocline at forcing amplitude F 0 = 3 is less than 0.01 and negligible as compared to the amplitude of the density oscillations in IWs.On the other hand, the forcing frequency must be much smaller than the inverse time step 1/ t = 100 to ensure sufficient resolution.Figure 3a shows the power spectrum, E(k), of the velocity field obtained in DNS for forcing amplitude F 0 = 3 at different z levels at t = 300.Each spectrum is obtained by FFT transform over the x coordinate and averaged over the y coordinate.The figure shows that the spectra of the velocity fluctuation obtained at different z levels are of the same shape and are reduced in amplitude as the distance (z f − z) increases.The spectra are characterized by local maximum at about k = 1 (in accordance with the forcing spectrum 9), a short inertial interval at k = 7/15 and viscous dissipation region at larger k's.
Using the spectra, the characteristic (integral) turbulence spatial scale was evaluated as and found to be almost independent of z, l t ≈ 0.4 in the region 8 < z < 12.
The mean vertical profiles of the velocity and density fields,< U i > (z) and < ρ > (z), were obtained by averaging over the horizontal (x, y) plane performed for each z.Fluctuations (dispersions) of the velocity and density were then obtained as  Figure 3b shows the vertical mean profiles of the fluctuation velocity, and the mean density obtained in DNS at different times.After a short (as compared to the IW period) transient time interval, the flow becomes statistically stationary.This transition is of no importance for the further process.Once turbulence becomes stationary, the transient effects can be neglected.The figure shows also that the turbulence fluctuation profile is asymmetric about the level z f = 10.This asymmetry of the fluctuation velocity profile is related to the effect of the pycnocline on turbulence diffusion.Downward turbulent diffusion is not effective since turbulent velocity fluctuations can not penetrate the pycnocline at the considered, relatively low, Froude number (F r = 1), whereas turbulence is free to diffuse upwards.Figure 3b also compares the mean density profile obtained in DNS at time t = 300 to the initial (reference) density profile ρ ref (z) in Eq. ( 7).The figure shows that the effect of turbulent mixing on the pycnocline remains negligible.On the other hand, as it will be shown below in Sect.5, the turbulence (with maximum u ≈ 0.06 at z = z f = 10 and u ≈ 0.06 in the vicinity of the pycnocline at z = z p = 8) is sufficiently intensive to affect the IWs propagation.

Damping of internal waves by turbulence
Consider now the effect of the turbulent layer on the internal waves propagating in the pycnocline.In this case, DNS was performed with both the turbulent forcing and the initial conditions (10)-( 12) at t = 0 with IW amplitude W 0 = 0.05 and different wavelengths.Figures 4 and 5 present the results obtained in DNS for λ = 4.The results obtained for wavelengths λ = 5, 8 and 10 are qualitatively similar to those in Figs. 4 and 5.
Figure 4a shows instantaneous distributions of the y component of vorticity, = ∂ z U x − ∂ x U z , and density (ρ + ρ ref (z)) obtained in DNS for IW wavelength λ = 4 at times t = 100 and t = 300.The figure shows that the amplitude of the isopycnal displacement in IW is reduced under the effect of turbulence and becomes quite substantial at late times.Figure 4b shows vertical profiles of the velocity fluctuation and mean density obtained in the same run.The figure also shows the profile of the velocity fluctuation induced by IW in the absence of turbulence (in blue color).As compared to the case with no initially excited IW (cf.Fig. 3b), we observe the effect of the IW on the u profiles, especially in the vicinity of the pycnocline level.Figure 4b shows that for z < 7 the velocity fluctuation is entirely due to IW.
As it was stated above, we neglect the diffusion effect on N ref < ρ > +ρ ref (z), is allowed to vary due to both viscous and turbulent diffusion effects.However, this change is found to be negligible in the present study (cf.Fig. 4b), obviously due to comparatively low Froude number (unity) (i.e. the turbulence is considered to be sufficiently weak).Of course, the turbulent mixing effect would be of greater importance if a stronger turbulence is considered.However, we chose the flow parameters so as to minimize the effect of turbulent mixing and molecular diffusion on the pycnocline and, at the same time, retain the effect of turbulence on IW propagation.Figure 5 compares temporal development of the dimensionless density ρ in the central point x = 20, y = 0, z = 8 (i.e. in the center of the pycnocline) obtained in DNS with the turbulent forcing on, and the results of DNS with the same initial conditions for IW but without turbulence forcing.The figure also shows the instantaneous density difference ρ in the same point obtained in DNS with and without turbulence forcing.The figure shows that turbulence causes a considerable damping of IWs.Note that since in the pycnocline center dρ ref (z)/dz = −1 (cf.Eq. 7), and the density disturbance is small as compared to the reference density, ρ ρ ref (z), the displacement of the isopycnal surface at the pycnocline center (with density ρ ref (z) + ρ = 1.5) from the undisturbed position is equal to (−ρ).
Figure 5 also shows that the damping effect due to turbulence is much stronger than the damping caused by viscosity.By the time t = 300, the amplitude of IW affected by turbulence becomes almost four times less as compared to the initial amplitude (decreasing from 0.8 at t = 0 to about 0.2 at t = 300).On the other hand, the IW amplitude decreases by about 30 % due to viscous damping, when no turbulence forcing is applied (decreasing to about 0.5 at t = 300).The results shown in Fig. 5 can be used to estimate the IW damping rate caused by turbulence.Since the perturbations of the reference state are small, it can be assumed that the decrease of the IW amplitude under the effect of turbulence in terms of the density deviation is described by where ρ t (t) is the amplitude of the instantaneous density deviation in the IW affected by turbulence and ρ nt (t) is the amplitude of the instantaneous density deviation in the IW in the absence of turbulence, and γ is the damping rate due to turbulence.Then, the variation of the IW amplitude is For γ t < 1 Eq.( 18) can be approximately rewritten as Thus, the density difference grows linearly at earlier times, where γ t < 1, and saturates at sufficiently large times, for γ t 1.The damping rate can be evaluated from the data for the density difference ρ in Fig. 5 by deducing the dependence of the IW amplitude on time, both without and with the forcing, and using function ( 19) as the best fit.
Figure 6 shows the value of γ obtained from DNS data in Fig. 5 within the time interval 100 < t < 300, and compares it with the theoretical estimate (Ostrovsky et al., 1996): where N = N ref (z) is the buoyancy frequency given by Eq. ( 5),  is the turbulent viscosity where b = 3u 2 /2 the turbulence kinetic energy.Here ν t is evaluated for the turbulent fluctuation profile u shown in Fig. 3b and the turbulence length scale l t = 0.4; W (z) is the IW vertical velocity distribution in Fig. 2a, and c = ω/k is the IW phase velocity evaluated for given wavenumber, k, and frequency, ω.The empirical factors in the turbulent exchange coefficients defined by Ostrovsky et al. (1996) are Figure 6 shows a good qualitative and quantitative agreement with the DNS results.The error bars at DNS points in Fig. 6 are caused by variation of the amplitude of ρ(t) from the linear increase (cf.Fig. 5 for IW wavelength λ = 4; similar variations are observed for λ = 5, 8 and 10).These variations are due to distortion of the IW by turbulent eddies and increase as the forcing amplitude F 0 increases (as discussed later in Sect.6).
As already mentioned in the introduction, in the theory developed by Ostrovsky and Zaborskikh (1996) it was assumed that the IW is sufficiently weak as compared to the turbulence fluctuation amplitude to allow linearization of the semiempirical equations.This assumption can be considered as applicable in the case of forcing amplitude F 0 = 3 discussed above.Indeed, the average IW velocity amplitude (which can be evaluated as r.m.s.(or dispersion) of Eqs. ( 10)-( 12) by averaging over x with W 0 = 0.05) is about 0.03, and can be regarded sufficiently small as compared to the turbulent velocity amplitude (u ≈ 0.06, cf.Fig. 4b).As shown below, a stronger IW (or weaker turbulence) can violate the applicability of the estimate, as in Eq. ( 20).

Dependence of the damping rate on turbulence forcing amplitude
In order to evaluate the applicability limits of the theoretical formula in Eq. ( 20), we performed DNS with different forcing amplitudes (F 0 = 1, 2, 2.5, 3.5 and 4) and IW wavelengths λ = 4 and λ = 8.The results are presented in Figs.7-9.
Figure 7 shows vertical mean profiles of the fluctuation velocity u and mean density, < ρ + ρ ref (z) >, obtained in DNS with forcing amplitudes F 0 = 2 and 4 and initially excited IWs.The figure also shows the profiles of the velocity fluctuations induced by IW in the absence of turbulence (in blue line).The figure shows that the turbulence is affected by the IW propagating in the pycnocline, especially when the forcing in Eq. ( 8) is relatively small, F 0 = 2.It is interesting to note that the vertical turbulence profile is split due to the IW action.The absolute maxima of u located at z = 10 are due to the choice of parameter z f in the forcing function of Eq. ( 7), and qualitatively analogous to that in Fig. 4b   quantitatively (u ≈ 0.04 for F 0 = 2 and u ≈ 0.08 for F 0 = 4), so that the increase in F 0 leads to a larger turbulence intensity, as expected.The maxima of u ≈ 0.02 at the pycnocline horizon (at z = 8) in both cases F 0 = 2 and 4 are due to IWs action.The figure shows also that at z < 7 velocity fluctuations are entirely due to IW, and in the case F 0 = 4u is reduced as compared to the case with no turbulent forcing.This reduction is due to the damping of IW by turbulence.
Note that under the conditions shown in Fig. 7 in the case F 0 = 2 (a weak forcing), the perturbations in the internal wave and the turbulence have comparable intensities and therefore they affect each other.Indeed, the average velocity amplitude in this case is of about 0.03 which is rather close to the turbulence velocity amplitude (about 0.04).Therefore, in this case the theoretical estimate in Eq. ( 19) is not applicable since it was derived from the hydrodynamic Reynolds equations linearized with respect to the IW.
Note also that the of the second maximum of u is in qualitative agreement with the results of Ostrovsky et al. (1989) where amplification of turbulence under the action of a stronger internal wave was observed.This provides an incentive for a future study of the phenomena of the amplification of small-scale turbulence by IWs.  the instantaneous difference ρ of the densities obtained in the same point in runs with and without turbulence forcing.These data are qualitatively analogous to those shown in Fig. 5. Using these data we evaluated the IW damping rate and compared it with the results of calculation according to the theoretical solution (20). Figure 9 compares the IW damping rate obtained from DNS with the theoretical estimate for different forcing amplitudes.
Figure 9 shows that in the case of a sufficiently weak turbulence (low forcing amplitude, F 0 < 3) and the IW velocity amplitude on the order of the turbulent velocity amplitude, the semi-empirical theory overpredicts the damping rate.As F 0 increases from 1 to 3, and the turbulent velocity amplitude becomes sufficiently large as compared to the IW velocity amplitude, the agreement between the DNS results and the theoretical estimate in Eq. ( 20) improves.
It should be pointed out, however, that as the amplitude F 0 further increases from 3 to 4, the accuracy of the theoretical estimate is again reduced.This reduction can be attributed to the influence of strong turbulent eddies on the pycnocline, when the turbulent mixing changes the pycnocline and, consequently, the IW mode structure.On the other hand, the error of the evaluation of the damping rate from DNS data (e.g., from data in Fig. 8a and b for IW wavelengths λ = 4 and λ = 8) becomes much larger for F 0 = 4.This error increases due to considerable distortion of the IW by turbulent eddies which become strong enough to cause the pycnocline erosion as F 0 increases.Figure 10 compares the instantaneous density contours in the center of the pycnolcine obtained in DNS at t = 300 for forcing amplitude F 0 = 3 and F 0 = 4.The figure shows that IWs become much less regular under the action of turbulence in the case F 0 = 4.This irregularity is also reflected in the behavior of the density in Fig. 8.This indicates a direct effect of small-scale turbulent fluctuations as wave "scatterers".The figure also shows that the instantaneous density contours sometimes overturn, leading to convective turbulent mixing.These overturnings are visible in the more turbulent region, for level ρ = 1.1 at forcing amplitude F 0 = 4 (where the isopycnal is locally almost vertical).

Discussion
Here we studied the damping of internal gravity waves by small-scale turbulence.IW damping rate was calculated from DNS data and compared with the results of the semiempirical theory.We also investigated applicability limits of the semi-empirical approach in terms of the ratio between the amplitudes of the IW and of the turbulent pulsations.Our results show that in the case where IW amplitude is of the order of turbulent pulsation velocity, the semi-empirical formula obtained for weak IWs overestimates the damping rate by the order of magnitude.However, the theoretical estimate becomes quite accurate in the case where the IW amplitude is sufficiently small as compared to the turbulent fluctuations amplitude.
The parameters of the present DNS study can, in particular, be matched to the flow parameters in the experiment by Ostrovsky et al. (1996).In this experiment, internal waves were generated in the pycnocline by a wavemaker, and smallscale turbulence was induced by an oscillating grid located 0.2 m above the pycnocline.The IW period in the experiment ranged form about 23 s to 45 s and IW wavelength varied from 35 cm to 170 cm.The thermocline (or pycnocline) thickness (half-width) was about L 0 ≈ 0.1 m and maximum buoyancy frequency was about N max ≈ 0.4 rad s −1 .Thus, for the internal wave with a wavelength λ I W = 1m and buoyancy period T IW ≈ 35 s, the dimensional quantities matching those in the experiment can be obtained by choosing the length and timescales of L 0 = 0.1 m and T 0 = 2.5 s.Then, IW with dimensionless wavelength λ = 10 and period T = 13 in DNS corresponds to the IW with dimensional wavelength λ IW = 1 m and period T IW ≈ 33 s which are quite close to the experimental parameters.In DNS, the dimensionless buoyancy frequency at the pycnocline center is N 0 = 1 which in this case corresponds to the dimensional value N max = 0.4 rad s −1 in the experiment.Then, the damping rate in Fig. 6, obtained in DNS, which is ranging from γ ≈ 0.008 for λ = 4 to γ ≈ 0.002 for λ = 10, can be recast in the dimensional variables as γ ≈ 3 × 10 −3 s −1 for λ = 0.4 m and to γ ≈ 8 × 10 −4 s −1 for λ = 1 m which is also close to the experimentally observed values of γ .As mentioned, the present DNS study considers the evolution of internal waves in time for a spatially periodic initial condition.This is a more convenient setting for DNS since it assumes x periodic boundary conditions which are easily implemented numerically.Note that the theoretical estimate (20) for the damping rate Im(ω) was obtained by Ostrovsky and Zaborskikh (1996) and Ostrovsky et al. (1996) for the geometry similar to that used in DNS.At the same time, in the laboratory experiment by Ostrovsky et al. (1996) the spatial wave damping was observed.These data were linked in Ostrovsky et al. (1996) to the theory developed in the same paper by letting Im(k) = Im(ω)/v g , where v g is group velocity of the internal wave.This substitution gives a reasonably good prediction for the measured amplitude attenuation (ratio A 1 /A 0 of IW amplitudes after vs. before the turbulent region).In the present paper we make the comparison which, as mentioned, refers to the same spatially homogeneous configuration, whereas the range of parameters (turbulence intensity) in DNS goes beyond those used in the experiment.
In order to link the IW damping observed in DNS to that observed in the experiment by Ostrovsky et al. (1996), let us first evaluate the time interval in which it would take for the IW to travel the distance of 2 in the dimensionless length units in DNS which corresponds to 20 cm in the experiment.This gives (e.g., for IW with wavelength λ = 4 with the group velocity of about 0.1) the dimensionless time interval t = 20, so that for the damping rate γ ≈ 0.008 observed in DNS we obtain the ratio of IW amplitudes A(t = 20)/A(t = 0) = 1/ exp(γ t) = 1/ exp(0.16)= 0.85.This value is in a good agreement with the IW damping observed in the experiment for the considered dimensional IW wavelength 0.4 m (wave number 0.16 rad cm −1 ).
Note that in the present study the distance between the forcing level, z f , and the pycnocline level, z p , is fixed and chosen to ensure that the effect of turbulent mixing on the pycnocline remained sufficiently small for the considered forcing amplitudes in the range of 1 < F 0 < 4. On the other hand, the effect of turbulence on the IW propagation is sufficiently large to be detected in DNS.It would be of interest to reduce the distance between the turbulence forcing level and the pycnocline.We expect that this would result in an enhancement of turbulence intensity at the pycnocline level and, as a consequence, a larger damping rate of IWs for the same amplitude of forcing F 0 , provided the mixing effects remain negligible.An interesting relevant problem is a possibility of turbulence modification in case the energies of the internal wave and turbulence are comparable as implied in the paper by Ostrovsky et al. (1989).We plan to continue the DNS work in this direction.
Note also that since DNS is performed in dimensionless variables, the dimensional quantities can be obtained using similarity with the corresponding use of velocity, length, and timescales.Thus, the above results can be extrapolated to oceanic pycnoclines by an appropriate scaling.If, for example, density stratification can be approximated by the expression in Eq. ( 6) with the vertical scale of 100 m corresponding to N max ≈ 0.01 rad s −1 (Phillips, 1977), then one can use the scales U 0 = 0.1 m s −1 , L 0 = 10 m, and T 0 = 100 s.Then, the scaled turbulence velocity (of the order of O(1 cm s −1 ) is close to that in the typical near-surface turbulence.For the IW with wavelength λ IW ≈ 100 m (period T IW ≈ 1300 s), the damping rate is estimated to be on the order of γ = O(10 −5 s −1 ) which is in good agreement with the estimates made by Ostrovsky and Zaborskikh (1996).
In the present study, we prescribed the IW amplitude to be sufficiently small (and the wave slope ka < 0.1) when the IW nonlinearity can be neglected.This allows us to compare the DNS results with the corresponding theoretical and experimental results of Ostrovsky et al. (1996).If the IW mode is nonlinear it still does not significantly perturb turbulence, since it is expected that the higher wave harmonics will attenuate faster than the main harmonic.The case when the amplitudes of motion in IW and turbulence are comparable is planned to be considered in the future.

Fig. 1 .Δ
Fig. 1.Schematic of the numerical experiment: x,y,z are the Cartesian coordinates; 0 ρ is the

Fig. 1 .
Fig.1.Schematic of the numerical experiment: x, y, z are the Cartesian coordinates; ρ 0 is the density above the pycnocline; ρ 0 the density jump across the pycnocline; g the gravity acceleration, N m the buoyancy frequency in the pycnocline center; L 0 the pycnocline thickness; z p and z f the locations of the pycnocline center and the forced turbulence layer center, respectively.
Fig.2a.Distribution of the vertical velocity W(z) for wavelengths 4,5,8 and 10 (left) and the dispersion relation ) (k ω (right) for the first IW mode.The wavenumbers corresponding to the selected wavelengths λ = 4, 5, 8 and 10 are shown by symbols.

Fig. 2a .
Fig. 2a.Distribution of the vertical velocity W (z) for wavelengths 4, 5, 8 and 10 (left) and the dispersion relation ω(k) (right) for the first IW mode.The wavenumbers corresponding to the selected wavelengths λ = 4, 5, 8 and 10 are shown by symbols.

Fig. 2b .
Fig. 2b.The instantaneous contours of the density deviation obtained in the central (x,z)-plane at different time moments in DNS with initial condition (3.2) prescribed for IW, propagating from left to right with wavelength λ = 10 (frequency ω = 0.489, phase velocity c ≈ 0.78).The turbulence forcing is off.Density contours are 1.3, 1.5, 1.7.Contour 1.5 marks the position of the pycnocline center.

Fig. 2b .
Fig. 2b.The instantaneous contours of the density deviation obtained in the central (x, z) plane at different time moments in DNS with initial condition (3.2) prescribed for IW, propagating from left to right with wavelength λ = 10 (frequency ω = 0.489, phase velocity c = 0.78).The turbulence forcing is off.Density contours are 1.3, 1.5, 1.7.Contour 1.5 marks the position of the pycnocline center.

Fig
Fig.3a.Turbulence power spectrum obtained in DNS with the forcing on and no initially excited IWs at t = 300 at different z-levels.Dashed line shows the Kolmogorov's spectrum asymptotics.Here and in Figs.4-7 below, the forcing amplitude is

Fig. 3a .
Fig. 3a.Turbulence power spectrum obtained in DNS with the forcing on and no initially excited IWs at t = 300 at different z-levels.Dashed line shows the Kolmogorov's k 3/5 spectrum asymptotics.Here and in Figs.4-7 below, the forcing amplitude is F 0 = 3.

Fig. 3b .
Fig. 3b.Vertical profiles of the velocity fluctuation and mean density obtained in DNS with turbulence forcing on and IW with wavelength λ = 4. Blue, green, magenta and black lines are for velocity fluctuation at times t = 1, 6, 9 and 300.Red line is for density at t = 300.The reference (initial) density profile, ρ ref (z), is shown in dashed line for comparison.
Fig. 4a.Instantaneous distribution of the vorticity y-component Ω (in grey scale) and density contours obtained in DNS at t =100 (top) and t = 300 (bottom) with turbulence forcing on and initially excited IW with wavelength λ = 4. Density contours are 1.3, 1.5, 1.7.

Fig. 4a .
Fig. 4a.Instantaneous distribution of the vorticity y component (in grey scale) and density contours obtained in DNS at t = 100 (top) and t = 300 (bottom) with turbulence forcing on and initially excited IW with wavelength λ = 4. Density contours are 1.3, 1.5, 1.7.

Fig. 4b .
Fig. 4b.Vertical profiles of the velocity fluctuation and mean density (in red line) obtained in DNS with turbulence forcing on and off (in black and blue line, respectively) and IW with wavelength λ = 4 at t = 300.The reference (initial) density profile, ρ ref (z), is shown as a dashed line for comparison.

Fig. 5 .
Fig. 5. Temporal development of the density deviation in center of the pycnocline, at a point with coordinates x = 20, y = 0, z = 8, obtained in DNS for IW wavelength λ = 4, with and without turbulence forcing, f ρ and ρ (in black and grey line, left) and the difference

Fig. 5 .
Fig. 5. Temporal development of the density deviation in center of the pycnocline, at a point with coordinates x = 20, y = 0, z = 8, obtained in DNS for IW wavelength λ = 4, with and without turbulence forcing, ρ f and ρ (in black and grey line, left) and the difference ρ ρ f (right).

Fig. 6 .
Fig.6.Damping rate of IW γ caused by turbulence obtained from the data in Fig. 5 vs. IW wavelength λ and the analytical estimate (20).

Fig. 6 .
Fig. 6.Damping rate of IW γ caused by turbulence obtained from the data in Fig. 5 vs. IW wavelength λ and the analytical estimate (20).

Fig. 7 .
Fig. 7. Vertical profiles of the fluctuation velocity u and mean density /rho (red line) obtained in DNS at t = 300 with turbulence forcing on and off (in black and blue line, respectively) and initially excited IW with wavelengths λ = 4 (left coloumn) and λ = 8 (right coloumn) for the forcing amplitudes F 0 = 2 (top) and F 0 = 4 (bottom).The reference (initial) density profile, ρ ref (z), is shown in dashed line for comparison.

Fig. 8a .
Fig. 8a.Temporal development of the density deviation in center of the pycnocline, at a point with coordinates x = 20, y = 0, z = 8, obtained in DNS for IW wavelength λ = 4, with and without turbulence forcing f ρ and ρ (in black and grey line, left column) and the difference

Fig. 8a .
Fig. 8a.Temporal development of the density deviation in center of the pycnocline, at a point with coordinates x = 20, y = 0, z = 8, obtained in DNS for IW wavelength λ = 4, with and without turbulence forcing ρ f and ρ (in black and grey line, left column) and the difference ρ − ρ f (right column) for forcing amplitudes F 0 = 2 (top) and F 0 = 4 (bottom).