## Abstract

The time-dependent boundary layer induced by a weakly nonlinear solitary internal wave in shallow water is examined through direct numerical simulation. Waves of depression and elevation are both considered. The mean density field corresponds to that typical of the coastal ocean and lakes where the lower fraction of the water column is subject to the stabilizing effect of a diffuse stratification. Sufficient resolution of the “inviscid” dynamics of the boundary layer is ensured through use of a Legendre spectral multidomain discretization scheme in the vertical direction. At higher Reynolds numbers, where the simulations become underresolved, because of restrictions in available computational resources, spectral accuracy and numerical stability at the scales of physical interest are preserved through use of a penalty scheme in the vertical and explicit spectral filtering. Thus, a highly accurate description of the qualitative dynamics of the wave-induced global instability is possible and finescale physical mechanisms critical to the appearance of this instability are not smeared out by the high artificial dissipation inherent in lower-order finite-difference schemes. Results indicate that, for a wave amplitude exceeding a critical value, the global instability occurs in regions near the bottom boundary where the wave induces an adverse pressure gradient. The structure of the associated separation bubble is modified through the establishment of coherent and synchronous dynamics, characterized by elevated levels of bottom shear stress and a periodic shedding of coherent vortex structures. Although details of the vortex shedding depend on the particular wave forcing involved, these vortical structures always ascend high into the water column. All findings suggest that this global instability is a potent mechanism for benthic turbulence, mixing, and possible sediment resuspension in shallow waters, presumably even more intense than the nominal turbulent boundary layer.

## 1. Introduction

### a. Long internal waves in shallow waters

It is now firmly established that packets of long internal waves appear frequently and widely on the continental shelves (Apel 1981; Ostrovsky and Stepanyants 1989; Stanton and Ostrovsky 1998; Bogucki et al. 2005). Their presence in lakes and estuaries is also well documented (Thorpe et al. 1972; Hunkins and Fliegel 1982; Farmer 1978; Saggio and Imberger 1998; Antenucci et al. 2000; Boegman et al. 2003). Furthermore, at least a mechanistic understanding of processes underlying their generation exists in most contexts. Their origin at a number of locations in the oceanic case can be linked to interaction of the barotropic tide with topographic features on and along the shelf break. In such instances a prominent thermoclinic depression forms in the lee of the shelf during the ebb phase of the barotropic tidal cycle. Later, during the flood stage, this depression moves onto the shelf and deforms, giving rise to formation of a long-wave packet as the “initial” depression evolves under the action of weak dispersion and nonlinear steepening (Maxworthy 1979; Maxworthy and Lansing 1984; Lamb 1994). It should be noted that conversion of energy in the baroclinic tide along the shelf break might also give rise to the appearance of long waves on the shelf, but a clear understanding of the mechanism underlying their generation through this type of forcing is still lacking. It is also evident that strong convergence zones along the shelf break, formed by prominent topographic irregularities (e.g., canyons penetrating into the shelf), can lead to localized sources of particularly energetic packets (Bogucki et al. 2005).

The appearance of long internal waves in fjords and estuaries can be related to tidally forced currents that pass over sills or other topographic features (Farmer and Armi 1999). Alternatively, in some circumstances where the temporally varying tidal flow reaches strengths commensurate with the phase speed of a local internal wave mode, a topographic resonance may arise (Grimshaw and Smyth 1986; Redekopp and You 1995). In the case of lakes and reservoirs, the primary forcing enters via the action of a surface wind stress. When a wind stress persists for times approaching a significant fraction of the seiche (pendulum) period for mode-1 internal oscillations in the closed basin, a basin-scale baroclinic pressure gradient oriented along the wind axis is invariably produced (Mortimer 1974; Imberger 1985; Horn et al. 2001). Subsequently, when the wind relaxes, the tilted isopycnal surfaces begin to steepen and a long internal wave packet develops with counterpropagating structures within the closed basin (Monismith 1987; Horn et al. 2002).

### b. Benthic stimulation by long internal waves

Although the existence of packets of long internal waves in shallow, stratified bodies of water is well established and, as noted, some of the processes leading to their formation in shallow seas are in hand, comparatively little is known about the consequences of these wave packets for benthic processes: mixing, resuspension, vertical transport, dissipation, particle motions, etc. Our understanding of the mechanistic dynamics driving benthic processes stimulated by long internal wave motions in shallow environments, and some of their environmental consequences, is at best only at a formative stage. The focus of the present work is directed particularly toward revealing and clarifying issues related to the physics of benthic mixing stimulated by long internal waves, and to point the way for more detailed studies that might yield the important quantitative measures necessary to assessing the role of long internal waves for physical and biological processes.

When reference is made to long waves, the wavelength under consideration is to be long relative to the controlling dimension that limits the phase speed of infinitely long waves. This dimension in most cases corresponds to the mixed layer depth. However, use of the term “long wave” in shallow environments often implies that the wavelength is commensurate with, or even larger than, the total depth of the fluid column. In such cases the wave-induced particle motions penetrate down to benthic levels with only slight diminution in magnitude. Consequently, the presence of long waves in shallow seas implies that particle motions that are a significant fraction of the wave phase speed are induced in benthic regions of shallow seas. It is not uncommon for wave-induced particle velocities to exceed 10–30 cm s^{−1}, levels that impose stresses sufficient to stimulate significant benthic mixing.

An essential feature of wave-induced currents associated with a propagating wave packet is that they possess both spatial and temporal variations. Now, it is often assumed that the wave-induced vorticity generation at the bottom boundary and the resulting benthic stresses will vary in phase and scale with the local strength of the wave-induced current. However, any quasiperiodic variability of the induced current field under a propagating packet implies the existence of a variable pressure gradient impressed upon the wave-induced boundary layer. Furthermore, this pressure gradient is necessarily adverse in regions where the strength of the wave-induced current diminishes below a local extremum, and the presence of regions of adverse pressure gradient implies that these are also regions of enhanced boundary layer instability. Of even greater significance for benthic dynamics, local regions of separated flow can be expected in the footprint of a wave packet where the adverse pressure gradient exceeds a critical level, a level that depends on both environmental conditions and the character of the wave packet.

### c. Wave-induced global instability

The point of particular interest in this work is the possible formation of local regions of separated flow in the wave-induced boundary layer. The local boundary layer profile in such regions will contain pockets of reversed flow and, because of the presence of inflectional profiles, the existence of an essentially inviscid mechanism of instability. Now, it is established that the growth rates of inviscid instabilities (e.g., those associated with free shear flows) are substantially larger than those that arise through a viscous mechanism (e.g., boundary layer flows) (Hammond and Redekopp 1998). The difference in temporal growth rate is typically at least an order of magnitude. Moreover, the inviscid instability in shear layers can potentially become of absolute type while the viscous instability in attached boundary layers is always of convective type (Hammond and Redekopp 1998). When the slower stream in a shear layer moves counter to the faster stream with sufficient strength, the instability character changes from convective to absolute type (Huerre and Monkewitz 1985; Huerre 2000).

The implications of these different characteristics for the flow in different portions of the wave-induced boundary layer under a long-wave packet are quite profound. To the point, the possible appearance of local regions of boundary layer separation where the inflectional profile is absolutely unstable may give rise to the spontaneous onset of global instability. Since the induced flow associated with long internal waves varies slowly and since the growth rate of (inviscid) instabilities in a flow with inflectional profiles is quite rapid, unstable global dynamics can likely grow to saturation levels before the local flow state tends toward relaminarization when the pressure gradient reverses and becomes favorable. Global instability is characterized by a well-defined frequency leading to a dynamical state that is coherent in space and time (Chomaz 2004; Huerre 2000). In the context of benthic layers under long internal waves, this dynamics will be manifest as a rollup of the boundary layer vorticity into coherent patches that, by arguments based on image vortices below the boundary, will tend to travel away from the boundary. Also, the dynamics giving rise to the formation of concentrated vortex structures will lead to accompanying regions of locally elevated bottom stress and enhanced mixing and dissipation.

### d. Previous numerical investigations

Several recent numerical studies have examined the temporal evolution of the wave-induced boundary layer under internal solitary waves (Bogucki and Redekopp 1999; Wang and Redekopp 2001; Stastna 2001; Stastna and Lamb 2002b). The simulations by Bogucki and Redekopp (1999) and by Wang and Redekopp (2001) considered a weakly nonlinear solitary wave of elevation in a two-layer stratification with an oncoming linearly sheared flow. Further, the solitary wave of elevation was implicitly assumed to be one formed by a local topographic resonance between an upstream-propagating wave and the sheared current. This means that the wave is almost stationary (at least on time scales for instability of a free shear layer) with respect to the bottom and any consequent dynamics has ample time to act over the length of a solitary wave. The particular choice of wave state coincided with the observed flow state reported by Bogucki et al. (1997) where the appearance of such waves of elevation was shown to stimulate pronounced episodes of resuspension with significant levels of sedimentary material detected at 12 m above the bottom surface in a fluid depth of 60 m.

The global instability arising in the direct numerical simulations reported by Bogucki and Redekopp (1999) and Wang and Redekopp (2001) emerged in a case where the mean vorticity of the oncoming flow dominated the strength and vertical scale of the vorticity generated in the (separated) wave-induced boundary layer. The simulations revealed a sudden onset of global dynamics as the wave amplitude exceeded a threshold level, and had the appearance of coherent vortex structures rising into the water column and being advected with the flow. The simulations by Stastna and Lamb (2002b), on the other hand, considered a fully nonlinear solitary wave of elevation in a slightly diffuse thermocline situated close to the bottom. The wave-induced boundary layer competed somewhat with an assumed form of a weak, tidal flow boundary layer. The separated flow, formed in the adverse pressure gradient region near the front of the wave, again appeared to be globally unstable, with the formation of vortex structures that moved up and into the deformed thermocline under the wave crest, and appeared to alter the wave structure. The same flow configuration was investigated in greater detail in the unpublished work of Stastna (2001). Stability bounds were derived for a phase space defined by the strength of the oncoming current *U _{S}* and the wave Reynolds number Re

