1 Introduction

Internal waves are ubiquitous in a stratified ocean, where they are generated by tidal flows over the bottom topography, which features continental shelf slopes, subsurface ridges, and seamounts. They have received much attention in recent decades, as they have been recognized to play an important role in mixing the deep ocean, consequently having huge impacts on large-scale circulation and global climate (Munk and Wunsch, 1998; Garrett and Kunze, 2007). The literature features numerous studies that have theoretically investigated the generation of internal waves (Llewellyn Smith and Young, 2002, 2003; Petrelis et al., 2006; Echeverri and Peacock, 2010; Mathur et al., 2016), numerical simulations (Niwa and Hibiya, 2004; Carter et al., 2008; Kang and Fringer, 2012; Buijsman et al., 2012, 2014), and in situ observations (Carter and Gregg 2006; Nash et al., 2006; Alford et al., 2011; Vic et al., 2018). However, these different methods have many restrictions. First of all, theoretical analysis is subject to different assumptions, e.g., linear and inviscid dynamics. Numerical simulations of large-scale systems inevitably have to parameterize mixing and are subject to other numerical uncertainties. Moreover, in situ observations can only be done at sparse spatial locations and are usually insufficient for obtaining an overall picture of the wave generation process.

In recent years, well-designed laboratory experiments have been performed, which provide another approach for quantifying the internal wave generation process. Gostiaux and Dauxois (2007) and Wang et al. (2012) studied the generation of internal wave beams over supercritical shelf slopes (i.e., topographic slope s is greater than the internal wave beam slope α), and demonstrated the presence of a third wave beam (transverse to the topography) at an abrupt shelf break. Zhang et al. (2008) and Jia et al. (2014) investigated the generation of internal waves over a critical (s = α) continental shelf slope and a triangular ridge with both critical slopes, respectively. They found that strong boundary currents were resonantly generated over the critical slope. Wang et al. (2015) showed that the resonant generation of internal waves can also occur for a subcritical (s < α) triangular ridge under particular conditions, which is a bit unexpected. Peacock et al. (2008) studied the generation of internal waves for both supercritical and subcritical ridges, and their results compared well to theoretical predictions. The aforementioned laboratory experiments were performed in the laminar, or small tidal excursion (see Eq. (3) below), regime. Echeverri et al. (2009) also studied the situation for both supercritical and subcritical ridges, but the tidal excursion was extended from small (1% of the topographic width) to relatively large (15% of the topographic width) values. The results showed that, for the large excursion case, in spite of the nonlinear phenomena that break down the wave beam structure and generate harmonics and inter- harmonics, linear theory can still provide a reasonable estimate of the energy conversion rate.

For the lab-scale models, direct numerical simulations can be performed, which avoid mixing parameterization and can reduce many other numerical uncertainties (Echeverri et al., 2009; Gayen and Sarkar, 2010; Jalali et al., 2014; Chen et al., 2017). Gayen and Sarkar (2010) extended the lab work of Zhang et al. (2008) in the laminar to the turbulent flow regime, and found that convective and shear instabilities caused the transition to turbulence. Jalali et al. (2014) investigated the effects of tidal excursion, which was systematically increased from small to a value equal to the topographic width, and also demonstrated the characteristics of wave energetics and turbulence. Chen et al. (2017) studied the generation of internal waves in the laminar flow regime, and demonstrated that the wave beam direction was determined by the geometric constraint of the topography.

As mentioned above, on a critical continental shelf slope, the transition of the strong boundary flow to turbulence was identified by Gayen and Sarkar (2010) to be caused by convective and shear instabilities. However, is there any other flow instability that can contribute to the transition to turbulence? For a critical triangular ridge, the present direct numerical simulations show that the internal wave beams emitted from the ridge top become unstable by subharmonic instability, which is an important nonlinear resonant interaction mechanism that is responsible for transferring energy from high to low frequency waves in the ocean interior (Onuki and Hibiya, 2015). Baroclinic energy conversion and radiation rates up to the appearance of turbulence and the relation between energy and energy flux are also quantified.

The rest of this paper is organized into sections. Section 2 describes the physical model studied and the numerical setups. Section 3 presents results and discussion, which include spatial and temporal evolutions of the flow field, the relation between baroclinic energy and energy flux, and the baroclinic energy conversion and radiation rates up to the appearance of turbulence. Finally, the conclusions are summarized in Section 4.

