The effect of strong shear on internal solitary-like waves

. Large amplitude internal waves in the ocean propagate in a dynamic, highly variable environment with changes in background current, local depth, and stratification. The Dubreil-Jacotin-Long, or DJL, theory of exact internal solitary waves can account for a background shear, doing so at a cost of algebraic complexity and a lack of a mathematical proof of algorithm convergence. Waves in the presence of shear that is strong enough to preclude theoretical calculations have been reported in 5 observations. We report on high resolution simulations of stratified adjustment in the presence of strong shear currents. We find instances of large amplitude solitary-like waves with recirculating cores in parameter regimes for which DJL theory fails, and of wave types that are completely different in shape from classical internal solitary waves. Both are spontaneously generated from general initial conditions. Some of the waves observed are associated with critical layers, but others exhibit a propagation speed that is very near the background current maximum. As such they are not freely propagating solitary waves, and a DJL 10 theory would not apply. We thus provide a partial reconciliation between observations and theory.

Large amplitude internal waves, often referred to as internal solitary-like waves or ISWs, are a well-studied coherent, nonlinear phenomenon accessible via field measurements, laboratory experiments, and simulations of density stratified fluids.Historically, ISWs were described by approximate, perturbation theories that lead to an equation of Korteweg-de Vries (KdV) type 15 for the temporal-horizontal portion of the object (Helfrich (Talipova et al. (1999); Helfrich and Melville (2006); Talipova et al. (2006)).The vertical component is described by solving a linear ordinary differential eigenvalue problem.This is generally referred to as the weakly nonlinear description, or WNL.While aspects of WNL continue to be actively developed in recent literature (e.g.Talipova et al. (1999); Grimshaw et al. (1997); Grimshaw et al. (1997); Talipova et al. (1999); Caillol and Grimshaw (2012)), ( 2012)) with an excellent survey available in Ostrovsky and Stepanyants (2015), it has been known for some time that WNL is lacking as a quantitative descriptor of ISWs.The KdV theory itself predicts that 20 wave amplitude is bounded above by the onset of breaking, in contrast 20 to observations of wave broadening and the formation of waves with flat crests (i.e., tabletop waves; Lamb and Wan (1998); Rusås and Grue (2002)).Exact ISWs are solutions of the nonlinear elliptic eigenvalue problem referred to as the Dubreil-Jacotin-Long Dubreil-Jacotin-Long (DJL) equation.This equation is formally equivalent to the full stratified Euler equations (in a frame moving with the wave), and thus this theory is referred to as fully nonlinear, or exact.
25 Comparisons of WNL with exact ISWs (Lamb and Yan (1996); Lamb (1998); Stastna and Peltier (2005)) typically yield poor results for the vertical structure, especially for complex stratifications.Higher order WNL theory, with the Gardner or mKdV 1 equations as prominent examples, does a better job of qualitatively matching the type of upper bound on wave amplitude 1 (Grimshaw et al. (1997)), but there is a significant gap between what is observed in simulations and what WNL can describe (Stastna and Peltier (2005)).In the context of ISWs with strong shear, this gap is well illustrated by the study of (Caillol 30 and Grimshaw (2012)), which develops a WNL theory for weakly stratified critical layers.The theory is algebraically quite 30 complex with significant ties to work on Rossby wave critical layers.However, the theory is not compared to simulations or experiments, and thus it is unclear to what extent it applies to a situation with a dynamically evolving wave field.Interestingly it is very similar in style to Maslowe (1980) which considered a closely linked topic, some 30 years earlier.
While field measurements are often taken in a dynamic, complex environment with changes in the stratification structure 35 and a non-trivial background current, much of the understanding of ISWs has been developed based on relatively simple stratifications (typically quasi-two layer, or exponential stratification).The presence of a background shear current complicates the algebra needed to derive theory, whether WNL or the DJL equation.In the case of WNL, the numerical methods are unchanged, but for the DJL equation the iterative algorithm used to obtain solutions requires an ad hoc modification (see Stastna and Lamb (2002) and the Methods section below).In practice, publicly available software for the DJL equation (Dunphy et al. 40 (2011)) can handle many situations, but there is no a priori proof of convergence.It has been shown that the presence of a background shear current can modify the type of upper bound on wave amplitude, and can also change the polarity of exact ISWs (again see Stastna and Lamb (2002) for details).The detailed oriented reader is cautioned that Figure 3 in Stastna and Lamb (2002) is based on ISWs of elevation which have an opposite sign of wave-induced voritcity vorticity to ISWs of depression.For ISWs of depression a background current with positive vorticity tends to decrease the limiting wave amplitude.In his study of 45 ISW energetics in the presence of a background shear current current, Lamb (2010) thus employed a negative background current.The DJL equation has also been extended to model ISWs past breaking, or waves with trapped cores.Such waves are known to form during shoaling (Lamb (2002)), and the question of whether the core is completely trapped (Xu et al. (2016)), quiescent: (Derzho and Grimshaw (1997);Luzzatto-Fegiz and Helfrich (2014); Lamb and Wilkie (2004)), Lamb and Wilkie (2004); Luzzatto-Fegiz and Helfrich (2014)), or found immediately adjacent to the boundary or at depth (He et al. (2019)) have all been the subject of recent studies.Any ISW can be described as a 50 propagating baroclinic vortex, but trapped cores provide a more complex wave-vortex coupling that will be revisited below.In the classical theory of hydrodynamic stability the gradient Richardson number (Ri) is the standard necessary condition; when Ri > 0.25 the flow is linearly stable.However, Ri < 0.25 is not a sufficient condition for instability, instability (Hazel and Hazel (1972)).Large amplitude ISWs, both with and without a background shear current, can induce strong shear currents near the wave crest.This situation has drawn interest over several decades (Bogucki and Garrett (1993); Barad and Fringer (2010); 55 Lamb and Farmer (2011); Xu et al. (2019)) first as a pure conjecture, then with two-dimensional and finally threedimensional situations.It has been found that shear instability within ISWs is quite complex, with onset dependent on the strength of upstream perturbations, the Reynolds number, and duration of time spent below Ri= 0.25 as a perturbation passes through the wave.Three-dimensionalization, at least on scales for which direct numerical simulation (DNS) has been accessible, occurs preferentially at the rear of the wave and in some parameter regimes instability is episodic with bursts and quiescent periods.
60 Walter et al. (2016) measured large amplitude ISWs in northern Monterey Bay, California (USA) as part of a complex interplay between a coastal upwelling front and and local wind forcing.These waves were large amplitude (up to a 10 m maximum isopycnal displacement in water approximately 20 m deep) with strong background shear currents.Indeed the 2 background shear currents were so strong that solutions of the DJL equation were impossible to compute.Some analysis of large ISW-like waves in the presence of strong background shear currents is provided in Stastna and Walter (2014).These authors considered the 65 resonant generation of ISWs by flow of a background current with shear over isolated topography.They found large amplitude wave trains, as well as waves trapped behind the topography with vortex rich cores that were quasitrapped quasi-trapped near the surface (i.e., both the upstream propagating wavetrains and the trapped waves exhibited vortex-rich cores).However, the connection between simulations and field measurements is incomplete since it is unclear if the waves observed in Walter et al. (2016) were resonantly generated.70 In recent work, Zhang et al. (2018) considered the effect of a background current on mode-2 waves.The waves were generated via stratified adjustment, similar to the methodology we follow below.The stratification used was centered at the mid-depth, so that in the absence of a shear current initial perturbations have the dominant part of their energy deposited into the mode-2 wave field.The authors document the manner in which the adjustment process, as well as the resulting leading mode-2 wave and trailing mode-1 tail are modified by the presence of shear.They also compare their results to the observations profile and shear layer in Zhang et al. (2018) are centered at the mid-depth makes a detailed comparison impossible.
The interaction of ISWs with background currents has also been explicitly demonstrated in the context of frontal generation.Bourgault et al. (2016) discussed observations and simulations of a case in the Saguenay fjord, Quebec, Canada.The authors found that convergence near the front was important in numerical generation experiments, though the parameter space they explored 80 explored was fairly limited.In all cases discussed by the the authors, rightward and leftward propagating wavetrains propagating away from the front differ in amplitude, with larger waves observed in stronger shear.
In a related non-ISW context, Lamb and Dunphy (2018) considered the generation of internal waves by tidal flow in the presence of shear.They concentrated on characterizing energy flux, and found profound asymmetries between the waves propagating with the shear and against the shear.Interestingly, whether the downstream or upstream side of the ridge that 85 generated the internal waves dominated the energy flux was found to depend on the slope of the flanks of the ridge.It is difficult to ascertain implications of this study for the problem of ISWs with shear, since these authors only considered a linear stratification.Nevertheless, the study clearly suggests that shear can profoundly influence wave energetics.
In this manuscript, we report on high resolution simulations of stratified adjustment and the resulting internal wave field in the presence of strong shear currents.The structure of the manuscript is organized as follows.The Methods section presents 90 the governing equations, the DJL theory and the design of, and parameter values for, all numerical experiments we present.The Results section presents the major findings, divided into three parts: i) waves modified by shear but without critical layers, ii) waves with critical layers, iii) situations in which the stratification is in the opposite half of the water column from the shear layer.The manuscript concludes with a Discussion section and a brief set of Conclusions.
∂t 100 where u is the velocity, P is the pressure, ρ is the density, and ρ0 is some reference density of the fluid.The physical parameters 100 are the shear are the molecular kinematic viscosity ν (set to equal 1× 10−6 m s −2) and scalar diffusivity κ (set to equal 2× 10−7 m s −2).The unit vector in the vertical direction is denoted by k ̂.The computational domain is rectangular with the x-axis running left to right along the bottom of the domain, so that 0< x < Lx and 0< z < Lz .
Simulations were carried out with the pseudospectral code SPINS (Subich et al. (2013)).The code has been thoroughly 105 validated in a number of different configurations (e.g.shear instabilities, internal wave generation, internal solitary wave propagation) and is available for download through its online manual: https://wiki.math.uwaterloo.ca/fluidswiki/index.php?title=SPINS_User_Guide All simulations reported employed regularly spaced grids.Approximately 30 exploratory simulations were carried out, prior to identifying a parameter space of interest.Table 1 provides a list of cases discussed in this manuscript, along with their key 110 physical and numerical parameters.As a general comment, simulations were much higher resolution than many comparable studies in the literature, and the pseudospectral aspect of SPINS means that the model has very little numerical dissipation.This proved vital for some of the results reported below, and is discussed further in the Discussion section.
All simulations were initialized with a background shear current The total depth was fixed as Lz = 20 m for all simulations, but the horizontal extent of simulations was varied on a case by 115 The total depth was fixed as Lz = 20 m for all simulations, but the horizontal extent of simulations was varied on a case by case basis.The thickness of the shear layer was fixed at dU = 3 m, or 0.15 of the total depth.Umax was varied as described in Table 1.The form of the background current is such that, in the absence of stratification, it is linearly stable (i.e.Fjortoft's criterion is not satisfied).
The initial density field was specified as ρ(x,z) = ρ̄(z− η) with While other values were used in preliminary experiments, all numerical experiments reported below set η0 =−5 m (0.25 of the total depth) and wd = 100 m.
The background density field was specified as ∆ρ 6) 2 4 125 The dimensionless density difference was set as ∆ρ= 0.001 for all simulations, while the pycnocline center and thickness, (z0,d), were varied as shown in Table 1.
For no background current the DJL equation takes the more familiar form ∇2 N2(z− η) η+ η = 0 c2 135 with the strong non-linearity encapsulated in the evaluation of the buoyancy frequency squared profile at its upstream height (i.e. at z−η).While several different techniques for the solution of the DJL equation without a background current are described in the literature, the only publicly available approach to the DJL equation with a background current that the present authors are aware of employs an ad hoc extension of the direct variational method of Turkington et al. (1991), as described in Stastna and Lamb (2002).This consists of an iterative algorithm, which while generally stable for weak currents, leads to 'wandering' 140 when shear is strong.The simulations reported below specifically concentrate on parameter regimes in which the DJL was found not to converge (unless otherwise noted in the discussion).Table 1.Cases considered throughout the manuscript, including case label, and their dimensional parameters.