*. The interaction of the solitary wave with the bottom boundary layer exhibited an increased tendency for global instability with increasing Re*

_{W}*,*

_{W}*U*and/or wave amplitude. Implications for sediment resuspension were considered through investigation of the dependence of bottom shear stress and vertical velocities on Re

_{S}*and*

_{W}*U*. These quite distinct simulations suggest that the benthic excitation under waves of elevation takes on a different character depending on the scale and strength of vorticity in the ambient environment into which the wave propagates.

_{S}From a physical standpoint, all of the above studies restrict consideration to solitary waves of elevation and to waves moving in environments where the benthic boundary layer forms in an essentially unstratified fluid. The case of a wave of elevation in the absence of a mean current was considered only by Stastna (2001) and Stastna and Lamb (2002b) where no global instability was observed for a wave amplitude (i.e., maximum isopycnal displacement) of *α*_{0} = 0.19 and Re* _{W}* as high as 7.4 × 10

^{6}. Waves of depression (of either mode 1 or 2) were not considered in any of these studies, presumably because of the prevalent perception that such waves are responsible for interfacial, but not benthic, mixing owing to their proximity to the top surface (Moum et al. 2003).

From a computational perspective, two outstanding issues exist in regard to the aforementioned numerical work. First, Bogucki and Redekopp (1999) and Wang and Redekopp (2001) provide only a brief description of the problem formulation. In a similar vein, Stastna (2001) does not elaborate on how the choice of a frame of reference moving with the phase speed of the wave affects the implementation of the bottom no-slip boundary condition and the advective operator in the governing equations. Second, one may wonder why Stastna (2001) and Stastna and Lamb (2002b) did not observe a global instability in the absence of a mean current. A possible answer may be that the wave amplitudes were simply not sufficiently large to initiate any instability in the wave-induced boundary layer. However, such a possibility may be highly unlikely for *α*_{0} ≈ 0.2 and Re* _{W}* ≈

*O*(10

^{7}) (see section 7).

The issues brought forth above motivate further consideration of the role of the accuracy of the numerical discretization scheme. Although the unsteady boundary layer forming in the wave footprint is necessarily viscous, the appearance of a global instability in the separated boundary layer is essentially “inviscid” and characterized by a length scale that decreases with increasing Re* _{W}* (Hammond and Redekopp 1998). The viscous sublayer is not of primary importance here and its detailed resolution would only result in unnecessary high computational cost. Therefore, a physically underresolved simulation, which uses a

*highly accurate*numerical method to capture only the inviscid dynamics of the separated flow, is the optimal choice for investigating such dynamics over as broad a range of Re

*values as possible. Nevertheless, the accurate enforcement of the no-slip boundary condition is absolutely critical to the correct reproduction of the physics of interest. Ensuring a smooth transition from such a strongly enforced boundary condition to the interior of the flow and the sharp gradients of the boundary layer does pose its own set of numerical challenges (Hesthaven and Gottlieb 1996). In principle, in this case, given that interpolating polynomials are the foundation of most every numerical discretization scheme (Boyd 2001), the number of interpolating modes dictated by sufficient resolution of the inviscid physics may simply not be enough to preserve numerical stability (Deville et al. 2002) (see also section 3). High-order accuracy (e.g., spectral) methods (Deville et al. 2002) are those most imminently faced with such a challenge because of their weak or nonexistent artificial dissipation. In contrast, a low-order accuracy scheme, such as the second-order finite-difference technique used by Stastna (2001), will produce a smooth solution at very high Re*

_{W}*(10*

_{W}O^{7}), as was the case with Stastna (2001) and Stastna and Lamb (2002b). This deceptive numerical stability results from the well-known inherent artificial dissipation of low-order finite-difference schemes, which could nonetheless generate significant error in the long-time integration of the governing equations (Abarbanel and Ditkowski 1999). In addition,

It is in no way guaranteed that the number of grid points used with such a low-order scheme provides sufficient accuracy in the resolution of the desired physics.

The artificial dissipation inherent in such a scheme, which cannot be quantified explicitly and spans a much broader range of scales than its counterpart associated with a high-order spectral filter, is quite likely to “smear out” critically important physics at the smallest resolved scales (e.g., the onset of global instability). Such a degradation in the accuracy of the representation of sharp localized features has been demonstrated for shocklike features (Garnier et al. 1999) and is most likely to be operative in low-order discretizations of high Reynolds number boundary layers.

An increase in resolution of the low-order scheme, beyond the conventionally used factor of 2, would be needed to determine the independence of resolution in terms of accurate reproduction of the physics against the damping effect of the artificial dissipation (Boyd 2001). In contrast, in a similar situation, a spectral scheme captures accurately the desired physics with at least a factor of 4 less resolution (Boyd 2001). The only obstacle for a spectral method would be to preserve its numerical stability at high Re* _{W}* without degrading its high accuracy. For incompressible flow simulations in vertically finite domains this capability was only recently made possible through the development of spectral multidomain penalty methods (Diamessis et al. 2005), which will be used to determine whether a global instability is indeed generated in the footprint of both waves of depression and elevation.

### e. A detailed numerical study of solitary wave–induced global instability

Motivated by the discussion of section 1d and the issues it gives rise to, both physical and numerical, the present study aims to explore further via the use of direct numerical simulations the character of benthic excitation by internal solitary waves for different wave polarities and modal structure in the presence of a diffuse bottom stratification (which could act as a damper of any boundary layer instability). In particular, we explore first the instability formed under waves of depression. In this case, the wave-induced vorticity has a sign opposite to that generated under a wave of elevation. The wave-induced current at the bottom under a wave of depression moves downstream relative to the wave crest, and the region of separation (when the wave amplitude is sufficiently large) occurs aft of the crest (i.e., in the region of adverse pressure gradient). The sign of the wave-induced vorticity at the bottom surface suggests, using an argument based on image vortices positioned below the surface, that any instability-generated vortex structure will move upward into the water column and downstream from the wave crest. We also explore further the case of instability under waves of elevation. The wave-induced vorticity in this case is generated by a relative flow (relative to the wave crest) moving in the direction of wave propagation (i.e., upstream). The global instability of the separated flow now develops from the tail of the wave forward toward the front of the wave, and the rollup and release of vortex structures will interact baroclinically with the front of the wave. As a result, any instability-generated flow state will exhibit quite different dynamics depending on the wave polarity and the base stratification.

The problem definition, governing equations, and boundary conditions are given in section 2. A summary of the numerical methodology is found in section 3. Given that this is only the second time such a scheme is applied to the study of small-scale geophysical flows at high Reynolds number, a slightly more expanded discussion is available in appendix B for the more interested reader. An overview and justification of wave forcing is provided in section 4 The calculation of the vertical eigenfunctions and the horizontal amplitude function of the wave is outlined in detail in appendix A. Details of the specific simulations performed are given in section 5. Results for waves of depression and elevation are presented in sections 6 and 7, respectively. A summary and thoughts for future work are given in section 8.

## 2. Model formulation

### a. Problem definition

The basic flow considered in this investigation is a solitary internal wave propagating with phase velocity *C* along a waveguide containing a structured thermocline in a region of uniform finite depth *H* (Fig. 1). When desired, a weak background shear corresponding to the mean barotropic tide may also be introduced. Observations of thermoclinic structure in coastal regions and in lakes motivate the choice of a Brunt–Väisälä frequency profile given (in dimensional form) by

where

Therefore, the ambient stratification consists of a top mixed layer, a seasonal thermocline with peak Brunt–Väisälä frequency *N*_{1} centered at level *z*_{0} and having a thickness scaled with *δ _{T}*, and a diffuse bottom stratification having a Brunt–Väisälä frequency of (1 + Δ

*)*

_{N}*N*

_{2}. The parameter Δ

*(0 ≤ Δ*

_{N}*≤ 1) allows the introduction of a weaker stratification in the upper layer but has been assigned the value Δ*

_{N}*= 1 in all cases considered herein. Other parameter values selected for the present study are (*

_{N}*α*,

*β*,

*δ*,

_{T}*z*

_{0}) = (0.8, 0.2, 0.2, 0.75

*H*) for a wave of depression and (

*α*,

*β*,

*δ*,

_{T}*z*

_{0}) = (0.9, 0.1, 0.4, 0.35

*H*) for a wave of elevation.

The computational domain is rectangular with dimensions *L _{x}* ×

*H*. The vertical configuration of the computational domain corresponds to an idealization of the coastal oceanic waveguide with an impermeable no-slip bottom and an impermeable stress-free top surface. Periodicity is assumed in the

*x*direction. The use of longitudinal periodicity does not allow for the examination of a spatially evolving flow such as a propagating internal wave in a laboratory-fixed reference frame. To overcome this restriction and investigate a temporally evolving flow, a steady wave-fixed frame is used. Use of a moving reference frame necessitates modifications of the governing equations, which are discussed in detail in section 2b. The assumption of longitudinal periodicity is valid provided the solitary wave possesses a wavelength sufficiently short with respect to the domain length

*L*to prevent any spurious interactions between the leading and trailing edges of the wave. Note that, if three-dimensional physics are of interest, a periodic spanwise direction may easily be introduced.

_{x}### b. Governing equations