2 Physical Model and Numerical Setups

The two-dimensional (2D) model studied extended from x = -2.5 m to 2.5 m, with an upper surface at z = 0 m and a maximum water depth of H = 0.25 m. The model topography is an isosceles triangle with its three vertices located at (x, z)= (0, -0.14), (-0.11, -0.25), and (0.11, -0.25) m, respectively (see Fig. 1, in which only the center part of the domain is shown). Thus, the triangular model ridge has two flanks with the same slope s = 1 (the ridge height 0.11 m is identical to the horizontal extent of the slope). The filled water was linearly stratified, with the surface and flat bottom density being 1005 kg m^-3 and 1025 kg m^-3, respectively. This results in a constant buoyancy frequency of

$N = \sqrt { - \frac{g}{{{\rho _0}}}\frac{{\partial \rho }}{{\partial z}}} = 0.88{\rm{ rad }}\;{{\rm{s}}^{ - 1}}, $

(1)

in which g is the gravitational acceleration, and a reference density ρ₀ = 1015 kg m^-3 is used. At the left and right open boundaries (x = ± 2.5 m), the barotropic velocity in the form of u =u₀sin(σt) was imposed, in which u₀ and σ = 2π/T represented the velocity amplitude and the frequency, respectively. The oscillation period was set to T = 10 s, resulting in a frequency of σ= 0.628 rad s^-1. The internal wave beam slope is then given by

$\alpha = \sqrt {\frac{{{\sigma ^2}}}{{{N^2} - {\sigma ^2}}}} = 1, $

(2)

in which the rotational effect is neglected. Given that α = s, the triangular model ridge has two slopes that are both critical. We increased the barotropic velocity amplitude u₀ from 0.001 to 0.005 m s^-1 by a step of 0.001 m s^-1. Thus, a total of five cases of simulation were performed. This corresponds to a tidal excursion parameter given by

${\varepsilon} = \frac{{{u_0}}}{{\sigma L}}, $

(3)

in the range of 0.014≤ ε ≤ 0.072, in which L= 0.11 m is the horizontal extent of the ridge slope. Even though ε remains much smaller than one, we see that turbulence has occurred in the range of ε being considered.

The above physical model is similar to that used in Wang et al. (2012), Jia et al. (2014), Wang et al. (2015), and Chen et al. (2017). However, in these latter references, u₀ was so small that the flow was always in the laminar regime, corresponding to the generation of linear and well-organized internal wave beam patterns. Our present work goes beyond this laminar to the turbulent regime, which corresponds to the generation of strongly twisted and disorganized beam patterns (see Fig. 1).

Direct numerical simulations were performed using an in-house numerical code (Chen et al., 2013, 2017) that solved the full Navier-Stokes equations. The kinematic viscosity and mass diffusivity were set to be molecular values, namely, ν = 10^-6 m² s^-1 and κ = 10^-7 m² s^-1, respectively. The upper surface was approximated as a rigid lid, where the horizontal velocity was allowed to be stress- free. At the bottom boundary, a no-slip condition was applied. Near the left and right open boundaries, sponge layers were added to avoid wave reflections. The horizontal direction was discretized using 4000 grid cells, with Δx = 0.8 mm between x = ±0.3 m. Then, the spacing was gradually stretched to Δx = 2.5 mm near the two open boundaries. The vertical direction was discretized using 250 grid cells, with Δz = 0.5 mm near the flat bottom and Δz = 0.28 mm near the ridge top. Then, Δz was stretched upwards by a factor of 1.005. The boundary layer thickness was estimated to be $\sqrt {2\nu /\sigma } = 1.78{\rm{mm}}$, in which there were at least 3 grid cells. The time step was Δt = 0.005 s, yielding 2000 steps per tidal period. Starting from an initial motionless state, 20 periods of simulations were performed, and the results from the last ten periods were selected for analysis. For each simulation case, a companion run with homogeneous fluid was also performed to obtain the barotropic flow field, as was done in Wang et al. (2012) and Chen et al. (2017).

In the convergence tests, doubling the spatial and temporal resolutions changed the computed velocities by less than 1% and halving the resolutions changed the velocities by less than 3%.

3 Results and Discussion 3.1 Spatial and Temporal Evolutions