Results
Sample initial conditions are shown in Figure 1.The horizontal component of velocity is shaded and five isopycnals are superimposed in black.Only a portion of the computational domain is shown.The velocity range is chosen to be symmetric in 5 Figure 1.Typical initial conditions.Horizontal component of velocity (m s−1) shaded (case NO DJL 1 shown as a representative sample), 5 isopycanls in black.
145 order to demonstrate the the background current is always greater or equal to zero.Simulations are reported for combinations of background profiles shown in Figure 1. 2. In all cases, the region of highest shear occurs away from the peak in the buoyancy frequency.For cases NO DJL 2 and CL 2 this is due to the thin pycnocline, and for cases with the sub-label B this is due to the fact that the pycnocline center is below the mid-depth.5 Figure 1.Numerical experiments were designed to compare three qualitative categories of results: 150 1.Those with wavetrains that have a DJL description propagating both with (rightward) and against (leftward) the background shear current; Cases DJL and DJLB.
2. Those with wavetrains that have a DJL description propagating against the background shear current (leftward) but no DJL description for wavetrains propagating with the background shear current (rightward); Cases NO DJL 1, NO DJL 1B, NO DJL 1BF, and NO DJL 2. 155 3. Cases with possible critical layers; Cases CL 1, CL 1B, CL 2.
For Cases DJL and DJLB, the background shear current was not large enough to preclude DJL solutions for either the rightward or leftward propagating wavetrains.Figure 2  In both cases a rank ordered wavetrain is formed, with the leading wave described by solutions of the DJL equation.When 160 the pycnocline is in the upper half of the domain, waves of depression form, and the asymmetry between steeper and taller leftward propagating (against background shear) waves waves, and broader, shorter rightward propagating (with background shear) waves is clearly evident.When the pycnocline is in the lower half of the domain, waves of elevation formed, form, the left-right asymmetry is less pronounced, and the fissioning of a rank ordered wave train takes place more slowly.