The goal of the present study, as already outlined, is to examine the dynamics of the temporally evolving, wave-induced boundary layer under a propagating internal wave. It is convenient for computational purposes to examine this wave-induced boundary layer in a wave-fixed frame. To this end we consider the computational base state to consist of a waveguide translating with a uniform speed *C* to the right (i.e., a flow state with exact correspondence to a physical waveguide with the same internal wave propagating with speed *C* to the left). That is, the computational base state consists of a prescribed internal wave that is stationary in a uniform flow moving with speed *C* to the right. Further, the waveguide boundaries in this computational base state are also translating with speed *C* to the right. This is necessary in order to ensure that a boundary layer develops only when and where a wave-induced mismatch in horizontal velocity appears at the lower boundary of the waveguide.

The governing equations are the Navier–Stokes equations under the Boussinesq approximation (Winters et al. 2004; Diamessis et al. 2005):

and

where

The nonlinear term in the momentum equation in (3) is written in the skew-symmetric form to minimize aliasing effects in the numerical solution (Deville et al. 2002). In addition, for the purpose of time advancement of the discrete equation in (3), the buoyancy term has been absorbed into **N**(**u**). The quantities *p*′ and *ρ*′ are the perturbations of the pressure and density from their respective (mean) reference values, which are in hydrostatic balance (Diamessis et al. 2005).

The forcing terms *F*^{u}_{w} and *F ^{ρ}_{w}* appearing in the equations are determined by employing a particular decomposition of the fields for the velocity and the state variables. The disturbance velocity field (

*ũ*,

*w̃*) forced by the wave-induced boundary layer is measured relative to the specified wave field (

*u*,

_{w}*w*) and the mean current [

_{w}*C*+

*U*(

*z*)], where

*U*(

*z*) is a possible background shear flow satisfying, in general, the condition

*U*(0) = 0. The disturbance fields (

*ρ̃*,

*p̃*) for the density and pressure are measured relative to the hydrostatically balanced base state plus the fields (

*ρ*,

_{w}*p*) associated with the prescribed wave. Hence, the dependent variables appearing in the governing equations are decomposed into the following parts:

_{w}We note that all wave fields (*u _{w}*,

*w*,

_{w}*ρ̃*,

_{w}*p̃*) are stationary because of the particular frame of reference selected for the computational problem. Also, the specific expressions for these wave fields that were used in generating the simulations described later are defined in appendix A. The pressure and density disturbance fields (, ) are measured relative to the hydrostatically balanced base state plus those disturbances (

_{w}*p̃*,

_{w}*ρ̃*) consistent with the specified permanent wave. Substituting (7) into the system (3)–(5) yields the modified form of the Navier–Stokes equations:

_{w}and

The advective operator 𝒜 is defined as

Note that satisfaction of the solenoidal character of the wave velocity field has been taken into account. In addition, no wave forcing terms appear in (8) and (9) on account of the requirement that the wave velocity–density field be a solution of the Navier–Stokes equations in the interior of the computational domain in the absence of any perturbations. Specifically, the forcing terms appearing in (3) and (4) are equal to

The formulation as given assumes that the flow is planar. If a fully three-dimensional disturbance field is allowed, a spanwise velocity component and the lateral momentum equation must be introduced in (3)–(5) and included in (8)–(10). We specifically restrict our consideration of plane waves whereby the wave-induced boundary layer will, at least in early stages of evolution, also be two-dimensional. Subsequent evolution of the boundary-generated vorticity may exhibit three-dimensional structure, but characterization of issues such as benthic mixing and dissipation are left for future studies.

### c. Boundary conditions

The boundary conditions used in the numerical model correspond to the idealized waveguide description of the computational domain given in section 2a. In the horizontal direction, periodic boundary conditions are employed:

The top boundary is a free-slip nondeformable surface:

Choice of the bottom boundary condition for the disturbance velocity *ũ* is crucial to obtain a correct simulation of the unsteady, wave-induced boundary layer. Since, as noted earlier, the boundaries of the waveguide are moving with speed *C* to the right, the application of the no-slip condition on the bottom boundary leads to the condition:

It is the above mismatch in *u* velocity at the bottom boundary that leads to the formation of the wave-induced boundary layer and subsequent global instability. Last, the density perturbation is subject to a Dirichlet boundary condition at both vertical boundaries:

The boundary conditions for the pressure are of purely numerical nature and their discussion is thus deferred to appendix B.

## 3. Numerical method: Overview

The numerical methodology used is a recently developed spectral multidomain penalty method for the simulation of high Reynolds number incompressible flows in vertically finite domains. A full discussion of the numerical scheme and its validation (through comparison of simulations of stratified turbulent wakes with nonzero net momentum with corresponding laboratory data) may be found in Diamessis et al. (2005). appendix B provides some details of the temporal and spatial discretization, spectral penalty methods, and spectral filtering. Only a brief overview is provided here.

The temporal discretization of the governing equations combines third-order stiffly stable and backward-differentiation schemes with a dynamic high-order boundary condition for the pressure (Karniadakis et al. 1991). Thus, maximum temporal accuracy is attained and splitting errors at the vertical boundaries are minimized, as *O*(Δ*t*^{2}) accuracy is achieved in both velocity and pressure. In the periodic horizontal direction, Fourier spectral discretization is used with *N̂ _{x}* Fourier modes. In the vertical direction, the computational domain is partitioned into

*M*subdomains of variable height

*H*and order of polynomial approximation

_{k}*N̂*(

_{k}*k*= 1, . . .

*M*) (Fig. 2). The total number of vertical grid points is

*N̂*=

_{z}*M*(

*N̂*+ 1) + 1. Within each subdomain, Legendre spectral discretization (Boyd 2001) is used and, for the specific problem under consideration,

*N̂*is fixed and equal to a fixed

_{k}*N̂*in all subdomains. Subdomains communicate with their neighbors via a simple patching condition (Boyd 2001). Among others, the primary advantage of the multidomain discretization (Hesthaven 1997) is flexibility in local resolution, which allows positioning of a sufficient number of grid points in the boundary layer (and to a lesser degree in the seasonal thermocline) combined with minimal resolution of the less active ambient fluid outside the boundary region (Fig. 3).

For the resolutions considered in this paper, dictated by available computational resources and the need for rapid run turnaround, the simulations become underresolved at higher Reynolds numbers. Although the vertical resolution of the bottom boundary is carefully chosen to capture all phenomena critical to the dynamics of the inviscid region of the unsteady boundary layer, the viscous sublayer is not captured. As a consequence, if an orthogonal polynomial-based spectral discretization is used with no auxilliary stability-ensuring methodology (e.g., filters, penalty schemes), viscosity is not felt by the resolved scales and numerical instabilities develop because of lack of sufficient interpolating polynomial modes and subsequent aliasing effects (see also the discussion in section 1d). Spurious energy with increasingly higher and higher frequency content is then generated, producing a catastrophic effect on the long-term integration of the governing equations (Gottlieb and Hesthaven 2001). To ensure stability of the numerical solution while preserving its spectral accuracy, penalty techniques (Hesthaven and Gottlieb 1996; Hesthaven 1997) are used in the vertical direction along with strong adaptive interfacial averaging (Don et al. 2003). As a final safeguard against numerical instability, explicit spectral filtering is used in all spatial directions (Levin et al. 1997).

In summary, in regard to the flow under consideration, the above numerical scheme provides two distinct advantages. First, the temporal discretization, when combined with a spectral spatial discretization, provides maximum temporal accuracy at the boundaries (Guermond and Shen 2003) and the solution is not contaminated by the formation of splitting-induced spurious numerical boundary layers (Boyd 2001). Second, the penalty scheme allows the simulation of the inviscid physics of the wave-induced boundary layer [provided sufficient grid points are inserted therein (Boyd 2005)] without requiring resolution of the thin viscous sublayer, which would require a prohibitively high number of grid points. Thus, the internal (internal with regard to the viscous sublayer) inviscid high Reynolds number physics of the flow are captured with spectral accuracy without having to worry about underresolution-driven numerical instabilities. No enhanced numerical viscosities, which will bias the internal flow dynamics when using a spectral scheme (Boyer et al. 2004), are required. In addition, the excessive artificial dissipation of a low-order finite-difference scheme, which can damp features at the smallest resolved scales essential to the global instability (see the discussion of section 1d), is no longer an issue.

## 4. Wave forcing: Overview and justification

The vertical eigenfunction of the internal solitary wave in a general waveguide of depth *H* is given by the Taylor–Goldstein equation. For the purpose of illustrating some of the qualitative features of wave-induced global instability, only long internal waves described by weakly nonlinear theory are considered. The details of specific forcing calculations may be found in appendix A.

The choice of a weakly nonlinear [i.e., Korteweg–de Vries (KdV)] representation for the wave field (*u _{w}*,

*w*,

_{w}*p̃*,

_{w}*ρ̃*) entering the forcing functions in (12) and (13) is not essential to our approach and might be criticized from several perspectives. It is known that more accurate characterizations are necessary in order to capture the state of larger-amplitude observations in the ocean (Stanton and Ostrovsky 1998; Stastna and Lamb 2002b). Even fully nonlinear waves obtained as solutions of the Euler equations could be employed (Turkington et al. 1991). However, the primary intent of this work is to expose the possible presence of global instability in the footprint of long internal waves of depression in shallow seas and to illustrate the qualitative dynamics associated with the mature state of this wave-induced instability. The objective is not to generate parameterizations of resuspension, dissipation, and mixing arising from this instability for subsequent use in larger-scale ocean models. The generation of such parameterizations requires much more extensive and ambitious simulations than are envisioned here. It is our perspective that before such studies are initiated, a qualitative insight into the likely appearance of globally unstable eruptions of the wave-induced boundary layer, and of typical dynamics associated with such eruptions, is needed. The qualitative understanding that we seek requires principally a wave forcing containing a local extremum in the isopycnal displacement in order to drive a wave-induced boundary layer and its concomitant pressure field that is impressed on the developing boundary layer. Now, it is certainly desirable that the scaling characteristics of the imposed wave field mimic with reasonable fidelity those associated with long internal waves in typical environments. Weakly nonlinear theory provides analytical expressions for the wave field that captures at least the salient features outlined here, and does so in a relatively simple and straightforward manner.

