Observations

The essence of our echo-sounder, current and density measurements is summarized by two repeated transect lines carried out in a region of the Saguenay Fjord32 (Fig. 2) known to host internal solitary waves of unknown origin33. These data were acquired ∼1 h before the time of low water (20:30 UTC on 10 July 2015 at Tadoussac) during a neap ebb tide of range Δη=2.9 m. The first transect (Fig. 2a–c) shows a prominent solitary wave of amplitude a=7.4 m and maximum vertical velocities w s =±0.2 m s−1 propagating northward along a pycnocline of buoyancy frequency N pyc =0.25 s−1 located at depth h=7.6 m in the undisturbed state ahead. The total depth H along this transect varies between 100 and 150 m (not shown). Idealized numerical simulations (presented below) indicate that such a solitary wave in this stratified and sheared background environment has horizontal lengthscale λ=50 m and induces a perturbation mechanical energy E′=1.9 MJ m−1 (see Methods for details). This wave is therefore highly nonlinear3 (a/h=1), long relative to the surface layer thickness (λ/h=7) but fairly short relative to the bottom layer thickness (H−h) that varies over the transect in the range 92–142 m such that 0.4<λ/(H−h)<0.5. The fact that the wave is short relative to the bottom layer thickness will become important below for interpreting these observations. This wave was approximately the fourth wave seen in this area, and the previous ones seen had greater amplitude, up to a=10 m (Fig. 3).

Figure 2: Field observations. (a,d) Echograms, (b,e) northward velocity v (ms−1) and (c,f) upward velocity w (ms−1). The black contours on b and c are isopycnals 2 kg m−3 apart based on sorted density profiles. The ‘Z’ symbols on a mark the centre of detected density overturns and their size is proportional to the logarithm of the available potential energy of the fluctuation P found in the range 4 × 10−6−3 × 10−4 W kg−1. The ‘+’ symbols between a and b mark the positions of the 20 temperature–salinity profiles collected. The inset highlights Kelvin–Helmholtz billows within the dashed box. These transects correspond to the grey track presented in Fig. 1. Full size image

Figure 3: Echogram of a large internal solitary wave. Echogram of the largest internal solitary wave (a=10 m) seen during this sampling. This northward propagating wave was observed ahead of the intrusion, roughly at the same position (centred around y′=0) as the wave on Fig. 2a seen ∼45 min later. Full size image

The wave of Fig. 2a is followed by what can be interpreted from the echogram as being an intermediate turbulent layer, centred ∼10 m depth and being roughly 20 m thick, populated with Kelvin–Helmholtz billows near the top edge. Eleven density overturns were detected (see Methods for details) throughout this transect (Fig. 2a), of which seven were within the billowing layer behind the solitary wave. These seven overturns are characterized by mean Thorpe scales =0.92 m (s.d.=0.36 m). Assuming a 1:1 ratio between Thorpe and Ozmidov scales34 gives an order of magnitude for the mean dissipation rate within the overturns, consistent with intense ocean turbulent conditions17 and supporting the echogram interpretation that this layer is turbulently mixing.

These observations may suggest that the turbulent layer was produced by the passage of the wave in a manner comparable to that seen in idealized numerical simulations of unstable solitary waves where turbulent fluid is left behind as the breaking wave keeps moving forward and restabilizes. This can happen for example when the underlying bottom topography shallows sufficiently to cause wave breaking35,36. However, this hypothesis is eliminated here because the wave remains short along this transect relative to the bottom layer thickness (0.4<λ/(H−h)<0.5) such that it is virtually unaffected by the underlying changing bottom topography. This is supported by a numerical simulation we have carried out (not shown) in which the propagation of a comparable wave over the changing bottom topography was simulated. As expected, the result showed negligible change in wave shape and behaviour relative to a flat-bottom simulation.

Another possibility to explain the overturning and billowing fluid seen behind the wave could be that it was produced by shear instability that had previously developed within the solitary wave37. Qualitatively, parallels could well be made between the observations of Fig. 2a and the numerical results presented in Lamb and Farmer37 (their Figs 13,14 and 18). However, there are reasons to also reject this hypothesis. First, the vorticity seen in the observed billowing layer is opposite to that expected from shear instability that would arise from such a northward propagating internal solitary wave (Fig. 4). For example, the billow centred around y′=−275 m and at 6 m depth rotates anti-clockwise, while billows arising from wave-induced shear instability would be expected to rotate clockwise37. Another reason to reject this hypothesis is that the reduced stratification within the wave is N2−S2/4>0 (equivalent to N2/S2>1/4)38 such that there is no indication that this wave may have been dynamically unstable.