6
Figure 2. The asymmetry of waves that emerge from the initial conditions is consistent with both KdV and DJL theories (Stastna and 165 blue curve in Figure 1c. Figure 3 Lamb ( 2002)).In KdV theory the presence of the background current changes the coefficients of nonlinearity and dispersion, while in DJL theory the properties of the solitary waves (such as half width) can be computed from the solution itself.In principle it would be possible to use the inverse scattering theory for the KdV equation to predict the number of solitary waves that emerge from the initial conditions for the leftward and rightward wave trains, but at the time of writing the authors are unaware of a practical implementation of this technique.170 2.2 Solitary-like wavetrains without DJL theory We now turn to the more dynamically interesting cases for which DJL theory proves incomplete.We first consider the case labelled NO DJL 1.The stratification corresponds to the black curve in Figure 1a, 2a, while the background shear current corresponds to the blue curve in Figure 1b.This case yields coherent wavetrains propagating both with and against the background 2b.This case yields coherent wavetrains propagating both with and against the background shear current, but the iterative DJL solver is not able to converge to steady solutions for the rightward propagating waves going 175 s.Panel (a) shows the N2(z) field with the background shear current.This is almost certainly due to the low values of Ri in the shear layer, as indicated by the blue curve in Figure 2c. Figure 4 shows the density field in panel (a), and the vertical component of velocity, or w, in panel (b).The response of the density field at this advanced time (t= 21,600 s) is clear: two wavetrains form, with propagation distance altered by the Doppler shift due to the vertical mean of the background current.A large amplitude wave train propagates against the shear 180 remnants and rightward propagating wavetrains both yield regions current (waves are leftward propagating).Its constituent waves are rank-ordered with no sign of instability.The portion of the domain between 8< x < 12 km is dominated by the remnant of the adjustment process, or the portion of the initial disturbance that is deformed by the background shear current to form an overturning region.The wavetrain propagating with the shear 7 current (waves are rightward propagating) can be seen for x > 14 km.Again the waves appear rank-ordered, but are smaller in amplitude and much wider than the waves that make up the leftward propagating wavetrain moving against the shear current.
Figure 4 (moving against the shear current).
185 Figure 5 revisits the above comment on stability, but at an two earlier point points in the evolution of the wavetrains, namely at t= 7200 3600, and 7200 s.We show Ri using a symmetric color range, whereby the blue regions denote static instability (N2 < 0).It can be seen that remnants of the initial adjustment yield a very large region of overturns, but interestingly some overturning extends to very near the front of the rightward (or with the background shear current) propagating wavetrain.Panel (b) shows Ri(z), saturated over a color range−0.25<Ri< 0.25.It can be seen that the leftward propagating (or against the background shear current) wave train only yields regions with Ri < 0.25 very near to the upper boundary.In contrast, the adjustment as the initial density perturbation is strained by the background current.Interestingly, even at early times, some overturning extends to very near the front of the rightward (or with the background shear current) propagating wavetrain.The leading rightward propagating wave is evident 190 at t= 3600 s, and at t= 7200 s three rank ordered waves have separated from the region of instability.At 7,200 s a Rayleigh Taylor type instability is observed at the bottom of the shear layer (17< z < 19 m) for 8.6< x < 9.1 km.All three rightward propagating waves exhibit overturns (blue regions) in their core, while the adjustment remnants yield a region ofRi < 0.25 over a region spanning several kilometers near the upper boundary.In contrast, the leftward propagating (or against the background shear current) wave train only yields regions with Ri < 0.25 very near to the upper boundary with no overturns observed.
195 There are thus three different types of instability observed in these cases: i) Rayleigh Taylor type instability that occurs beneath the background shear current, ii) a stratified shear instability triggered by the stratified adjustment of the rightward propagating wavetrain, and iii) the generation of strong vortex cores in the leading ISW-like waves of the rightward propagating wave train.
The detailed evolution of the wave-trains, as expressed in the density field, is shown in Figures 5 and 6. 6 and 7. Three times are 200 shown (t= 7,200, 14,400 and 21,600 s), with a domain that extends 1 km in the horizontal, but shifts locations as the wavetrain propagates.If the reader was not told of its presence in Figure 5, 6, they would have no reason to suspect that a background shear current was present as the leftward propagating (against the shear) waves take the form of classical "bell-shaped" solitary waves of both WNL and DJL theories.This is in sharp contrast to rightward propagating (with the shear) waves in Figure 6. 7.Here the region near the upper boundary is clearly active, with short wave and "roll" activity evident over a significant portion of 205 the domain.The leading waves take a very different form from classical solitary waves, with sharp crests more akin to Stokes waves.Figure 7 waves, or cnoidal waves (Rubenstein (1999)).Figure 8 reconsiders the near upper boundary region (top 2 m of the domain) for the wavetrain propagating rightward with the shear current.Forty density isolines that range over ρmin < ρ < ρmin+0.2∆ρρmin + 0.2∆ρ are shown.By concentrating the view on this region, it can be seen just how active the near boundary region is.Not only do the rightward propagating waves exhibit quasi-trapped, recirculating cores that remain active over a long time period, it can be 210 seen that a tumultuous region of overturns and rolls initially nearly reaches the two leading waves, but gradually lags behind.This figure clearly illustrates that the most important feature of internal solitary-like waves is their ability to outrun localized regions of overturns and instability, instability (what we labelled as instability ii) above), before such regions drain significant energy from the wave.In the language of stability theory, nearly all instabilities that are spontaneously generated are convective (as opposed to global) in nature.
However, because it takes some time for the solitary-like waves to fission, there is ample time to trigger a stratified shear instabilities over a wide spatial region.
215 In Table 2, we show estimates of the propagation speed of the leading rightward propagating wave going with the background current.It is immediately evident that all cases discussed above have propagation speeds (cest) that are very near the background current maximum (Umax).Moreover the scaled propagation speed (cest/Umax) tends toward 1 as Umax increases.The DJL theory assumes that streamlines (in a frame moving with the wave) connect to infinity both upstream and downstream of the wave.This is not the case for the NO DJL cases.We hypothesize that this is why the DJL solvers fail to converge in this 220 parameter range.Essentially, the wave trains that forms cannot be fully decoupled from the adjustment region for a long enough time so that a vortex dominated leading wave (or several waves) form.
To summarize, in all 'NO DJL' cases shown in Table 1, rightward propagating (with the background shear current) waves with trapped cores and near-surface regions of overturns and local instabilities were consistently observed.These waves cannot be described by presently available solvers for the DJL equation.In all simulations they were observed to be long lived, with 225 active core regions that could transport mass over a considerable distance in the field.