_{w}Higher-order theories relevant to larger amplitudes nominally yield solitary waves with flatter crests (or troughs in the case of waves of depression) and possess steeper fronts and tails than those defined by KdV theory (Lamb and Wan 1998). As such, the pressure gradients will generally be stronger than those predicted by KdV theory as the amplitude increases. As a consequence, the adverse pressure gradients would increase more rapidly with amplitude and the onset of boundary layer separation and global instability could then be expected at even lower amplitudes than predicted here. Of course, Lamb and Wan also show that profiles of strongly nonlinear waves vary somewhat depending on the structure of the Brunt–Väisälä profile, and in some specific cases can even mitigate against the formation of steeper fronts and tails in larger amplitude waves.

It is true that wave amplitudes giving rise to global instability based on weakly nonlinear theory and the Brunt–Väisälä profile considered in this work are not necessarily in consonance with rigorous upper bounds for solitary waves, especially those without recirculating cores (Stastna and Lamb 2002a; Lamb 2002). Nevertheless, the presence or absence of a region of local overturning is only relevant to the present objectives through its possible influence on the shape of isopycnal surfaces in the lower portion of the water column, that is, regarding how the near-bottom wave structure influences the formation and dynamics of the wave-induced boundary layer. We believe that weakly nonlinear theory provides a minimal representation of any adverse pressure gradient region aft of a wave of depression and does so via a rather simple scaling relation in terms of the wave amplitude. Further, since the simulations reported here are based on the time-dependent equations with full viscous and diffusive effects, any simulation containing a trapped core is, in fact, a fully adjusted solution, albeit one subject to a rather specific forcing function. The simulations, therefore, are expected to capture the true dynamics of the lower part of the water column where the impressed pressure gradient is described approximately by that prescribed by KdV theory. We also note that the vorticity in any trapped core appearing in the simulations presented here is always found to be much weaker than any vorticity generated in the separating wave boundary layer and that, for a mode-1 wave of depression, it is located significantly above the boundary layer (e.g., see Figs. 4a and 4b). Hence, any velocities induced on the separating flow by a trapped wave core are expected to have at most only a marginal effect on the breakdown dynamics of the separation bubble. Last, the presence of a vortex core does shift the location of the wave front but has a negligible effect on the steepness of the front (Derzho and Grimshaw 1997) and, therefore, the induced adverse pressure gradient.

Given the above arguments, our choice of wave amplitudes is not unrealistic in terms of the corresponding values observed in wave packets in the coastal ocean: Bourgault and Kelley (2003) observe waves of depression with amplitudes *α*_{0} ≈ 0.25 in an estuary and Scotti and Pineda (2004) have identified waves of elevation that lift the near-bottom-situated seasonal thermocline in Massachussetts Bay to about *half* the depth of the oceanic waveguide. Moreover, the recent northern South China Sea study of Duda et al. (2004) reports waves of depression that, prior to shoaling, achieve a near-permanent form with amplitudes of 160 m or more at 350-m depth. Recent unpublished measurements off the coast of North West Australia (G. Ivey 2005, personal communication) have revealed trains of waves of depression where the leading wave displaces the thermocline by about 80 m in a 125-m-deep ocean!

The simple choice of wave forcing used here can also be motivated from another perspective. Virtually all long internal waves appearing on the shelf arrive as evolving wave packets instead of isolated, equilibrium structures (solitary waves). Furthermore, either wave packets or isolated waves generally undergo continuous deformation arising from the space-varying nature of the waveguide (e.g., shoaling waves due to depth changes). Such evolving waves, at least insofar as their induced velocity and pressure near the bottom of the waveguide are concerned, may not differ so significantly from a KdV description. In any case, some approximate representation of these fields is required, and the present choice is a reasonable first-order description that should elucidate the dynamics of globally unstable separations induced in the footprint of long internal waves and allow the ensuing dynamics to be linked to typical threshold wave conditions.

## 5. Simulation description

Simulations of the wave-induced unsteady boundary layer are carried out using the spectral multidomain penalty method model of Diamessis et al. (2005) summarized in section 3 and appendix B. A summary of the parameter values, grid resolutions, and spectral filter orders considered in each simulation is provided in Table 1. The governing nondimensional parameters of the problem are the wave Reynolds number Re* _{W}* =

*C*

_{0}

*H*/

*ν*, the nondimensional wave amplitude

*α*

_{0}[see Eq. (A8)], and the Prandtl number Pr =

*ν*/

*κ*, where

*ν*and

*κ*are the kinematic viscosity and diffusivity of the fluid, respectively. Re

*spans the range [10*

_{W}^{4}, 10

^{5}]. Note that the Re

*used to characterize runs with either mode-2 waves, which exhibit a reduced phase speed*

_{W}*C*

_{0}relative to their mode-1 counterpart, is that defined in terms of the value of

*C*

_{0}for the corresponding mode-1 wave. All “production” runs use Pr = 1. Higher Pr values, representative of fresh or saline water, would have a considerable impact on the development of the small scales of the density field. However, these small scales, where mixing is most important, are expected to have a negligible impact on the development of the global instability, regardless of Pr value, and are not captured by the resolutions considered in this study. Given that the purpose of this study is not the study of small-scale mixing, the choice of Pr = 1 is deemed sufficient.

The computational domain of height *H* has a length *L _{x}* = 10

*H*, which allows the exponential decay of all soliton waveforms considered here. A uniform grid is used in the horizontal direction whereas in the vertical a spectral multidomain discretization (shown in Fig. 2) is used with

*M*= 10 nonuniform subdomains of uniform order of approximations

*N̂*=

_{k}*N̂*. The positioning of the subdomains (Fig. 3b) is made 1) to ensure sufficient resolution of the boundary layer (Boyd 2005), a region critical to the dynamics of the problem, and 2) to provide smooth but not excessive resolution of the seasonal thermocline. An estimate of the thickness of the boundary layer may be obtained using simple scaling arguments (Kundu 2002), where the characteristic horizontal length scale and velocity are the solitary wave wavenumber

*k*

^{−1}

_{s}(see appendix A) and the maximum wave-induced velocity [given in Eq. (A1)], respectively. In the meantime, minimum resolution is used for the nearly nonactive top mixed layer. Higher Re runs utilize a subdomain configuration that is more closely clustered in the boundary layer in comparison with the Re

*= 2 × 10*

_{W}^{4}cases. In addition, a higher order of approximation is employed in each subdomain (

*p*refinement). When moving up to Re

*= 5 × 10*

_{W}^{4}, doubling the horizontal resolution is inevitable. At even higher Re

*,*

_{W}*p*refinement is employed in the vertical while a reduced filter order is used in the horizontal to suppress any underresolution-driven numerical noise at the highest Fourier modes. In general, the orders (

*p*,

_{F}*p*) of the spectral filters used in both vertical and horizontal directions are chosen to minimize aliasing effects but to also leave a sufficient number of modes in the numerical solution (i.e., at least 50% of the modes in the horizontal direction and 8–9 modes in each subdomain) not directly affected by the filter. Thus, the spectral accuracy of the numerical scheme is preserved (Gottlieb and Hesthaven 2001). A more detailed discussion of the role of explicit spectral filtering with respect to the effective Re

_{L}*and the specific damping it imposes on the numerical solution may be found in appendix B.*

_{W}The grids used in each simulation are chosen to provide sufficient resolution of the inviscid dynamics of the global instability. Full resolution of the scales where dissipation and mixing are important is beyond the scope of this paper and, as already indicated in section 3, any numerical instability that may develop as a result of neglecting these scales in the spectral discretization is damped through the penalty scheme and spectral filtering. Figure 4 demonstrates the adequacy of resolution used at Re* _{W}* = 2 × 10

^{4}. Keeping in mind that the primary objective of this study is to provide a qualitative description of global instability dynamics, only minor differences are observed in the vortical structures formed in the separation bubble (soon after the onset of global instability) when resolution is doubled. These minor differences, most evident in the exploded view of the separation bubble (Figs. 4c and 4d), such as more “rolled up” vortices in the higher-resolution case or the presence of a third, more prominent, vortex at the tailend of the bubble, may be due to differences in temporal resolution between regular- and high-resolution runs and the availability of a greater range of scales in the latter. The configuration with

*N̂*= 384 and

_{x}*M*= 10 subdomain of

*N̂*= 16 configuration employed at Re

*= 2 × 10*

_{W}^{4}is used as a template for the choice of resolutions/subdomain positioning at all other Re

*values. To verify the adequacy of resolutions chosen, in all runs, horizontal Fourier spectra and vertical (in each subdomain) Legendre spectra at all Re*

_{W}*were examined and revealed no signs of aliasing, that is, energy accumulation at the smallest resolved scales (Boyd 2001). An example of such spectra, corresponding to the flow fields in Figs. 4a and 4b, is shown in Fig. 5. No qualitative difference between the two spectral curves is observed. The spectrum of the regular-resolution run is a little less smooth possibly because the vertical averaging is performed over roughly half the number of grid points than in the high-resolution case. For the purpose of a qualitative analysis, a more robust indicator of resolution independence is that both curves in Fig. 5 exhibit a common distinct peak at a wavenumber that corresponds to the characteristic wavelength of the global instability and the separation between like-signed vortices in Figs. 4a and 4b.*

_{W}The computational time step Δ*t* was chosen as such that the CFL stability criterion be obeyed in both spatial directions for a third-order stiffly stable scheme (Karniadakis et al. 1991). The following requirements are imposed:

In all runs, the initial value of Δ*t* is set by (18). Immediately before the appearance of the global instability, and throughout its development, the energy in the separation region, and thus the vertical velocity, increase drastically. Successive reductions in Δ*t* must then be made. Efficient handling of the decrease in time step at the required time levels is provided by an adaptive time-stepping technique (Diamessis et al. 2005) based on a third-order variable Δ*t* Adams–Bashforth and backward differentiation schemes.

As indicated in section 2b, a third spanwise direction (chosen periodic for facility of implementation) may be easily introduced in the problem formulation and, thus, the *υ* velocity may be solved for. Although three-dimensional effects are not the focus of this study, test simulations at supercritical Re = 10^{4} and 2 × 10^{4} were performed with the addition of a spanwise dimension of width *L _{y}* = 0.5

*H*and

*N̂*/

_{y}*L*=

_{y}*N̂*/

_{x}*L*. Nonnegligible

_{x}*υ*velocities, and thus three-dimensional behavior, were observed at approximately three nondimensional time units (

*tC*

_{0}/

*H*= 3) after the onset of global instability. Thus, a two-dimensional simulation does provide a reliable description of the onset of global instability and the formation of coherent concentrations of vorticity that are predominantly parallel to the wave crest.

## 6. Global instability of a propagating internal solitary wave of depression

### a. Fundamentals and stability bounds

As indicated in section 1e it is the aim of this section to fill in some gaps existing in the types of wave forcing considered in previous numerical investigations (Bogucki and Redekopp 1999; Wang and Redekopp 2001; Stastna 2001; Stastna and Lamb 2002b) and to provide a more complete description of the qualitative dynamics of global instability. A global instability (Hammond and Redekopp 1998; Huerre 2000) arises spontaneously in a spatially developing flow when a pocket of local absolute instability (i.e., a pocket wherein a band of instability waves exists with group velocities oriented both upstream and downstream) of sufficient streamwise length is allowed to form. Synchronous dynamics of strong spatiotemporal coherence are observed to arise spontaneously over the entire span of the unstable region at marginally supercritical conditions. The disturbance field associated with the instability imposes a discrete set of frequencies and wavelengths on the base flow and behaves as an oscillator with intrinsic dynamics, in contrast to a convectively unstable flow (e.g., a coflowing shear layer), which may be viewed as a noise amplifier with extrinsic dynamics (Huerre 2000). Thus, the unstable region is quite insensitive to any disturbances advected into its interior.

In the unsteady boundary layer induced by a propagating solitary internal wave, a pocket of absolute instability develops when the reversed flow generated by the adverse pressure gradient (APG) associated with the longitudinal gradient of the horizontal wave velocity becomes sufficiently deep. The APG, which scales as *A*^{3/2}_{0} under a KdV solitary wave, gives rise to a relatively symmetric separation bubble and the formation of a reversed flow as the amplitude rises. Temporally averaged profiles may be obtained at the center of the separation bubble, with averaging performed from the beginning of the separation (the point where the maximum *u* velocity in the backflow region is about 10% of the phase speed *C*_{0}, say) until the onset of global instability (the point where oscillatory behavior is observed in any temporal signal sampled in the separation region). Any such profiles (such as those shown in Fig. 6 for a solitary wave of depression) may then be analyzed using linear instability theory [specifically the methodology outlined in Hammond and Redekopp (1998)] to determine conditions for the onset of local absolute instability in the center of the separation bubble. Of course, the base-state boundary layer in the present formulation is not stationary but grows progressively in time as the simulation is set up as an initial value problem. Figure 6 clearly shows that an increase in wave amplitude *α*_{0} corresponds to a more rapidly appearing separated backflow of increasing magnitude.

A stability boundary (Hopf bifurcation) diagram may now be constructed in the *α*_{0} – Re* _{W}* plane (Fig. 7). The construction of such a stability curve calls for an appropriate definition of the onset of global instability. In this study, a different definition is used than that proposed by Stastna (2001). It is assumed that the solitary-wave-induced boundary layer is unstable if distinct shorter wavelength (with respect to the horizontal extent of the bubble) disturbances are observed in the separation bubble before time

*tC*

_{0}/

*H*= 8. In all simulations, whenever such spatial oscillations are observed, the separation bubble is guaranteed to fragment and shed vortices into the water column within no more than two nondimensional time units. Although an increase of the arbitrarily selected cutoff time of

*tC*

_{0}/

*H*= 8 would likely lead to a slight reduction of critical amplitude, it must be emphasized that global instability is a threshold phenomenon: increasing subcritical wave amplitudes do not produce an ever-weaker global instability. To the contrary, no global instability is observed and the separation bubble remains intact, even after a long time has elapsed, when the wave amplitude is subcritical. In Fig. 7, a marked reduction in critical wave amplitude as a function of increasing Re

*is observed. This drop in critical amplitude may be attributed to both a thinner boundary layer and a stronger backflow resulting from the higher Re*

_{W}*. As indicated in Fig. 7, higher Re*

_{W}*, more relevant to the ocean (e.g., in the range from 10*

_{W}^{5}to 10

^{7}), are anticipated (Wang and Redekopp 2001) to asymptote to a constant critical value

*α*

_{0cr}not vastly different from solitary wave amplitudes observed in the coastal ocean (Stanton and Ostrovsky 1998; Bourgault and Kelley 2003; Duda et al. 2004; G. Ivey 2005, personal communication). Specifically, an exponential best fit to all the available data points on the stability curve yields the relation

Consequently, values of Re* _{W}* = 10

^{6}or 10

^{7}would exhibit critical values of

*α*

_{0cr}equal to 0.29 and 0.22, respectively. It is reasonable, therefore, to believe that the description of the global instability process presented here is directly pertinent to the ocean.

### b. Visualizations of global instability

A visualization of the formation and subsequent destabilization of the separation bubble induced by a solitary wave of depression at Re* _{W}* = 2 × 10

^{4}at critical and at a strongly supercritical wave amplitude (

*α*

_{0}= 0.46 and 0.55, respectively) is given in Fig. 8. Prior to the onset of global instability, a regular separation bubble has formed in the footprint of the solitary wave with a vorticity structure consisting of a “tongue” of negative perturbation vorticity overlying a thin region of positive value (Figs. 8a and 8d). The different signs of perturbation vorticity arise from the different levels of the separated boundary layer. A backflow with negative vorticity forms adjacent to the boundary and a layer of positive vorticity forms associated with the inflectional velocity profile in the shear layer along the dividing streamline.

A completely spontaneous dynamics occurs as the global instability imposes its own reduced characteristic length scale on the separation bubble. This characteristic length scale manifests itself through the creation of distinct concentrations of coherent vorticity in the interior of the bubble and the emission of weaker finescale waves away from the bubble in the downstream (positive *x*) direction (Figs. 8b and 8e). The formation of such distinct oscillatory features fore and aft of the trailing edge of the bubble is the physical manifestation of an internal excitation source located in the bubble, a unique signature of local absolute instability as stated in the beginning of this section. Eventually, the bubble breaks up into a sequence of coherent eddies of smaller wavelength (Figs. 8c and 8f). This reduction in horizontal scale of the coherent dynamics leads to a sharpening of the horizontal pressure gradient and, thus, enhanced bottom shear stress. At the same time, a discrete frequency pulsation is imposed upon the coherent dynamics and a process of vortex shedding is initiated (shown in greater detail in Fig. 10, described below). The downstream emission of weaker instability waves continues until the ejection of the first vortex (Figs. 8c and 8f). Thereafter, the whole structure of the separated region is significantly altered as time progresses and the global instability approaches a saturated nonlinear state. The weak downstream-propagating instability waves tend to “break” farther downstream with negligible consequences on the overall flow dynamics.

A sequence of contour plots, similar to that in Fig. 8, for a higher value of Re* _{W}* = 10

^{5}is shown in Fig. 9. The wave amplitudes displayed are chosen as such to exhibit the same degree of criticality for each Re

*. The global instability is established in a manner similar to that observed in the lower Re*

_{W}*simulations (Fig. 8). However, distinct differences are also noticeable. As expected, the boundary layer is much thinner. The downstream-emitted weak instability waves are now characterized by a considerably shorter wavelength as illustrated by Figs. 9b and 9f. The emitted instability waves propagate into the convectively unstable boundary layer downstream of the separation region (Figs. 9c and 9g) where they are destabilized in a more vigorous fashion than their lower Re*

_{W}*counterparts. The upstream-propagating waves quickly establish an imprint on the separation bubble and cause it to fragment into coherent vortices of a smaller length scale than the Re*

_{W}*= 2 × 10*

_{W}^{4}case, presumably because of the shorter wavelength of the instability and the thinner boundary layer.

Values of Re* _{W}* less than 2 × 10

^{4}were also examined (Table 1). In this case, viscous diffusion is responsible for larger separation bubbles and the subsequent formation of shed vortices of larger size but weaker magnitude. Although not shown here, no global instability was observed in exploratory runs at the lowest Reynolds number of Re

*= 5 × 10*

_{W}^{3}and at significantly high amplitudes (e.g.,

*α*

_{0}= 0.55). Instead, the separation bubble was fragmented into smaller scales that were unable to leap into the water column, evidently because of viscous effects. The fragmented separation bubble remained in place as the vortices scoured the bottom, but only to a height approximately equal to the top of the bubble. In conclusion, one may reasonably assume that global instability does not occur below a certain critical Reynolds number.

An investigation of the global instability at Re* _{W}* = 10

^{5}and

*α*

_{0}= 0.55 would have been desirable to assess the role of the Reynolds number at a given wave amplitude. However, it was found that the resolution typically used for lower amplitude Re

*= 10*