The internal wave beam patterns at the end of simulations (t = 20T), which we obtained by contouring the vorticity $\omega = \partial w{\rm{/}}\partial x - \partial u/\partial z$, are shown in Fig. 1. Only the cases for u₀=0.001, 0.003, and 0.005 m s^-1 are shown. At a small u₀ (Fig. 1(a)), very clear and well-defined beams are emitted from the ridge top and are reflected at the upper and lower boundaries to form classic X-patterns. As both ridge slopes are critical, two downward propagating beams lie exactly on the slopes. When u₀ is increased to 0.003 m s^-1 (Fig. 1(b)), the beams become slightly twisted and blurred. Furthermore, apparent signals (secondary beams) of different frequencies (making different angles with the horizontal direction) are emitted from near the ridge top. When u₀ is further increased to 0.005 m s^-1 (Fig. 1(c)), the flow near the ridge top becomes turbulent. Given that this is the source region where the wave beams are emitted, the primary beams are now strongly twisted and lose their coherent structures. We do not increase u₀ any further as turbulence has occurred.

The time series of vorticity at the three positions labeled A, B, and C in Fig. 1(a) were extracted for the power spectrum analysis in Fig. 2, wherein we normalized the horizontal axes by the forcing frequency σ. For the case with the smallest forcing (black lines), the significant peak only appears at the forcing frequency, indicating that the dynamics is in the linear regime, which is consistent with the very clear and well-defined beam pattern in Fig. 1(a). For larger forcing (red and blue lines), significant subharmonic peaks also appear. The peak at the forcing frequency for the black line in Fig. 2(a) is very small and almost indistinguishable, which is due to the weak linear response and the fact that point A is located away from the primary beams. Moreover, as point A is located in a region where secondary beams pass through, one of the two subharmonic peaks for larger forcings has an intensity comparable to that of the primary peak. For the blue line, these two subharmonics each have a frequency of 0.4σ and 0.6σ, respectively, satisfying a resonant relation of 0.4σ + 0.6σ = σ, thus suggesting that they are generated by a subharmonic instability of the primary beams (e.g., Lamb, 2004). At point B located in the left-upward- going primary beam (Fig. 2(b)), the intensity of the primary peak for the blue line is only slightly larger than that for the red line. This is because, at the largest forcing, the primary beam is strongly twisted and loses its coherent structure. The two subharmonics for the blue line also satisfy a resonant relation of 0.4σ + 0.6σ = σ. At point C located right above the ridge top (Fig. 2(c)), the two subharmonics for the blue line now satisfy a resonant relation of 0.3σ + 0.7σ = σ, which is slightly different from the situations at points A and B. Instabilities in the present work do not produce subharmonics at exactly half the forcing frequency (i.e., 0.5σ + 0.5σ = σ). This is also the situation of the subharmonic instability in Lamb (2004), in which two subharmonics of 0.65σ and 0.37σ are produced near the ridge top. The most apparent feature of Fig. 2(c) is that, for the red and blue lines, the broad bands of higher harmonics with the nonvanishing intensities appear, indicating the appearance of turbulence. This result is consistent with the spatially chaotic structures of the flow field near the ridge top (see Fig. 1).

The vertical dashed lines in Fig. 2 denote the position of the buoyancy frequency N in the spectrum. Given that point C is in the region where turbulence is generated, the broad bands of frequencies higher than N can exist in Fig. 2(c). However, for points A and B away from the ridge top, the power falls off quickly crossing the vertical dashed lines (Figs. 2(a) and 2(b)). This is because the background stratification does not support freely propagating waves with frequencies larger than N.

3.2 Baroclinic Kinetic Energy and Flux

The baroclinic kinetic energy density is calculated as E =ρ(u'²+w'²)/2, in which u' and w' are the horizontal and vertical velocity disturbances obtained by subtracting the corresponding barotropic values (from companion barotropic runs) from the full ones. The time-averaged E for different cases are shown in Fig. 3. Given that both ridge slopes are critical, strong boundary currents are resonantly generated along the slopes, where high energy densities appear. This was also shown in the experiment of Zhang et al. (2008), in which a critical shelf slope model was used. In the present model setup with a triangular ridge, high energy density also occurs above the ridge top, where turbulence appears as u₀ increases. Moreover, from Fig. 3 the two upward going beams from the ridge top are clearly much stronger than the two with the same directions and from the slope bottoms. This is consistent with the experimental results shown in Fig. 2 of Jia et al. (2014), in which the authors demonstrated that the upper beams of velocity disturbances are much stronger than the lower ones. Sections of the energy density in Fig. 3 along the left magenta dashed line in Fig. 1(c) are shown in Fig. 4, where the horizontal axis starts from the point (x = -0.2 m, z = -0.25 m). Clearly, the intensity of the upper beam (near l= 0.2 m) is much stronger than the lower one (near l = 0.05 m). Their maximum intensity ratios are 3.29, 3.35, and 4.23 for u₀ = 0.001, 0.003, and 0.005 m s^-1, respectively.