Wave trains with critical layers
The waves described in the previous section (all 'NO DJL' cases) were found to have propagation speeds that were greater, if only marginally, than the maximum background shear current.When Umax > c > 0 a critical layer may form.In order to get some sense of how robust the waves described above are to the presence of a critical layer we present the results of two cases with Umax = 0.5 m s−1.In Figure 8 230 with U −1 max = 0.5 m s .In Figure 9 we show the vertical velocity and density fields over the top four meters for the wave train propagating rightward with the background current at three different times.The leading and second wave are shown, with both exhibiting a prominent quasi-trapped, recirculating core.Interestingly, while the quasi-trapped core of the leading wave appears relatively quiescent, a region of instability is observed ahead of it it, near the top boundary.
The majority of the rightward propagating wavetrains described above involved trapped cores.In contrast, the cases with 235 DJL theory yield non-breaking waves.In order to explore the effect of stratification in the high shear region we performed a series of experiments with a narrower pycnocline, and report one representative example, case CL 2. We found that the leading waves no longer matched the bell-shaped internal solitary waves of DJL theory, and that the short length scale behavior in the high shear region required a substantial increase in horizontal resolution (∆x= 0.12 m for this case).In Figure 9 11.It can be seen that the leading wave has a substantial vortex that is trapped near the upper boundary.The wave train behind the leading wave is associated with a train of vortices that is gradually descending form from the upper boundary into the main water column.Thus what is referred to as 'the leading wave' above is, to a substantial degree, a deformation of the pycnocline by the leading vortex.2. Estimated wave speeds for select cases.