Figure 4: Shear and reduced stratification. Details of the birth of an internal solitary wave from the head of the interfacial gravity current. (a) Echogram; (b) the logarithm of the shear squared S2 (colour) with contoured isopycnals (black) 2 kg m−3 apart; and (c) the reduced stratification N2−S2/4 (colour) with contoured isopycnals (black). Full size image

An alternative and more plausible interpretation is that the Kelvin–Helmholtz billows arose from the vertical shear associated with an intrusive layer that lies just above the pycnocline and characterized by northward velocity of up to v=0.5 m s−1 (Fig. 2b,e). This interpretation is supported by the fact that the reduced stratification averaged over the length of the intrusion is negative near the top edge of the sheared interface (Fig. 5). The Reynolds number Re=5 × 106 is around three orders of magnitude higher than typical laboratory settings of intrusive gravity currents in two-layer fluids22. This near-surface intrusion differs from deeper intrusions found in other fjords39,40,41,42 and in the Saguenay32 by the fact that it intrudes through the sharp pycnocline where internal solitary waves are most susceptible to being excited.

Figure 5: Background conditions. Observed (solid) and idealized (dashed) background currents (a) and density (b) profiles on the north (black) and south (grey) side of the front. (c) Reduced stratification profile spatially averaged within the intrusion in the [−600, −100 m] interval (Fig. 2). Full size image

The observations suggest that this dynamically unstable intrusion was forced by a small-scale convergent front found at the southernmost end of the transect (y′=−680 m; Figs 2 and 4). This front is characterized by steep isopycnals angled at 25°, strong backscatter intensity indicative of bubble entrainment typical of such fronts43, a large downward vertical velocity w front =−0.2 m s−1 within the first 10 m or so of the water column and of width L≈20 m, a surface density jump Δρ=4 kg m−3 and Rossby number R O =v/(fL)=2 × 102≫1, where f=1.08 × 10−4 s−1 is the Coriolis parameter and v≈0.5 m s−1. This front falls into the type of submesoscale fronts that are ubiquitous in the coastal ocean44,45.

By the second transect (Fig. 2d–f), the leading solitary wave previously seen at 19:20:52 UTC had detached from the head of the intrusion and was captured again 635 m northward at 19:45:12 with an amplitude a=4.4 m. This wave had travelled against the surface current at an average ground speed c=0.4 m s−1. The reason for its reduced amplitude may be due to radial spreading loss as there is no evidence that this wave had lost amplitude due to wave breaking. Evidence of internal solitary wave radial spreading can be seen in this environment from shore-based cameras (for example, Fig. 1). A few other internal solitary waves of smaller amplitudes were also seen trailing behind, as if these waves had somehow been ejected from the head of the intrusion (Fig. 2d), an hypothesis that will be supported below with numerical simulations. During this second transect, the intrusion was still clearly defined and more Kelvin–Helmholtz billows were visible near the top of its sheared interface (Fig. 2d, inset).

Shore-based photography suggests that the phenomenon was phase-locked with the semi-diurnal tide as surface patterns indicative of the concomitant and proximate existence of fronts and internal solitary waves were seen on the 6, 7, 9 and 10 of July in the same area and around the same tidal phase, that is, an hour or so before the time of low water at Tadoussac. July 8 was too windy to capture any interpretable surface patterns. The clearest camera and underwater evidence of this process in action was captured on the 6 (Fig. 1). Note that at the time these data were acquired, we were not aware that this process was taking place. It is by chance that we happened to carry out an echo-sounder transect at that place and time, and it is only during the subsequent data analysis phase of this research that we realized we had sampled something unusual. This also explains why the echo-sounder transect does not extend across the front on that day. At the time of sampling, we were not preoccupied by fronts. Fronts and internal solitary waves will however be the focus of future field experiments.

Numerical simulations

Next, we carried out numerical simulations to determine whether this unusual series of concomitant phenomena (Kelvin–Helmholtz billows, internal solitary waves, intrusion and front) seen on Fig. 2 was simply coincidental or whether they were tied together by a single process, possibly related to intrusive gravity currents. Given the number of different nonlinear and nonhydrostatic processes involved, we based our analysis on the two-dimensional incompressible Euler equations. These are the simplest equations that can explicitly deal with all salient features of this problem (see Methods for details).

A first numerical simulation was carried out in a 150 m deep flat-bottom open channel initialized with a three-layer idealization of the observed front (Figs 5 and 6a,d,g; Methods for details). Importantly for the rest of this demonstration, the ambient is not initially still but is characterized by a convergent flow at the front, and the background environment is sheared, a complex natural situation not examined before in laboratory experiments. This configuration is equivalent to a partial-depth intrusion because the intrusive fluid of density ρ 1 does not initially occupy the entire water depth31. For partial-depth and initially motionless intrusions, the degree of non-equilibrium can be characterized by the following non-dimensional density ratio31