_{W}^{5}runs was insufficient to capture the thinner boundary layer and shorter instability wavelengths and the required increase in resolution implied prohibitively high computational cost. Nevertheless, in terms of the potential of global instability for resuspension and its depedendence on Reynolds number, significant insight may be gained by examining simulations at Re

*= 2 × 10*

_{W}^{4}and 10

^{5}, critical and strongly supercritical wave amplitudes at the

*same level of criticality*, and times well after the onset of global instability.

^{1}In Figs. 10a–d a string of oppositely signed vortices has been ejected into the water column, each vortex mixing the local density field and likely transporting material in its core. The vortices rise high into the water column until they reach a trapping depth and/or radiate their energy away in the form of small internal waves. The maximum height of ascent of the vortices appears to be independent of

*α*

_{0}for Re

*= 2 × 10*

_{W}^{4}, although

*α*

_{0}= 0.55 exhibits a shorter streamwise separation between successive vortices because of the higher intrinsic wavenumber of the global instability (see Figs. 11a and 11d). A fundamental difference of the Re

*= 10*

_{W}^{5}runs, with respect to their lower Re

*counterparts, is the release of smaller-diameter vortices at higher frequencies into the water column. Eventually, as shown in Fig. 10d and best illustrated by animations (not shown here), the separation bubble at high Re*

_{W}*is annihilated, although it never ceases to emit downstream instability waves (i.e., it remains globally unstable). Under the persisting effect of the APG of the large-scale internal wave the bubble appears to reconstruct itself, which could lead to a new cycle of vortex shedding. Such behavior could imply that, at higher Re*

_{W}*, the vortex shedding associated with the global instability is not continuous but rather consists of temporally intermittent bursts of near-vertically ejected vortices, the duration of which is commensurate with the lifespan of the separation bubble.*

_{W}### c. Quantitative results

Vertically averaged streamwise Fourier wavenumber spectra (where the averaging is performed over the temporally evolving height of the separation bubble) are shown in Fig. 11. Along with the visualizations of Figs. 8 and 9, these spectra allow the identification of different growth regimes during the lifetime of the separation bubble. Initially, the energy in the boundary layer grows slowly. However, as soon as the first signs of shorter wavelengths are apparent in the separation bubble (Figs. 8b,e and 9b,c,f,g) a rapid increase in the energy level of the Fourier spectra is observed. Although the initial appearance of such wavelengths in physical space is not accompanied by a distinct signature in Fourier space (most probably due to the filtering effect of the averaging), as soon as the separation bubble begins to fragment into distinct coherent vortices, a visible peak is established at intermediate wavenumbers in the Re* _{W}* = 2 × 10

^{4}simulations. This spectral peak is directly associated with the separation between vortices prior to the initiation of the shedding process (see also Fig. 5). That such distinct peaks are not as evident in the Re

*= 10*

_{W}^{5}spectra is possibly due to the more vigorous breaking of the bubble and the downstream-emitted instability waves.

The spectra in Fig. 11 indicate that the downstream-propagating global instability waves are not an artifact of aliasing, which would manifest itself in “spectral blockage” (Boyd 2001), that is, an accumulation of energy at the highest resolved wavenumbers. To the contrary, the spectra drop off smoothly by three orders of magnitude to the value of the wavenumber where the spectral filter has a direct effect (see Fig. 17). Although the instability waves are probably filtered out in the vertically averaged spectra, spectra sampled at distinct vertical locations where these waves first appear (not shown here) do exhibit associated maxima, which, for all Re* _{W}* considered, are sufficiently removed from the range of modes directly influenced by the spectral filter. In addition, visualizations for the higher resolution (768 × 250) runs at the lower Re

*= 2 × 10*

_{W}^{4}value (not shown here) do not display any observable difference in the structure of these instability waves. As is also evidenced by Figs. 9b and 9f, upon their appearance in the Re

*= 10*

_{W}^{5}runs, the instability waves exhibit a separation between peaks in physical space equal to

*x*/

*H*= 0.1, that is, roughly eight grid points, which is sufficiently larger than wavelengths equal to the grid spacing, a classical symptom of aliasing (Boyd 2001). Thus, we are confident that the instability waves are a physical effect.

The effects of the spontaneous coherent dynamics on the bottom perturbation shear stress for the supercritical Re* _{W}* = 2 × 10

^{4}and 10

^{5}cases are illustrated in Fig. 12. At both Re

*values, upon onset of global instability, a sudden growth of an order of magnitude in ∂*

_{W}*ũ*/∂

*z*and

*w̃*is observed, associated with the formation of shorter-length-scale horizontal vortices in the separation bubble's interior. As soon as the shedding of vortices into the water column is initiated, ∂

*ũ*/∂

*z*is relaxed to a smaller value. The shear stress in Fig. 12 is normalized by Re

^{1/2}

_{W}(associated with the boundary layer thickness). The similar levels of normalized ∂

*ũ*/∂

*z*at both Re

*demonstrate the scaling of the shear stress with boundary layer length scales. Last, the vertical velocities remain significantly high and comparable to the phase velocity*

_{W}*C*

_{0}of the solitary wave upon and after the onset of global instability and appear to be independent of Re

*.*

_{W}The temporal signature of discrete frequency pulsation associated with the global instability dynamics is obvious in Fig. 13, which shows time series of *ũ* and *w̃* velocities at two different streamwise locations in the separation bubble of a weakly supercritical solitary wave at Re* _{W}* = 2 × 10

^{4}. A dominant frequency, characteristic of the global instability, is evident in the time signals. In addition, the in-phase oscillations at the two different locations demonstrate that the coherent dynamics are synchronous over the entire length of the bubble.

### d. Mode-2 waves

Mode-2 internal solitary waves, though less energetic than mode-1 waves, are frequently observed in lakes and the coastal ocean (Imberger 1998; Bogucki et al. 2005). These waves exhibit similar global instability characteristics as their mode-1 counterparts (see Fig. 14) at corresponding Re* _{W}*. However, because of the sharper gradient of the eigenfunction at the bottom (Fig. 16) and their inherently shorter wavelength and sharper horizontal pressure gradient, mode-2 waves are characterized by a significantly lower critical amplitude (at Re

*= 2 × 10*

_{W}^{4},

*α*

_{0cr}= 0.28) and more a rapid onset of global instability.

## 7. Wave of elevation

This section provides a qualitative description of the evolution of global instability for the case of a mode-1 solitary wave of elevation. As discussed in section 1e, the generation and evolution of wave-induced vorticity is a function of wave polarity. Thus, different behavior is expected than in the case of a wave of depression. Figure 15 displays the global instability due to a propagating solitary wave of elevation. The stronger gradient of the vertical eigenfunction at the bottom of the waveguide, arising from the lowered position of the seasonal thermocline, is responsible for a reduced value of critical wave amplitude (for Re* _{W}* = 2 × 10

^{4},

*α*

_{0cr}= 0.42). In this case, the wave-induced APG is located near the front of the wave and the center of the separation bubble is located immediately below the center of the wave. Since the wave-induced velocity is upstream relative to the wave crest, positive-signed vortices are formed in the bubble and shed upstream of the wave (not downstream as for the case of a wave of depression). Eventually, the vortices reach the strong seasonal thermocline and produce significant density perturbations. Baroclinically generated vorticity then forms and imposes significant distortions to the horizontal wave structure and the associated APG that controls the formation of the separation bubble. Hence, the dynamics of the separated region are altered by the interaction of the global instability with the solitary wave.

Assuming that a wave of elevation exhibits a stability boundary similar to that in Fig. 7, a value of *α*_{0cr} = 0.42 at Re* _{W}* = 2 × 10

^{4}indicates that at Re

*≈*

_{W}*O*(10

^{7}) and wave amplitudes comparable to

*α*

_{0}= 0.19 a global instability is very likely to occur. The role of a vortex core could have a more direct influence on the dynamics of global instability. However, at

*α*

_{0}= 0.19 the vorticity in the core is quite weak relative to that induced in the boundary layer. In addition, as in the case of a wave of depression (Derzho and Grimshaw 1997), the presence of a vortex core may induce a shift in the position of the adverse pressure gradient but will not necessarily diminish its magnitude. Thus, it is legitimate to believe that the absence of global instability in the work of Stastna and Lamb (2002b) may be a consequence of the smoothing effect of the high artificial dissipation inherent in a low-order finite-difference scheme (see the discussion of section 1d).

## 8. Summary and concluding remarks

The boundary layer induced by a solitary internal wave in shallow water has been examined through direct numerical simulation (DNS) of the incompressible Navier–Stokes equations under the Boussinesq approximation. Both waves of depression (modes 1 and 2) and elevation (mode 1 only), derived through weakly nonlinear theory, are considered. The mean density field corresponds to that observed in the coastal ocean and lakes, consisting of a top mixed layer and a seasonal thermocline overlaying a diffuse stable stratification. A high-order splitting scheme of the Navier–Stokes equations allows maximum temporal accuracy in the numerical solution and minimizes any splitting errors at the boundaries. Fourier spectral discretization is used in the periodic longitudinal direction. In the vertical direction, a highly flexible Legendre spectral multidomain discretization scheme provides spectral accuracy and adequate resolution of the inviscid dynamics of the boundary layer. At higher wave Reynolds number Re* _{W}*, for the resolutions considered, the use of a spectral multidomain penalty scheme and explicit spectral filtering suppresses any numerical instabilities driven by underresolution. The use of a numerically stable and spectrally accurate scheme enables the reliable representation of any finer features in the smallest resolved scales of the flow that can be responsible for the development of a physical instability. Such features would be smoothed out if a low-order finite-difference scheme were used. A highly accurate description of the qualitative dynamics of the wave-induced boundary layer and its consequent global instability is therefore possible for values of Re

*as high as 10*

_{W}^{5}.

Results of filtered DNS show that, for all Re* _{W}* > 5 × 10