Wavetrains with unstratified shear layers
The results in Figure 10 255 Figure 12. Figure 12a The results in Figure 11 suggest that stratification in the shear layer region has a profound influence on the type of wave that forms via adjustment.In the absence of a shear layer, a pycnocline located in the bottom half of the water column yields ISWs of elevation.We thus performed a number of numerical experiments with a pycnocline located in the botom half of the water column to see if wave trains of elevation behave similarly to those reported above.The results of our numerical experiments indicate that waves propagating against the background current are far less affected when the pycnocline is well removed from the shear layer.In contrast, ISWs propagating with the shear are almost completely 11 destroyed in all but the weakest currents that were tried.Figure 11 12 shows the vertical component of the velocity field shaded with a saturation with 8 density isolines superimposed.The strength of the background shear current increases from top to bottom.In the weakest shear current case [panel (a)] the trademark updraft-downdraft pattern of a DJL ISW can be seen for both the rightward and leftward propagating wavetrains.At the time shown, the wave train has not fully fissioned into rankordered rank-265 ordered ISWs.The reason for selecting such an early time in the evolution, t= 5400 s, can be seen observed in the bottom two panels.Panels (b) and (c) show U0 = 0.3 and U −1 0 = 0.5 m s , U0 = 0.5 m s −1, respectively.In both cases there is no trace of a coherent rightward propagating wave train, and the perturbation vertical velocity field is dominated by relatively short length scale perturbations.In panel (b) an argument could be made for a weak wavetrain, but unlike Figure 6 7 there is no clear separation between a leading wave form and the trailing waves.In the highest velocity case the remnants of the rightward propagating 270 wavetrain trail the disturbances near the surface by at least 1000 m.Since background shear has been found to possibly affect the polarity of ISWs, we performed numerical experiments with an initial condition with reversed polarity (a depression).The results of polarity reversal are compared for U0 = 0.3 m s−1 in Figure 13. Figure 13a reproduces Figure 11b, while Figure 12b 12b, while Figure 13b shows the polarity reversed case.It can be seen that the leftward propagating wave in panel (b) takes the form of an undular bore, consistent with the predictions of WNL theory.The 275 undular bore is much slower than all ISWs previously simulated.Interestingly, the rightward propagating wave train also yields an undular bore, an undular bore, however for later times (not shown) the undular bore largely fades away to irrelevance, and the disturbances of the high shear region dominate the rightward propagating response.
The cases described in this subsection, are perhaps the clearest demonstration of the difference between the weak shear regime discussed in Stastna and Lamb (2002), for which ISWs are modulated in form by the background shear current, and the 280 strong background shear regime.