The horizontal and vertical baroclinic energy fluxes are calculated as F_H= p'u' and F_V= p'w', respectively, where p' is the pressure disturbance obtained by subtracting the corresponding barotropic value (from the companion barotropic run) from the full one. The time-averaged spatial structures of these two fluxes are similar and only the horizontal ones are shown in Fig. 5. The positive and negative values denote the rightward and leftward fluxes, respectively. As can be seen, two strong downward propagating beams on the slopes do not start from the ridge top, but from points (indicated by the magenta arrows in Fig. 5(a)) on the slopes and close to the ridge top. These points are the so-called 'amphidromic points', from which the internal wave beams are emitted to propagate in opposite directions (Gostiaux and Dauxois, 2007). The laboratory experiments of Gostiaux and Dauxois (2007) showed that, for a supercritical shelf slope topography, this amphidromic point is located at the shelf break. Chen et al. (2017) inferred that it is located on the slope and is at some distance to the shelf break when the slope is critical. However, these latter authors did not calculate the energy flux to directly demonstrate the position of this point. Now, for the present model ridge with two critical slopes, this energy flux is calculated, directly showing the presence of one amphidromic point on each slope.

One interesting feature of Fig. 5 is that the in ternal wave beams emitted upwards from the amphidromic points are of much weaker intensity compared with the ones emitted downwards and reflected at the slope bottoms. This is very different from the situation shown in Fig. 3, in which the intensities of the upper wave beams are much stronger than those of the lower ones. The same sections as those in Fig. 4 but for the energy flux fields in Fig. 5 are shown in Fig. 6. The ratios of maximum intensity of the upper beam to that of the lower one are 0.29, 0.37, and 0.71 for u₀ = 0.001, 0.003, and 0.005 m s^-1, respectively (compared with the corresponding ratios of 3.29, 3.35, and 4.23 for the energy density fields). Thus, the strong energy density does not imply a strong energy flux. Dettner et al. (2013) obtained similar results by direct numerical simulations and inferred that the more easily measured kinetic energy density cannot be used as a proxy to characterize the conversion of tidal energy into radiated internal wave power. However, in the infinitely deep ocean model of these latter authors, it was only in the case with supercritical topography that the lower beam of energy flux was stronger than the upper one. In the critical topography case, it was still the upper beam that was stronger than the lower one. This finding is different from the present work with a finite water depth model (the bottom topography height occupies almost one half of the total water depth).

Thus far, the disparity between kinetic energy density and energy flux (in the critical case of the present work or the supercritical case of Dettner et al. (2013)) has not been explained. Fig. 7 shows u' and p' at the end of simulation for u₀ = 0.001. As can be seen, while u' forms classic X-patterns (Fig. 7a), as in Fig. 1(a), it is not the case for p' (Fig. 7b). Due to the blocking effects of the topography, strong pressure disturbances exist near the bottoms of the slopes, thus corrupting the X-pattern. The magnitude of u' in the upper beam is about 1.5 times larger than that in the lower one (similar situation for w', but not shown). Given that the kinetic energy density is calculated based on these velocity disturbances, the upper beams in Fig. 3 are stronger than the lower ones. However, p' near the slope bottoms is more than 5 times stronger than that in the upper beams. As the energy fluxes are calculated by multiplications of pressure and velocity disturbances, the lower beams in Fig. 5 are stronger than the upper ones. This result is consistent with the in situ observations of Nash et al. (2006) on the flanks of Kaena Ridge, Hawaii. In their work, they demonstrated that much of the velocity signal (and hence energy) was in high modes (and thus beams), but the pressure, and hence the energy flux, was concentrated in the low modes.

3.3 Energy Conversion and Radiation

The barotropic-to-baroclinic energy conversion rate is calculated by the following equation (Khatiwala, 2003):