^{3}and a wave amplitude exceeding a critical value, a global instability occurs in regions near the bottom boundary where the wave-induced current diminishes below a local extremum and an adverse pressure gradient develops, leading to a sufficiently strong reversed flow. If this reversed flow is of sufficient strength, a pocket of local absolute instability appears, giving birth to a spontaneous (requiring no external perturbation) dynamics. The instability imposes its own reduced characteristic length scale on the separation bubble and the dynamics is coherent in the sense that it extends along the entire extent of the bubble. As the instability develops, the bubble is quickly fragmented into a sequence of coherent eddies of reduced wavelength. This decrease in length scale leads to a sharpening of the horizontal pressure gradient near the bottom and, consequently, enhanced bottom shear stress. Once the global instability has been established, the separation bubble pulsates synchronously (i.e., the entire unstable region oscillates in phase) at a discrete frequency characteristic of the instability and a periodic shedding of vortices into the water column is initiated. The vertical velocities of the vortical structures are comparable to the phase velocity

*C*

_{0}. At Re

*= 10*

_{W}^{5}or higher, this periodic vortex shedding is not continuous in time, but appears to manifest itself in temporally intermittent events, each consisting of a number of successive vortex ejections. Although details of this periodic train of vortices do depend on the sign of the wave-induced vorticity, itself a function of the wave polarity, what is always observed is an ascent of the shed vortices upward into the water column, sometimes approaching a level as high as 30%–35% of the depth of the computational domain In a 60-m coastal ocean, for example, this would correspond to a height approximately 20 m above bottom, perhaps carrying resuspended particulate matter to the same level.

It is stressed again that global instability is a threshold phenomenon: increasingly subcritical wave amplitudes do not generate global instability. However, this study has established that the critical wave amplitude *α*_{0cr} decreases with increasing Re* _{W}* and, therefore, it is reasonable to expect the mechanisms of global instability described in this paper to appear frequently at values of Re

*and*

_{W}*α*

_{0}

_{cr}typically encountered in the coastal ocean and lakes. Given the highly energetic content of the associated coherent vortical structures and the maximum height above bottom that they attain, it is most likely that this specific instability, and not the bursts and sweeps typical of the nominal turbulent boundary layer (Head and Bandopadhyay 1981; Moin and Mahesh 1998), is the driving mechanism behind sediment resuspension observed in shallow waters (Bogucki et al. 1997, 2005). Such resuspension phenomena of bottom depositions of sand, nutrients, and chemicals have profound implications for the optical and acoustical properties of the water column, water quality and hygiene, and the balance of local ecosystems. It is also natural to expect that three-dimensional secondary instabilities and turbulence subsequent to the onset of global instability are what dominate benthic mixing and dissipation. A next step would thus be to utilize well-resolved three-dimensional DNS at as high a Re

*as possible and Pr ∈ [1, 7] (to attain values of this parameter that are at least of limnological relevance) to elucidate and parameterize the details of the wave-induced boundary turbulence, mixing, and resuspension. The spectral multidomain penalty scheme is ideally suited for such a pursuit; the increased resolution required for DNS automatically transforms the penalty scheme into a regular spectral multidomain discretization where both boundary and patching conditions are enforced strongly (Hesthaven 1997; Diamessis et al. 2005).*

_{W}To our best knowledge, there appears to be a glaring paucity of available field measurements on internal solitary-wave-induced boundary layers in shallow waters. Certain studies have documented the bulk features of particulate resuspension in the passage of internal solitary waves (Bogucki et al. 1997; Dickey et al. 1998; Bogucki et al. 2005). Lake measurements (Lemckert et al. 2004) have reported “active turbulent benthic boundary layers” generated and maintained by basin-scale waves. Coastal data obtained from a shoaling wave of elevation (Klymak and Moum 2003) suggest, but do not resolve, a “highly sheared bottom layer and stress-driven turbulence.” However, no existing field datasets resolve details of the global instability and the associated dissipation and mixing. In addition, the perception shared by observational oceanographers appears to be that it is simply the action of the wave-induced current sweeping over the bottom that generates such an energetic boundary layer. It is our hope that this study will provide the incentive for new detailed investigations of such energetic benthic flows and will provide additional insight for the interpretation of the basic physics.

The boundary layer induced by a propagating internal solitary wave in a waveguide of uniform depth is not the only benthic boundary flow energized by internal waves where global instability is of relevance. It is anticipated that a global instability in the shoaling of an internal solitary wave is intimately connected to the modal conversions experienced by the propagating wave, and that the benthic mixing and turbulence produced in such a case are comparable to or even stronger than the corresponding interfacial phenomena near the seasonal thermocline (Moum et al. 2003). Internal waves reflecting off a critical slope have been a popular subject of recent observational (Eriksen 1998; Nash et al. 2004), laboratory (De Silva et al. 1997; Dauxois et al. 2004), and numerical studies (Slinn and Riley 1998b; Tabaei et al. 2005). Although overturning and shear instabilities are commonly proposed as the driving mechanism behind such a flow, it is more than likely that the most energetic dynamics are driven by global instability in the wave-induced boundary layer. Utilizing the advantages of our spectral multidomain penalty solver (and any extensions of it to implement open boundary conditions in the longitudinal direction), our future plans include further simulations to investigate the role of global instability in both of the above flows.

## Acknowledgments

We are grateful to Jan Hesthaven and Andrzej Domaradzki for useful insight on the adaptation of the numerical model to the specific problem. Discussions with Leon Boegman and Kevin Lamb at the First ISCOR Conference on Ocean Mixing were particularly illuminating on the physics of the specific problem. We thank Greg Ivey for preliminary data on internal solitary waves off the coast of NW Australia. PJD is thankful to Otto Rehhagel for providing inspiration during the course of writing this paper, through the perfect blend of enthusiasm and focus. Geoff Spedding kindly provided access to additional computing resources at the USC Fluid Mechanics Laboratory. This work was supported by ONR Grant N00014-03-1-0476.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Wave Forcing

##### Wave flow field

Following the problem formulation described in section 2, the two-dimensional wave-induced velocity field (*u _{w}*,

*w*) and the pressure and density perturbations (with respect to the ambient density and hydrostatic pressure fields)

_{w}*p̃*and

_{w}*ρ̃*, respectively, of the solitary internal wave are given by their leading-order forms for a weakly nonlinear longwave solution of the Euler equations (Maslowe and Redekopp 1980):

_{w}and

In these expressions, *ρ*_{0} is a reference density, *g* is the acceleration of gravity, *ϕ*(*z*) is the vertical eigenfunction of the solitary wave, and *A*(*x*) is an amplitude function indicating the wave's horizontal structure. It is obvious from (A1) that the wave velocity field satisfies conservation of mass. Details of the computation of *ϕ*(*z*) and *A*(*x*, *t*) are given in the next subsection.

##### Details of calculation

For an internal solitary wave in a general waveguide of depth *H* with background velocity exhibiting at most only linear variation, the eigenvalue problem in the long wave limit is given by the Taylor–Goldstein equation:

Note that hereinafter we will work only in the long internal wave limit and thus the wave velocity will be denoted by the corresponding long-wave velocity *C*_{0}. In addition, because this specific study focuses on propagating internal solitary waves in the absence of a background current, *U*(*z*) = 0 in all equations that follow, which greatly simplifies their representation. Equation (A5) may be discretized on a one-dimensional single domain of height *H* with an *N̂*th-order Gauss–Lobatto–Legendre grid (Boyd 2001). Following some straightforward manipulations, the discretized equation may be written in the form of a generalized eigenvalue problem (Boyd 2001). The first *N̂* eigenvalues and eigenvectors are then obtained by inputting **A** and **B** into a generic *QZ* algorithm routine. The eigenfunctions *ϕ*(*z*) are normalized by their maximum dimensional value. Only the first two eigenmodes are retained. Different values of *N̂* have been tested and *N̂* = 64 seems to be adequate for the needs of this work. Last, the eigenfunction *ϕ*(*z*) and its two derivatives are interpolated on the final multidomain grid used in the main simulation. Figure A1 shows the eigenfunctions for the three basic cases considered in this work: modes 1 and 2 of a wave of depression and mode 1 of a wave of elevation (see also Fig. 1).

Since only weakly nonlinear waves are considered in this study, the horizontal structure of the solitary wave [given by the function *A*(*x*, *t*); see (A1)–(A5)] may be obtained by solving the KdV equation:

In a frame of reference moving with the wave, (A7) has the solution

where *A*_{0} = *α*_{0}*C*_{0}*H* and the wavenumber *k _{s}* is

It is straightforward to show that the nondimensional wave-induced streamline displacement is given by *ζ*(*x*, *z*) = *α*_{0}*A*(*x*)*ϕ*(*z*). Thus, following (A8), *α*_{0} is a direct measure of the maximum streamline displacement, which occurs at the location of the maximum value of *ϕ*(*z*). The coefficients *α _{k}* and

*β*, in the absence of a background current, are computed from the following formulas (Maslowe and Redekopp 1980):

_{k}where the normalization coefficient *I _{k}* is given by

Equations (A9)–(A11) clearly demonstrate the nonlinear nature of the solitary wave in that wavenumber *k _{s}* is a function of its own amplitude as well as environmental conditions [i.e., the stratification profile

*N*(

*z*)].

### APPENDIX B

#### Numerical Method: Some Specifics

##### Navier–Stokes splitting scheme

For the temporal discretization of (3)–(5), a high-accuracy pressure projection scheme is used (Karniadakis et al. 1991). According to this scheme, if one integrates (3)–(5) in time from level *t _{n}* to

*t*

_{n}_{+1}, the following semidiscrete equations, decomposed into three fractional steps for

**u**, are obtained (Diamessis et al. 2005):