Discussion
The numerical experiments reported above illustrate that when a strong background shear current is present the set of possible wave-vortex phenomena is considerably larger than the tidy ISWs described by the DJL equation.Indeed Indeed, essentially none of the phenomena we have simulated could be termed as truly steady.The main coherent feature from our simulations one could 285 expect to see in the field for strong background shear current field, for strong background shear current, are ISWs with a strong vortex core in the near surface, high shear region.Classical theory of ISWs with cores generally considers nearly quiescent cores (Derzho and Grimshaw (1997); Luzzatto-Fegiz and Helfrich (2014); Lamb and Wilkie ( 2004)), while recent theory and simulations that includes background shear currents yields so-called 'subsurface cores' (He et al. ( 2019)).Neither of these matches the strong wave-vortex coupling 14 Figure 9. Detail of vertical velocity (shaded and saturated at ±0.01 m s −1) and density fields in the core region for wavetrains propagating with the shear current critical layer case CL 1. 10 density isolines with −5× 10−4 < ρ <−4× 10−4 shown in black.(a) t= 7,200 s, (b) t= 9,000 s, and (c) t= 10,800 s that was observed above (e.g., Figure 7).8).In cases in which a critical layer is possible, the vortex cores coexist with smaller 290 disturbances that extend ahead of the leading ISW.
The simulations we performed required resolution that was much higher than that typically used in coastal ocean simulations (i.e., multi-kilometer domains with horizontal resolution of less than a meter and vertical resolution less than 5 cm).Even recent process studies on mode-2 waves in the presence of shear using the finite volume MITgcm (Zhang et al. (2018)) used a resolution of 5× 0.5 m (which was more than sufficient for the phenomena these authors set out to model).Moreover, the 295 pseudospectral model used has very low inherent numerical dissipation so that even phenomena with relatively few grid points across them are only weakly damped.The high resolution necessary makes extension to three dimensions a complicated task.However, 3D simulations initialized from a two dimensional simulation, i.e. using the early portion of the simulations described above, should be possible with present computational resources.It is likely that performing a detailed energy analysis, as in Lamb (2010) for ISWs propagating against the background shear current (hence ones with a DJL theory), would require 3D 300 simulations.However, a more serious issue when interpreting the above results is in the large gap between theory inspired simulations and those on the regional scale in coastal oceans.To the authors' knowledge there has been no systematic study of ISW behaviour in marginally resolved situations, especially when parameterizations for mixed layers (e.g., the KPP scheme) are employed.As  Even with the above caveats, the simulations we have conducted pose a very clear question for the theorist: "In deep water in which a pycnocline may be far removed from local shear layers can an ISW be destabilized by a shear that would be thought of as so far from the wave that it would be irrelevant to its evolution?"310 4 Conclusions We have provided one possible explanation for why the measurements in Walter et al. (2016) yielded large amplitude waves with no DJL based description.Namely, that for a portion of parameter space, mode-1 ISWs generated via stratified adjustment in the presence of a background shear current are not free waves, since they have a propagation speed very close to the maximum value of the background current.When some stratification is observed near the surface, these waves take the form 315 of ISWs with strongly recirculating cores.The simulations above were two dimensional, while the observations of Walter et al. (2016) were influenced by rotation, and likely include included some large scale three-dimensional structures (comments on small scale three-dimensionalization were made in the Discussion, above).Moreover, the observations take a different form from the wave trains that are spontaneously generated in our stratified adjustment simulations.This is not surprising, and a set of initial conditions more closely tied to the field observations remain 320 field observations remains a clear goal for future work.A different avenue for future work would be to more closely link our results to the work of Zhang et al. (2018) on mode-2 waves.Mode-2 waves lack a DJL theory due to the possible presence of a trailing tail of mode-1 waves, and the question of whether a particular shear current more strongly affects the leading mode-2 wave or its mode-1 tail remains largely open.
Finally we provided well resolved examples of stratified adjustment with critical layers (cases labelled 'CL').When a significant sig-325 nificant stratification was present near the surface, the waves formed consisted of strong vortex cores, with some evidence of fine scale structure on the leading side of the wave.In the absence of stratification near the surface, a train of vortices spontaneously spontaneously formed near the surface and this led to small amplitude deformations of the underlying pycnocline.Our simulations did not yield evidence for the common theoretical assumption that ISWs are the dominant object, with critical layers providing a slow drain of energy (Caillol and Grimshaw (2012)).This observation should be re-examined on scales for which threedimensional 3D DNS 330 is possible.
an incompressible fluid in the absence of rotation that obeys the Boussinesq and rigid lid approximations.The stratified Navier Stokes equations in the absence of rotation read ∂u 1 +u

Figure 6 .
Figure 6.Density fields for wavetrains propagating against the shear current for the case NO DJL 1.(a) t= 7,200 s, (b) t= 14,400 s, and (c) t= 21,600 s.

Figure 6 .Figure 7 .
Figure 6.260 when the pycnocline is well removed from the shear layer.In contrast, ISWs propagating with the shear are almost completely

15Figure 10 .
Figure 10.(a) Detail of density fields for wavetrains propagating with the shear current critical layer case CL 2 (the case with a narrow density profile) t= 4500 s, (b) the vertical velocity component in the top 2 m of the domain in the region between the thick black lines in