Figure 6: Numerical results. Numerical simulations of Kelvin–Hemlholtz billows and internal solitary waves generated by a frontally forced intrusion. (a–c) Density ρ, (d–f) northward current v and (g–i) upward current w with density contours superimposed as grey lines. The top row (a,d,g) represents the model initialization (t=0) with parameters ρ 0 , ρ 1 , ρ 2 , v 1 , v 2 , v 01 and v 02 as defined in the Methods section and given in Table 1. Only the top 50 m are shown but the model domain is 150 m deep. Full size image

where

Equilibrium situations correspond to ξ=0. For this simulation, ξ=−0.18 (Table 1, run #1). The situation can also be characterized by the ratio of layer thicknesses of the ambient22, that is, Δ=[h 2 −(H−h 2 )]/H=−0.89. Both interfaces on each side of the front are initially dynamically stable with minimum interfacial shear Richardson number Ri 1 =1.2 on the south side (that is, y<0) and Ri 2 =1.6 on the north side (y>0). Finally, both regions are also initially subcritical with corresponding composite Froude numbers Fr 1 =0.26 and Fr 2 =0.73 less than unity (Table 1, run #1).

Table 1 Parameter values (in meter–kilogram–second) used for each simulations carried out. Full size table

The result of this control run shows remarkable similarities with the field observations (Fig. 6). Such a convergent front forces a dynamically unstable intrusion comparable in thickness and velocity to the observations and characterized by roughly 5 m tall Kelvin–Helmholtz billows developing on the upper interface. Internal solitary waves are generated in both directions. A few southward propagating ones with amplitudes a<4 m are generated near the front during the first few minutes. The leading wave induces an energy perturbation E′=2.0 MJ m−1. These southward propagating waves may be similar to those generated by river plumes25 with the difference that, here, wave generation does not require the flow to be initially supercritical. We cannot confirm whether such smaller amplitude southward propagating waves were also emitted in the Saguenay Fjord as sampling did not extend far enough southward of the front.

The most remarkable solitary waves are generated at the head of the intrusion and propagate northward. The waves are rank-ordered with the leading wave having, at t=1,100 s and y=901 m, amplitude a 1 =14.1 m, horizontal lengthscale λ 1 =70 m, phase speed c 1 =0.91 m s−1 and energy perturbation E′=7.8 MJ m−1. The head of the intrusion continuously emits waves northward with one wave being emitted roughly every six buoyancy periods ( =2π/N 2 =25 s). This provides a continuous northward depth-integrated perturbation energy flux F′ that reaches an instantaneous value of max(F′)=79 kW m−1 for the leading wave and then decreases with time (Fig. 7). The total wavetrain perturbation energy induced during about 30 min , which is a lower bound estimate of the events duration observed in the field, is E′=37.4 MJ m−1.

Figure 7: Energy fluxes. Simulated southward (negative values) and northward (positive values) depth-integrated perturbation energy fluxes F′ at (dashed) y=−1 km and (solid) y=1 km for different parameter values. The blue curves represent the control run based on an idealization of the field observations (see Fig. 6 for a more explicit representation of this run). Other curves represent sensitivity tests carried out using the parameters listed in Table 1. Blue, control run #1; black, run #2; magenta, run #3; green, run #4. Full size image

Three other simulations were carried out to explore the sensitivity of the resulting internal solitary waves to whether the initial condition is in equilibrium or not and/or the stillness of the ambient (Table 1). The first of these tests examines the situation where the initial density field is identical as the control run discussed above (that is, ξ=−0.18) but where the ambient is still (run #2 in Table 1). This corresponds to a non-equilibrium situation where, according to laboratory experiments, large-amplitude solitary waves leading the intrusion head are expected to be generated. This is indeed what the numerical results show but the northward propagating waves generated are at least an order of magnitude less energetic than the control run (Fig. 7, black). The second of these tests examines the still ambient with the initial density field in equilibrium (that is, ξ=0, run #3 in Table 1). Virtually, no waves are generated, in agreement with laboratory experiments (Fig. 7, magenta). Finally, we examine the case where the initial density field is in equilibrium but the ambient fluid is convergent and sheared (that is, ξ=0, run #4 in Table 1). Interestingly, this equilibrium situation generates very large-amplitude internal solitary waves (Fig. 7, green), much larger than the non-equilbrium and still simulation (Fig. 7, black), and almost as large and as energetic as the non-equilibrium control run (Fig. 7, blue). These results show that intruding gravity currents, when forced, can generate large-amplitude internal solitary waves even when the initial state at rest is in equilibrium.