In the above, the advective operators on the left-hand side of (8) have been absorbed into the nonlinear operator **N**(**ũ**). The values of the coefficients *α _{q}*,

*β*, and

_{q}*γ*

_{0}for the third-order backward differentiation stiffly stable (BDF3-SS3) scheme of (B1)–(B3) may be found in Karniadakis et al. (1991). The quantity

*ϕ*

^{n}^{+1}given by

is an intermediate scalar field that ensures that the final velocity **ũ**^{n+1} is incompressible. In (B2), it is assumed **∇** · = 0 and the Poisson equation is solved for the pressure:

The boundary conditions for **ũ **(15)–(16) are enforced in (B3). Equation (B5) utilizes the high-order accuracy dynamic boundary condition (Karniadakis et al. 1991):

##### Spectral multidomain penalty methods

In underresolved simulations, the absence of any artificial dissipation in the spectral discretization scheme leads to instabilities due to aliasing effects (Gottlieb and Hesthaven 2001) closely linked to underresolved numerical/physical boundary layers (Deville et al. 2002) (see section 3 for the role of the viscous sublayer in this regard). Specifically, Gibbs oscillations develop, which are most pronounced at the physical boundaries and subdomain interfaces and may generate catastrophically interacting artificial internal waves in a stratified flow. The crux of penalty methods (Hesthaven and Gottlieb 1996; Hesthaven 1997) lies in the fact that numerical instabilities arise in underresolved simulations because boundary/patching conditions are explicitly enforced and there is no provision that the equation solved is satisfied arbitrarily close to the boundary–subdomain interface. By collocating a linear combination of the equation and boundary/patching conditions (the latter multiplied by a penalty term) at the corresponding spatial locations, the penalty method produces a smooth numerical solution with near-negligible error at the boundaries and interfaces. The spectral accuracy of the numerical scheme is negligibly impacted and one may compute stably the high Re internal dynamics of the flow without having to resolve the thin numerical/viscous physical boundary layers.

In terms of the splitting scheme outlined in section a of appendix B, the penalty method is applied at two different levels in the incompressible Navier–Stokes equations. The explicit nonlinear term advancement is treated as a hyperbolic equation, whereas the implicit viscous term treatment is as a parabolic equation [in this section all subsequent equations are written as a function of the *u*-velocity perturbation (dropping the tilde) without loss of generality]:

The temporal derivatives in (B7)–(B8) are only approximations to those appearing in (B1)–(B3). In both cases, the patching conditions are such that *C*_{0} and *C*_{1} continuity^{B1} is enforced only weakly since the interface of two adjacent subdomains corresponds to the same location in physical space but to two separate grid points. Thus, the penalty method is inherently discontinuous and it is this weak continuity on the numerical grid that allows for a more stable and smooth numerical solution at the subdomain interfaces.

Consider first (B7). The boundary conditions, imposed strictly in (B3), are not incorporated in the penalty treatment of the hyperbolic problem, which focuses solely on the subdomain interfaces. The penalty formulation of (B7), in *physical space*, for each subdomain of index *k* and uniform order *N _{k}* =

*N*is

where

Here *δ _{ij}* is the Kronecker delta function with subscript

*i*corresponding to the collocation point

*z*

^{k}

_{i}. The terms in brackets represent appropriate patching conditions (Diamessis et al. 2005). The values of the coefficients

*α*,

*γ*,

*τ*

^{k}

_{1}, and

*τ*

^{k}

_{2}are determined by treating each individual subdomain as a single domain whose interfaces support patching conditions that act as open boundary conditions through which information is exchanged with adjacent subdomains. Essentially, the patching conditions are treated as localized open boundary conditions where each subdomain experiences “inflow” or “outflow” depending on the value of the vertical interfacial velocities

*W*

^{k}

_{0}and

*W*

^{k}

_{N}at the previous time steps (Hesthaven and Gottlieb 1996; Hesthaven 1997; Diamessis et al. 2005).

The penalty formulation of (B8) is significantly different from that of (B7) because of its parabolic nature. After some manipulation (Diamessis et al. 2005), the discrete form of (B8) may be written in Fourier space for an individual longitudinal wavenumber *k _{x}* as

where *u* is the value of the velocity at time step *n* + 1. The small parameter *ε* is defined as

The penalty formulation of (B11) is now

where *Q*^{−}_{k}(*z*^{k}_{i}) and *Q*^{+}_{k}(*z*^{k}_{i}) are defined in (B10). Again, the terms in the brackets represent appropriate patching/boundary conditions determined, in this case, by the elliptic nature of (B11). At the physical boundaries, *g*^{1}_{1}(*t*) and *g*^{M}_{2} are determined by the the boundary conditions (16)–(17). The coefficients *α*, *β*, *γ*, and *δ* are set to 1 at each subdomain interface, and at the physical boundaries their values are set to satisfy the physical boundary conditions. The coefficients *τ*^{k}_{1} and *τ*^{k}_{2} for the discretized equations are determined by the type of the boundary/patching condition (Dirichlet, Neumann, Robin) under consideration (Hesthaven and Gottlieb 1996; Diamessis et al. 2005).

##### Spectral filtering, adaptive interfacial averaging

Penalty methods are most effective in the vicinity of subdomain interfaces. To ensure stability of an underresolved simulation in the subdomain interiors, spectral filtering is used. In each subdomain, any function *f* (*z*) may be approximated as (Boyd 2001)

where *f̃ _{j}* is the amplitude of the

*j*th Legendre mode and

*P*(

_{j}*z*) the

*j*th-order Legendre approximating polynomial (Boyd 2001). The filtered solution

*f*may now be expressed in terms of the modes of the numerical solution as

^{F}where *k _{j}* is the

*j*th discrete Legendre mode. In this study, an exponential filter function

*σ*(

*k*) is used (Gottlieb and Hesthaven 2001):

_{j}where *p* is the filter order, *k _{c}* is the filter lag, and

*α*= −ln

*ε*with

_{M}*ε*being the machine precision. In all simulations considered in this paper,

_{M}*k*= 0. Application of the filter function (B16) is equivalent to introducing a

_{c}*p*th-order hyperviscous operator in the governing equations (Boyd 1998; Gottlieb and Hesthaven 2001). The hyperviscous operator is multiplied by an artificial viscosity

*ε*≈ (1/

_{N}*N*

^{p}^{−1}) (Don et al. 2003). The advantage of the explicit application of a spectral filter is that it avoids the additional stiffness and subsequent time step limitations of a hyperviscous term.

An analogous filtering procedure is applied in Fourier space. In the incompressible spectral multidomain simulations presented in this paper, spectral filtering is applied at three levels when advancing the solution from time level (*n*) to level (*n* + 1), as described in detail in Diamessis et al. (2005). The order of the Legendre and Fourier filters is the same in all three levels.

Clearly, spectral filtering is inevitable in the simulations discussed in this paper. The scales where molecular viscosity is important, and can play an additional role as a damper, are much smaller than the numerical cutoff. Thus, the spectral filter acts as a needed dissipation mechanism that ensures numerical stability (Gottlieb and Hesthaven 2001). As shown by the spectral filter functions in Fig. B1, spectral accuracy is preserved because the direct effect of the filter never extends below the top 50% of the numerical solution modes in Fourier space and the first 8–9 Legendre modes in each subdomain. That is, the accuracy of the scheme remains at least *O*(Δ*z*)^{8} in the vertical direction, while exponential convergence of the solution is always guaranteed in the horizontal. Use of an explicit spectral filter does imply an unavoidable reduction of the effective Reynolds number. However, such a consequence is not only limited to higher-order methods (Levin et al. 1997; Slinn and Riley 1998a). Low-order finite-difference schemes are not damping free, even if no filter is used. In fact, for resolutions typically employed in high Reynolds number simulations, their inherent artificial dissipation affects a broader range of resolved scales (Domaradzki et al. 2002) than the spectral filter functions in Fig. B1, thereby reducing even more the effective Reynolds number. In addition, no explicit control is possible over the amount of artificial dissipation introduced into the numerical solution. On the contrary, spectral filtering does possess this capability of explicit control as the filter order may be adjusted to damp adequately, but not excessively, any numerical noise due to underresolution.

Caution is admonished in using the artificial viscosity coefficient, *ε _{N}*, suggested by Don et al. (2003) to compute an effective Reynolds number (Winters and D'Asaro 1989). Such an attempt would imply attaching a physical significance to the spectral filter beyond its numerical role as a damper, especially when it is not totally clear what units should be assigned to the hyperviscosity in the vertical direction where the form of the corresponding hyperviscous operator is quite complex (Boyd 1998; Don et al. 2003). Computing the actual effective eddy viscosity as a function of numerical mode is quite an involved process (Domaradzki et al. 2002) that extends beyond the scope of this paper.

The last stabilization technique employed in the flow solver is strong adaptive averaging at the subdomain interfaces (Don et al. 2003; Diamessis et al. 2005). Filtering is least effective at the subdomain interfaces and, on some rare occasions, the penalty method may be unable to damp singularities in these locations, which appear in the form of spikes of increasing steepness. If the jump in the inherently discontinuous numerical solution across a subdomain solution exceeds a certain threshold, a smoothing procedure is applied where the two interfacial values are replaced by the average of the grid points immediately above and below the interface. In this study, the averaging operation was never required for more than an *O*(0.1%) fraction of the subdomain interfaces.

## Footnotes

* Current affiliation: School of Civil and Environmental Engineering, Cornell University, Ithaca, New York

*Corresponding author address:* Peter J. Diamessis, School of Civil and Environmental Engineering, Cornell University, Ithaca, NY 14853. Email: pjd38@cornell.edu

^{B1}

Continuity of the solution and its normal derivative. Here, *C*_{0} should not be confused with the long-internal-wave phase speed.