$C = \int_x {p'(x, z = h(x), t)U(x, z = h(x), t)\frac{{{\rm{d}}h}}{{{\rm{d}}x}}{\rm{d}}x}, $

(4)

where h(x) is the bottom profile and U is the horizontal barotropic velocity. Fig. 8(a) shows the time average of this conversion rate and the total energy flux away from the topography (by integrating energy fluxes across the two magenta dashed lines in Fig. 1(c), together with the linear theory predictions of Petrelis et al. (2006). As u₀ increases, the energy conversion rates increase, and it is apparent that the results by our simulations (circle points) agree quite well with the theoretical predictions (square points). At u₀= 0.005 m s^-1, when the departure is the largest, our simulation agrees within 13% of the theory, indicating that the linear theory of Petrelis et al. (2006) is a very good prediction of the energy conversion rate, even when turbulence occurs near the ridge top. The total energy flux away from the topography (triangular points) behaves similarly to the total conversion rate, and the difference between them represents the energy dissipation rate within the region, which is enclosed by the ridge slopes and the two magenta dashed lines in Fig. 1(c). The total fluxes normalized by the conversion rates (circle points) are shown as an inset in Fig. 8(a), which represents the fractions of the converted baroclinic energy that are radiated. As can be seen in the figure, this faction increases first and then becomes saturated for the range of u₀ considered. The corresponding range of the tidal excursion parameter is 0.014 ≤ ε ≤ 0.072. Interestingly, the behavior in this inset is similar to that in Fig. 8 of Chen et al. (2013) for the range of small ε. The fraction of radiated energy in this latter work is obtained for a model of much larger scale. For a small ε and critical topography, the fraction in this latter work is about 80%, which is comparable to the present values in the inset. However, it should be noted that, while direct numerical simulations are used for the present lab-scale model, mixing parameterizations must be used in Chen et al. (2013) for the much larger scale model. Thus, the energy budget calculations in the present work are quite robust for models of different scales.

Finally, given that the linear theory of Petrelis et al. (2006) predicts that the conversion rate is proportional to $u_0^2$, the square roots of the values in Fig. 8(a) are plotted in Fig. 8(b). As can be seen, our calculated values (circle and triangular points) scale almost linearly with u₀, further validating the calculations in the present work.

4 Conclusions

Direct numerical simulations are performed for a labscale model with barotropic tidal flows over an isosceles triangular ridge with two critical slopes. This lab-scale model is identical to that used in the laboratory experiments of Jia et al. (2014). The tidal velocity amplitude u₀ is increased from 0.001 to 0.005 m s^-1, corresponding to a tidal excursion parameter in the range 0.014 ≤ ε ≤ 0.072, which is much less than one. However, the occurrence of turbulence near the ridge top produces broad bands of super-harmonic frequencies, and the part which is higher than the buoyancy frequency cannot propagate far away from the ridge top. Subharmonic instability occurs at different positions, which may contribute to the generation of turbulence, whereas a previous work has identified convective and shear instabilities as the mechanisms that generate turbulence on a critical shelf slope (Gayen and Sarkar, 2010).

Two amphidromic points are identified on the critical slopes, from which energy fluxes are emitted in the opposite directions. The beams of energy fluxes going upwards are much weaker than the one going downwards. However, the reverse is true in terms of the kinetic energy density. Even though the upper beams of velocity disturbances are stronger than the lower ones, strong pressure disturbances exist near the slope bottoms, which contributes to the stronger energy flux in the lower beams. Thus, a strong kinetic energy density does not imply strong energy flux, even when turbulence has occurred.

The tidally averaged energy conversion rate and the energy flux away from the topography are both proportional to $u_0^2$, as predicted by linear theory. For the largest u₀ considered, the conversion rate by our simulation agrees within 13% of the theoretical prediction. This implies that the linear theory can provide a very robust prediction of the conversion rate, even when turbulence occurs, which is a bit surprising. The fraction of radiated baroclinic energy increases first and then becomes saturated as u₀ increases. This finding agrees with the corresponding behavior in a much larger scale system with mixing parameterizations, indicating that the present calculations are robust for models of different scales.

Acknowledgements

This work was jointly supported by the Key Research Program of Frontier Sciences, CAS (No. QYZDJ-SSW- DQC034), the National Natural Science Foundation of China (Nos. 41430964, 41521005, 41776007, 41506005), the Pearl River S & T Nova Program of Guangzhou (No. 201610010012), the Youth Innovation Promotion Association CAS (No. 2018378), and No. ISEE2018PY05 from CAS.