A robust numerical method for the generation and propagation of periodic finite-amplitude internal waves in natural waters using high-accuracy simulations

The design and implementation of boundary conditions for the robust generation and simulation of periodic finite-amplitude internal waves is examined in a quasi two-layer continuous stratification using a spectral-element-method-based incompressible flow solver. The commonly used Eulerian approach develops spurious, and potentially catastrophic small-scale numerical features near the wave-generating boundary in a non-linear stratification when the parameter $A/\left(\mathit{\delta }c\right)$ is sufficiently larger than unity; A and δ are measures of the maximum wave-induced vertical velocity and pycnocline thickness, respectively, and c is the linear wave propagation speed. To this end, an Euler–Lagrange approach is developed and implemented to generate robust high-amplitude periodic deep-water internal waves. Central to this approach is to take into account the wave-induced (isopycnal) displacement of the pycnocline in both the vertical and (effectively) upstream directions. With amplitudes not restricted by the limits of linear theory, the Euler–Lagrange-generated waves maintain their structural integrity as they propagate away from the source. The advantages of the high-accuracy numerical method, whose minimal numerical dissipation cannot damp the above near-source spurious numerical features of the purely Eulerian case, can still be preserved and leveraged further along the wave propagation path through the robust reproduction of the non-linear adjustments of the waveform. The near- and far-source robustness of the optimized Euler–Lagrange approach is demonstrated for finite-amplitude waves in a sharp quasi two-layer continuous stratification representative of seasonally stratified lakes. The findings of this study provide an enabling framework for two-dimensional simulations of internal swash zones driven by well-developed non-linear internal waves and, ultimately, the accompanying turbulence-resolving three-dimensional simulations.

1 Introduction

Internal swash zones (ISZs) (Emery and Gunnerson, 1973; Woodson, 2018) are regions which develop along sloping oceanic boundaries through the action of periodically incident internal waves (IWs) in a manner analogous to a surface swash zone on the beach, albeit at slower timescales (O(10 min) or longer) and over longer wavelengths (O(1 km) or longer) (Cowen et al., 2003; Elfrink and Baldock, 2002). In ISZs, energy can effectively flux down the scale to turbulence through either shear or convective instabilities in the IW interior, similarly to spilling or plunging breaker waves on the ocean surface (Cowen et al., 2003; Ting and Kirby, 1996), or the turbulent boundary layer established through the interaction of the IW-induced current with the seafloor (Zulberti et al., 2022). In the latter context, particularly strong turbulence can be generated in the form of a near-bottom turbulent wake due to boundary layer separation associated with the along-bed wave-induced adverse pressure gradient induced by either internal bores or internal solitary waves (Hosegood et al., 2004; Boegman and Stastna, 2019). The above turbulence-generation mechanisms presumably conspire to drive a significant boundary–interior exchange (McPhee-Shaw and Kunze, 2002; McPhee-Shaw, 2006), i.e., the exchange of water between the boundary layer and the stratified interior, which effectively drives mixing in the relatively less active stratified waterbody interior (McPhee-Shaw and Kunze, 2002; Boegman and Ivey, 2009). The periodic shoaling and breaking of the above IWs in shallow environments, such as continental shelves or slopes, has a direct impact on the internal thermal equilibrium and biogeochemistry of the water column (Woodson, 2018). The periodically on-slope incident IWs are important in the transfer of mass, whether they transport nutrients and plankton toward the surface in the inner shelf (Omand et al., 2015) or eject bottom boundary layer sediments as high as 40 m into the water column during a strong vertical updraft event associated with the passage of a non-linear internal wave of depression over the slope (Cheriton et al., 2016). A similar class of long IW-driven phenomena of comparable biogeochemical importance also occurs on the slopes of lakes (Thorpe, 1998; Wuest and Lorke, 2003) and has served as the primary motivator of the research presented here.

The leading-order component of the periodic wave field forcing of an ISZ consists of a lower vertical-mode IW whose wavelength is O(50–100) longer than the water column depth – namely, in the form of oceanic internal tidal waves or the basin-scale internal seiche of a lake (Emery and Gunnerson, 1973; Nash et al., 2004; Stevens et al., 2005; Martini et al., 2013; Lemckert and Imberger, 1998). In the latter case, the internal seiche is further associated with a lower horizontal mode that is itself associated with the longer dimension of the lake. Such long waves are commonly expected to be represented with sufficient fidelity through the use of linear IW theory (Stastna, 2022) at, nonetheless, values of finite wave amplitude. Frequently, higher-frequency/shorter-wavelength highly non-linear features, such as turbulent/undular bores or internal solitary waves (Stastna, 2022), may be embedded within the longer incident IWs (Hosegood et al., 2004; Lucas and Pinkel, 2022; Thorpe et al., 1996).

The primary objective of this paper is the development of a robust numerical method for the generation and subsequent development of the longer component of the deep-water wave forcing at finite amplitude. The generated wave should have an amplitude that is not constrained by the limits of linear theory. Practically, this corresponds to wave-induced maximum isopycnal displacements that are at least 5 % of the total water column depth. The wave should also remain sufficiently robust near the source and, with an equal degree of robustness, non-linearly adjust its waveform as it propagates along the waveguide. In this regard, central to this paper's scope is that the background stratification extends beyond an uniform density gradient (Taylor, 1993) and is actually subject to variation in the vertical, as characterized by the presence of a distinct pycnocline which is commonly a close approximation of the in situ background profiles in the stratified ocean or lakes. Finally, an additional essential ingredient of this study, is that a high-accuracy discretization is used – specifically, a nodal spectral element method (Diamantopoulos et al., 2022). The particular discretization technique enables the optimal resolution of the generated waves, their non-linear adjustments away from the source and ultimately (though not explicitly considered here) the associated instabilities/turbulence upon encountering the waves with the slope.

In the laboratory, one approach to generating periodic long internal IWs is by tilting and releasing the actual laboratory tank (Boegman et al., 2005): the resulting horizontal standing wave is a lab-scale surrogate of the basin-scale internal seiche generated in a long stratified lake in response to a strong wind event (Boegman, 2009). An equivalent type of horizontal standing wave may be generated in a numerical simulation within a long rectangular computational domain using an initial condition consisting of a tilted pycnocline (Grace et al., 2019). One issue with the tilting-based wave-generation approach may be that it immediately produces finite velocities across the whole domain/tank when one would prefer waves propagating into an initially quiescent slope region. In many cases, the standing wave will break down into a propagating wave train (Grace et al., 2019).

An alternative, more flexible and effectively more controllable wave-generation approach involves introducing a form of deep-water (far from the slope) oscillatory wave excitation. Such an approach would ideally allow for a sufficiently long propagation distance in uniform depth waters, prior to the wave encountering the slope, which permits the generated IW to undergo any required non-linear adjustments. To this end, in the deep-water section of a laboratory tank, a horizontally oscillating paddle (Wallace and Wilkinson, 1988; Nakayama and Imberger, 2010; Ghassemi et al., 2022), a vertically oscillating semi-cylinder (Moore et al., 2016) or an array of plates vertically stacked on an eccentric camshaft (Mercier et al., 2010, 2013) have been used. It is worth noting that all the above experimental studies generated relatively short waves as represented by values of aspect ratio $\mathit{\lambda }/H$ and non-dimensional amplitude ${\mathit{\eta }}_{\text{max}}/H$; λ, η_max and H are the IW horizontal wavelength, IW-induced maximum isopycnal displacement and water depth, respectively. Reported directly, or inferred, values of $\mathit{\lambda }/H$ and ${\mathit{\eta }}_{\text{max}}/H$ lie in the range [2,12.5] and [0.0075,0.2], respectively, in the above laboratory studies wherever directly identifiable or inferable. Such a maximum pycnocline displacement range corresponds to a wave Froude number, $\mathit{Fr}=U_{\text{max}}/c$, with a value of [0.02,0.35]; U_max and c are the maximum wave-induced horizontal current and wave propagation speed.

The high-order-accuracy turbulence-resolving fully non-linear and non-hydrostatic three-dimensional simulations of Winters (2015) are one of the few computational studies so far which have considered the generation and incidence of a periodic long wave and on a relatively steep slope. The wave aspect ratio and wave Froude number can be inferred as $\mathit{\lambda }/H=\mathrm{48}$ and Fr=0.1385. Note that the work of Winters considered only a uniform background stratification. Moreover, Winters' generated waves were allowed to have a distance lower than one prescribed wavelength from the source to propagate until the slope most likely precluding any deep-water non-linear adjustments of the waveform.

To the authors' best knowledge, the only other computational study which has examined the periodic generation, the propagation away from the source over at least one wavelength and the incidence of long internal waves on a slope is the two-dimensional investigation by Dauhajre et al. (2021). The wave aspect ratios considered in this study are high and can be inferred to be residing in the range [50,400] while noting the very small ratio of computational domain depth to length. The wave-based Froude number values considered are in the range [0.05,0.4]. Furthermore, the subset of simulations that use a two-layer stratification (and not a linear one) have a thick pycnocline and focus on the aspect ratio of $\mathit{\lambda }/H=\mathrm{200}$ and a wave-based Froude number between 0.1 and 0.2. Note also that the curvature of the density profile at the base of the pycnocline is reduced by introducing a weakly stratified layer below. The numerical dissipation inherently built into the parameterizations of the regional ocean modeling code (Regional Ocean Modeling System, ROMS) used in this study could also effectively damp any spurious numerical features near the wave-generating deep-water boundary.

Note that the studies of Masunaga et al. (2015, 2016) and Walter et al. (2012) also considered a non-uniform stratification but positioned the wave-generating source only a fraction of the target wavelength from the slope. The generated waves, therefore, were not afforded an adequate propagation distance to undergo any non-linear adjustments before encountering the slope. Additionally, per this paper's focus on periodic IW simulation with high-accuracy methods and high resolution, the nesting-based robust mode-1 long internal tide generation within a regional-scale nonhydrostatic model (Rogers et al., 2019) is not pertinent to the scope of this study as it relies on low-pass filtering and sponge layers.

The laboratory and computational studies discussed above consider generated waves that may be deemed either short or long. Even when high-accuracy/resolution numerical methods are derived and efficiently implemented on a state-of-the-art high-performance computing platform, a computational study aiming to sufficiently resolve instability/turbulence formation due to sufficiently high-amplitude waves over a limited number cycles of an ISZ is practically limited to a wave aspect ratio in the range [40,50]. This is the aspect ratio regime accessed by the work of Winters (Winters, 2015), which is, however, limited to a linear stratification. The choice of a linear stratification effectively shielded this study from the challenges that emerge when forcing internal waves in a pycnocline-dominated stratification profile. As will be demonstrated later in this paper, the generation of high-amplitude periodic internal waves in more general, non-linear stratifications for waves operating in this intermediate aspect ratio range is confronted with non-trivial error if commonly used deep-water forcing approaches, such as those employed by Winters (2015), Dauhajre et al. (2021), are actually employed. The minimal numerical dissipation of a high-order-accuracy numerical method can allow this error to grow substantially. The stability of the simulation can thus be effectively undermined, and one can no longer leverage the high accuracy of the method for representing non-linear wave adjustments in deeper water and the finer-scale features once the slope is reached.

Thus, from a computational point of view, a relatively simple technique for generating larger-amplitude IWs for general stratifications in deep water is highly desirable. This is often achieved by choosing a form of boundary conditions at the boundary away from the slope region. For most field-relevant stratifications, a pycnocline dominates the stratification, and the vertical motion of the pycnocline is the clearest manifestation of internal waves. Historically, descriptions of internal waves typically built on a linearized theory and the literature have examples of two different choices for the vertical coordinate: one which uses the physical coordinate z and one which uses the upstream height of each isopycnal (and, more concretely, the upstream height of the dominant pycnocline) (Gear and Grimshaw, 1983; Yih, 1977). Since the former uses the physical coordinates, it is usually labeled as the Eulerian theory of linear internal waves. The latter, in contrast, is labeled the Euler–Lagrange theory because the horizontal coordinate is the physical coordinate x, while the vertical coordinate is the upstream coordinate, often written as $y=z-\mathit{\eta }(x,z,t)$, where η is the isopycnal displacement. Both the Eulerian (Lamb and Yan, 1996) and Euler–Lagrange theories (Gear and Grimshaw, 1983) have been used as a basis for multi-scale asymptotic expansions that extend the wave description to small but finite-amplitude waves (i.e., weak non-linearity) and waves of finite wavelength (i.e., weak dispersion). These lead to model equations in the Korteweg–de Vries family. The use of the upstream isopycnal height has found general use in the description of stratified flow in both the classical (Yih, 1977) and modern (Stastna, 2022) contexts. In the simulation context, the desire to generate finite-amplitude waves in a situation with a strong pycnocline implies that forcing methodologies based on the Eulerian linear wave theory may not yield robust results. The Euler–Lagrange theory offers an alternative, if algebraically more complex, development pathway.

In this paper, by following an Eulerian and an Euler–Lagrange approach (Turkington et al., 1991; Gear and Grimshaw, 1983), different types of time-dependent periodic wave-generating boundary conditions are derived with a particular emphasis on the subtleties associated with a continuous two-layer background stratification. The efficacy of each approach in generating a robust deep-water periodic finite-amplitude IW train and enabling any non-linear adjustments of the wave-train is thereafter assessed.

2 Problem setup and model formulation

2.1 Problem geometry

The canonical flow examined in this paper is the propagation of a two-dimensional finite amplitude periodically forced internal wave in a quasi two-layer continuous stratification. The computational domain is a two-dimensional rectangle of dimensions L×H and is stratified in vertical direction z with a vertically varying buoyancy frequency N(z), where

$$\begin{array}{}\text{(1)}& N^{\mathrm{2}}\left(z\right)\equiv -{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{\mathit{\rho }}}{\mathrm{d}z}}.\end{array}$$

Restricting one's focus to the Boussinesq approximation, the total density is decomposed as the addition of a reference density, ρ₀; a stratification, $\stackrel{\mathrm{‾}}{\mathit{\rho }}$; and a perturbation, ρ^′ (Kundu et al., 2015):

$$\begin{array}{}\text{(2)}& \mathit{\rho }(\mathit{x},t)={\mathit{\rho }}_{\mathrm{0}}+\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z\right)+{\mathit{\rho }}^{\prime }(\mathit{x},t),\text{with}{\mathit{\rho }}^{\prime }\ll \stackrel{\mathrm{‾}}{\mathit{\rho }}\ll {\mathit{\rho }}_{\mathrm{0}}.\end{array}$$

The quasi two-layer continuous stratification $\stackrel{\mathrm{‾}}{\mathit{\rho }}$ (Fig. 1a) is defined by

$$\begin{array}{}\text{(3)}& \stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z\right)=-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0}}{N}_{\mathrm{0}}^{\mathrm{2}}\mathit{\delta }}{g}}\tanh \left({\displaystyle \frac{z-z_{\mathrm{p}}}{\mathit{\delta }}}\right),\end{array}$$

where ${\mathit{\rho }}_{\mathrm{0}}{N}_{\mathrm{0}}^{\mathrm{2}}\mathit{\delta }/g$ is a measure of the density difference across the pycnocline, N₀ is a reference buoyancy frequency equal to the peak value of N(z) in the water column, δ is a measure of the pycnocline thickness and z_p is the position of the pycnocline's center.

.

Figure 1 (a) A two-layer continuous stratification, $\stackrel{\mathrm{‾}}{\mathit{\rho }}$, defined by its density jump, $\mathrm{2}{\mathit{\rho }}_{\mathrm{0}}\frac{N_{\mathrm{0}}^{\mathrm{2}}\mathit{\delta }}{g}$; its position, z_p; and its thickness, δ. (b) Corresponding vertical structure eigenfunction, W(z), which is the solution to the eigenvalue problem in Eq. (16).

The finite-amplitude internal wave of wavenumber k and angular frequency ω is generated through a forcing implemented within the left boundary conditions. Details on the exact derivation of the deep-water boundary conditions is covered in Sect. 3.

2.2 Governing equations

The governing equations for the problem are the incompressible Navier–Stokes equations (INSEs) under the Boussinesq approximation, written as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial \mathbf{u}}{\partial t}}=-\mathbf{u}\cdot \mathrm{\nabla }\mathbf{u}-{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{\rho }}^{\prime }\mathit{k}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{p}^{\prime }+\mathit{\nu }{\mathrm{\nabla }}^{\mathrm{2}}\mathbf{u},\text{(5)}& {\displaystyle \frac{\partial {\mathit{\rho }}^{\prime }}{\partial t}}=-\mathbf{u}\cdot \mathrm{\nabla }\mathit{\rho }+\mathit{\kappa }{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\rho }}^{\prime },\text{(6)}& {\displaystyle }\mathrm{\nabla }\cdot \mathbf{u}=\mathrm{0}.\end{array}$$

The simulations reported here will be limited to two dimensions on the (x,z) plane. Therefore, u will be limited to its two components, (u,w). Furthermore, ρ^′ is the density perturbation as defined in Eq. (2) and p^′ is the pressure perturbation defined as the deviation from the hydrostatic pressure. It is important to note that the hydrostatic balance between $\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z\right)$ and the corresponding background pressure field has been subtracted from Eq. (5). Here, k is the unit vector in the positive direction, ν and κ are the constant kinematic viscosity and mass diffusivity, and g is the gravitational constant.

At the deep-water wave-generating boundary condition, Dirichlet time-dependent boundary conditions are enforced as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathbf{u}(x=\mathrm{0},z,t)=f_{\mathbf{u}}(z,t)=\left(\begin{array}{c}{f}_u(z,t)\\ f_w(z,t)\end{array}\right),\text{(8)}& {\displaystyle }\mathit{\rho }(x=\mathrm{0},z,t)=f_{\mathit{\rho }}(z,t),\end{array}$$

where f_u, f_w and f_ρ are explicitly described in Sect. 3 for the different forcing approaches considered here. A free-slip boundary condition is prescribed for the velocity field along all other boundaries, and the density is subject to a boundary condition of zero diffusive flux:

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial \mathit{\rho }}{\partial n}}=\mathrm{\nabla }\mathit{\rho }\cdot \mathbf{n}=\mathrm{0}.\end{array}$$

2.3 Wave-based dimensionless parameters

The definition of the wave-based Reynolds number, Re_w, is given by

$$\begin{array}{}\text{(10)}& {\mathit{Re}}_{\mathrm{w}}\equiv {\displaystyle \frac{\mathit{\lambda }c}{\mathit{\nu }}},\end{array}$$

where λ is the associated wavelength and $c=\mathit{\omega }/k$ the wave-speed of the prescribed wave. The particular wave-based Reynolds number is independent of the wave amplitude and quantifies the strength of viscous effects during the time required for the wave to propagate a distance of one wavelength λ.

The wave aspect ratio, $\mathit{\lambda }/H$, is a metric used to describe how long the wave is in relation to the water depth. Additionally, in the context of this study, it effectively represents the upstream variation in wave-induced flow fields at the wave-generating deep-water boundary which then determines which variant of an Euler–Lagrange approach must be used.

The Froude number quantifies the strength of non-linear effects within the wave against the restoring effect of buoyancy and is defined as

$$\begin{array}{}\text{(11)}& \mathit{Fr}\equiv {\displaystyle \frac{U_{\mathrm{0}}}{c}},\end{array}$$

where U₀ is the maximum wave-induced horizontal fluid velocity.

2.4 Numerical method

The numerical method used to generate the simulation data sets examined in this paper is a high-order continuous Galerkin numerical method, originally developed for the simulation of non-linear, non-hydrostatic internal waves and turbulence in long computational domains with complex bathymetry. The time discretization is semi-implicit and relies on a third-order stiffly stable scheme (Karniadakis et al., 1991). The spatial discretization is based on the nodal spectral element method. Such a discretization enables robust wave propagation against numerical dispersion and diffusion effects, a highly accurate representation of complex geometries and a flexibility in localized resolution – namely, the across pycnocline. Details on the discretization of the Poisson pressure equation and its Laplacian operator (which are directly applicable to the viscous-term treatment) may be found in Appendix A. More details on the numerical model may be found elsewhere (Diamantopoulos et al., 2022).

Such a temporal discretization leads to a Poisson equation for the pseudo-pressure, $\stackrel{\mathrm{‾}}{p}$, at time level n+1:

$$\begin{array}{}\text{(12)}& {\mathrm{\nabla }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{p}}^{n+\mathrm{1}}=\mathrm{\nabla }\cdot \left(-{\displaystyle \frac{\widehat{\mathbf{u}}}{\mathrm{\Delta }t}}\right),\end{array}$$

where $\stackrel{\mathrm{‾}}{p}$ is defined as

$$\begin{array}{}\text{(13)}& \underset{t_n}{\overset{t_{n+\mathrm{1}}}{\int }}\mathrm{\nabla }{p}^{\prime }\mathrm{d}t=\mathrm{\Delta }t\mathrm{\nabla }{\stackrel{\mathrm{‾}}{p}}^{n+\mathrm{1}}.\end{array}$$

For the temporal discretization used in this work, the appropriate boundary conditions for the pressure Poisson equation are taken from Karniadakis et al. (1991). The particular boundary condition is augmented by a term which accounts for the time dependence of the wave-generating boundary condition the boundary-normal velocity:

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{p}}{\partial n}}=\\ {\displaystyle }& {\displaystyle }\mathbf{n}\cdot \left[{\sum }_{q=\mathrm{0}}^{J_e-\mathrm{1}}{\mathit{\beta }}_q\mathbf{N}\left({\mathbf{u}}^{\mathbf{n}-\mathbf{q}}\right)+\mathit{\nu }{\mathit{\beta }}_q\mathbf{L}\left({\mathbf{u}}^{\mathbf{n}-\mathbf{q}}\right)+{\displaystyle \frac{\partial f_{\mathbf{u}}(z,t)}{\partial t}}\right],\end{array}\end{array}$$

where coefficient β_q correspond to a third-order stiffly stable scheme and N and L are the non-linear and linear operators, respectively (Karniadakis et al., 1991).

3 Deep-water wave-generating boundary conditions

The generation of finite-amplitude periodic internal waves is a key component of this study. To this end, we examine the spatio-temporal structure of the generated wave as prescribed by linear theory, which is introduced into the computational domain in the form of time-dependent, vertically variable Dirichlet conditions at the deep-water boundary.

3.1 Internal-wave vertical structure: mathematical descriptions

The fluid's top and bottom boundaries naturally confine the propagation of internal waves so that it occurs in the horizontal direction, along a waveguide formed by the naturally occurring density stratification as shown in Fig. 2. This density stratification, which only varies in the z direction, is often dominated by a region of rapid change (the so-called pycnocline), and it is the up and down motion of this pycnocline that is essential for an accurate description of wave motion.

Figure 2 Schematic of the generation of a mode-1 wave in the waveguide analogy using time-dependent boundary conditions for a two-layer stratification, highlighting the pycnocline displacement, η.

Mathematically, internal waves can be represented as a separation of variable solutions, with a fixed or standing wave structure in the vertical and a propagating waveform (a plane wave in the linear theory) in the horizontal. If we choose the propagation direction to be from left to right and assume the waves to be periodic in x and t, the vertical component of velocity will have the form

$$\begin{array}{}\text{(15)}& w(\mathit{x},t)=W\left(z\right)\exp \left[i(kx-\mathit{\omega }t)\right].\end{array}$$

The equation governing the vertical structure may be derived by linearizing the stratified Euler equations under the Boussinesq approximation (which in turn result from dropping the viscous/diffusive terms in Eq. 4), performing a series of algebraic manipulations to leave an equation for w only and introducing the wave ansatz above. The buoyancy frequency profile, N(z), is assumed to be given, and for the scope of this paper, we are neglecting any form of background shear current. W(z) then becomes the solution of the following linear eigenvalue problem (Gerkema and Zimmerman, 2008) for either ω or k, with the other parameter assumed to be specified:

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}W}{\mathrm{d}{z}^{\mathrm{2}}}}+k^{\mathrm{2}}{\displaystyle \frac{N^{\mathrm{2}}\left(z\right)-{\mathit{\omega }}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}}}W=\mathrm{0}.\end{array}$$

Both top and bottom boundaries are assumed to be impermeable, such that

$$\begin{array}{}\text{(17)}& W\left(\mathrm{0}\right)=W(-H)=\mathrm{0}.\end{array}$$

For a given wave number k, an infinite number of eigenfunctions W_n(z) with their corresponding eigenvalue ω_n exist, each one representing a different vertical mode (i.e., mode-1 does not cross zero in the interior of the fluid, mode-2 crosses zero once in the interior of the fluid, etc.). Therefore, the general solution w(x,t) can be represented by a superposition of such modes, using an arbitrary constant, a_n∈ℂ:

$$\begin{array}{}\text{(18)}& w(\mathit{x},t)=\sum _n{W}_n\left(z\right)\left[a_n\exp \left(i\right(kx-{\mathit{\omega }}_{n}t\left)\right)\right].\end{array}$$

Details regarding the derivation of the solution in the linear stratification case, which are pertinent to the discussion in Sect. 4, are provided in Appendix B.

3.2 Eulerian approach

The above description computes w as a function of a fixed coordinate system. This is often called the “lab frame”, and the theory is labeled as Eulerian (Kundu et al., 2015). In this first approach to generate a finite-amplitude periodic IW, a two-dimensional perturbation field $(u_{\mathrm{E}},w_{\mathrm{E}},{\mathit{\rho }}_{\mathrm{E}}^{\prime })$ is constructed from the solution of the eigenvalue problem for W(z) via the following set of manipulations of the linearized, stratified Euler equations under the Boussinesq approximation:

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\partial {\mathbf{u}}_{\mathbf{E}}}{\partial t}}=-{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathit{\rho }}_{\mathrm{E}}^{\prime }\mathit{k}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{p}^{\prime },\text{(20)}& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{E}}^{\prime }}{\partial t}}=-w_{\mathrm{E}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\rho }}}{\partial z}},\text{(21)}& {\displaystyle }\mathrm{\nabla }\cdot {\mathbf{u}}_{\mathrm{E}}=\mathrm{0}.\end{array}$$

As mentioned above, such an approach can be considered Eulerian since we are looking at the evolution in time of the wave-induced velocity and density fields from a fixed frame of reference. Without a loss of generality, since the chosen equations are linear, only a mode-1 wave will be considered, corresponding to the smallest wave number possible. Following the stratified waveguide analogy (see Sect. 3.1) and multiplying the result by an arbitrary scaling factor A, which is effectively a measure of wave amplitude, the resulting w_E perturbation is

$$\begin{array}{}\text{(22)}& w_{\mathrm{E}}(x,z,t)=-Ak\cos (kx-\mathit{\omega }t)W\left(z\right).\end{array}$$

Using continuity, Eq. (21), an expression for u_E is derived accordingly:

$$\begin{array}{}\text{(23)}& u_{\mathrm{E}}(x,z,t)=A\sin (kx-\mathit{\omega }t){\displaystyle \frac{\mathrm{d}W}{\mathrm{d}z}}.\end{array}$$

Further appealing to the linearized form of the advection–diffusion equation, Eq. (20), the density perturbation, ${\mathit{\rho }}_{\mathrm{E}}^{\prime }$, is then

$$\begin{array}{}\text{(24)}& {\mathit{\rho }}_{\mathrm{E}}^{\prime }(x,z,t)=-{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{\mathit{\rho }}}{\mathrm{d}z}}{\displaystyle \frac{Ak}{\mathit{\omega }}}\sin (kx-\mathit{\omega }t)W.\end{array}$$

The result of the above derivation is a field $(u_{\mathrm{E}},w_{\mathrm{E}},{\mathit{\rho }}_{\mathrm{E}}^{\prime })$, shown in Fig. 3a and b at an arbitrary time of a propagating internal-wave solution of the linear Euler equations under the Boussinesq approximation. The approximate fields are two-dimensional in space, x,z, and also depend on time t. They exhibit a separable structure in x−t and the vertical direction, z.

Figure 3 Snapshots of the velocity fields over one wavelength in the Eulerian approach (a, b) and in the Euler–Lagrange approach (c, d) at an arbitrary time t=0 to highlight the spatial waveform. Panel (e) presents the corresponding wave vertical structures at different locations represented as vertical lines in panels (a) and (c). The displaced pycnocline is represented as a black line. Velocities are normalized with their maximum values.

In practice, the approximations are implemented through a Dirichlet boundary condition (see Eqs. 7 and 8) along a vertical boundary, which we assume to occur at x=0, without any loss of generality.

As a result, the deep-water boundary forcing functions for the Eulerian approach are defined as follows:

$$\begin{array}{}\text{(25)}& {\displaystyle }{f}_u^{\mathrm{E}}(z,t)=A\sin (-\mathit{\omega }t){\displaystyle \frac{\mathrm{d}W}{\mathrm{d}z}},\text{(26)}& {\displaystyle }{f}_w^{\mathrm{E}}(z,t)=-Ak\cos (-\mathit{\omega }t)W\left(z\right),\text{(27)}& {\displaystyle }{f}_{\mathit{\rho }}^{\mathrm{E}}(z,t)=-{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{\mathit{\rho }}}{\mathrm{d}z}}{\displaystyle \frac{Ak}{\mathit{\omega }}}\sin (-\mathit{\omega }t)W.\end{array}$$

3.3 Euler–Lagrange approach

When the Eulerian approach in Sect. 3.2 is applied to situations with a sharp pycnocline, a consistent error is observed for all but the smallest and shortest waves (see Sect. 3.4.2). This disintegration of the generated wave is due to the fact that the up-and-down motion of the pycnocline is not accounted for in the vertical mode description. A natural way to account for this motion is achieved by introducing the wave-induced displacement of the pycnocline. This is measured by the vertical displacement $\mathit{\eta }(x,z,t)$ of the isopycnals (Fig. 2), or isopycnal displacement. Such an approach is labeled as partially Lagrangian since it follows the vertical displacement of individual fluid parcels through time. When combined with the usual description in the horizontal, it is labeled as an Euler–Lagrange approach (Gear and Grimshaw, 1983; Turkington et al., 1991).

Introducing the vertical displacement $\mathit{\eta }(x,z,t)$ of the isopycnals allows for a different, more natural description of the perturbed stratification. Specifically, neglecting the uniform background density, ρ₀, for the sake of compactness, the wave-induced density field can be alternatively expressed as follows:

$$\begin{array}{}\text{(28)}& \mathit{\rho }(\mathit{x},t)=\stackrel{\mathrm{‾}}{\mathit{\rho }}(z-\mathit{\eta }(x,z,t\left)\right),\end{array}$$

which signifies that the density at a point is the same as that at an appropriate height far upstream. The density perturbation, ${\mathit{\rho }}_{\text{E-L}}^{\prime }$, can also be rewritten as follows:

$$\begin{array}{}\text{(29)}& {\mathit{\rho }}_{\text{E-L}}^{\prime }=\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z-\mathit{\eta }(x,z,t)\right)-\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z\right).\end{array}$$

Taylor-expanding the right hand side shows that ${\mathit{\rho }}_{\text{E-L}}^{\prime }$ is a polynomial in η, with the classical Eulerian description giving only the first term:

$$\begin{array}{}\text{(30)}& {\mathit{\rho }}_{\text{E-L}}^{\prime }\approx -\mathit{\eta }{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{\mathit{\rho }}}{\mathrm{d}z}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}\stackrel{\mathrm{‾}}{\mathit{\rho }}}{\mathrm{d}{z}^{\mathrm{2}}}}{\mathit{\eta }}^{\mathrm{2}}.\end{array}$$

By effectively keeping more terms in the expansion, the efficacy of the deep-water forcing is significantly improved. The first term carries a structure representative of the Eulerian approach in Eq. (24). However, from the second derivative of the background density profile in the second term, we can see that Euler–Lagrange effects can be expected to be important for finite-amplitude waves when the stratification exhibits a sharp pycnocline. Such is exactly the case for the study at hand and the continuous two-layer stratification it uses (see Fig. 1 and Sect. 3.4.1).

At this juncture, it is worth emphasizing that the inclusion of the second term in the Taylor expansion of Eq. (30) still only provides an approximate solution of the governing equations (both the linear and non-linear Euler equations) for the input wave field. The purpose of this study is not to provide an exact solution in this context, particularly in a non-linear sense, as enabled, e.g., by the internal-solitary-wave-generating algorithm of Turkington et al. (1991) as it was implemented in Dunphy et al. (2011). Instead, as is subsequently demonstrated, we are aiming for an approximate solution of the linear Euler equations that will drive the deep-water boundary forcing of finite-amplitude waves in a fully non-linear simulation such that the waves can remain robust both near the source and further along the propagation path when non-linear effects modify their waveform.

For the quasi two-layer continuous stratification case, an order-of-magnitude comparison between the two terms on the right-hand side of Eq. (30) may be obtained if one uses as characteristic density and length scales the density jump, Δρ, across the pycnocline and the pycnocline thickness, δ. Per Eqs. (34) and (40), as outlined in the next two sections, one may further write out the isopycnal displacement function as follows:

$$\begin{array}{}\text{(31)}& \mathit{\eta }(x,z,t)=A{\displaystyle \frac{k}{\mathit{\omega }}}r(x,z,t),\end{array}$$

where A is the previously introduced amplitude factor and $r(x,z,t)$ is a structure function, harmonic in x and t and determined by the eigenfunction W(z) in the vertical, which assumes values in the range $[-\mathrm{1},\mathrm{1}]$. Using the characteristic scales above and Eq. (31), one can show that the ratio of the magnitude of the second term to that of the first one on the right-hand side of Eq. (30) scales as $(A/\mathit{\delta })(k/\mathit{\omega })r(x,z,t)$. Therefore, the strength of the Euler–Lagrange effects becomes important when the parameter $(A/\mathit{\delta })/(k/\mathit{\omega })=A/\left(\mathit{\delta }c\right)$ is sufficiently larger than unity; here, $c=\mathit{\omega }/k$ is the linear-phase speed obtained by solving the eigenvalue problem outlined in Sect. 3.1. The particular condition is satisfied for high wave amplitudes, small pycnocline thickness and slow wave propagation speeds.

Restricting the scope to vertical mode-1 waves if $A/\left(\mathit{\delta }c\right)$ is sufficiently larger than unity; Euler–Lagrange effects also become important when the structure function $r(x,z,t)$ is O(1) over a long enough horizontal length scale. For very long waves, $\mathit{\lambda }/H\gg \mathrm{1}$, this is the case over effectively the entire wavelength, and $r(x,z,t)=r(z,t)$. For a wave with finite horizontal wavelength, λ, which is a finite multiple of the water depth (H), the along-wave variation in $r(x,z,t)$ needs to be retained. These two properties of the horizontal structure of $r(x,z,t)$ are at the crux of the formulations outlined in the next two sections.

Inserting the density perturbation, Eq. (29), into Eq. (20), a definition analogous to the free-surface kinematic boundary condition (Hodges and Street, 1999) arises for the isopycnal displacement, η:

$$\begin{array}{}\text{(32)}& {\displaystyle \frac{D\mathit{\eta }}{Dt}}=w.\end{array}$$

When deriving the actual velocity field using the Euler–Lagrange approach, the stratified waveguide analogy is still valid, but with a time-dependent stratification to account for the wave-induced changes to justify the dynamic nature of the pycnocline. The derivation is the same as above, only replacing the fixed frame of reference by the one tracking the displaced pycnocline. The decomposition presented in Eq. (15) is still valid as the change in the frame of reference considered only involves a vertical translation in W. Therefore, the vertical velocity w_E-L is now actually given by

$$\begin{array}{}\text{(33)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{w}_{\text{E-L}}(x,z,t)& {\displaystyle }=w_{\text{Eul}}(x,z-\mathit{\eta },t)\\ {\displaystyle }& {\displaystyle }=-Ak\cos (kx-\mathit{\omega }t)W(z-\mathit{\eta }).\end{array}\end{array}$$

Note that in Eqs. (29) and (33), η is defined using Eqs. (34) or (40), as elaborated in Sect. 3.3.1 and 3.3.2.

3.3.1 Long waves

As a first approximation, since the waves of interest are long ($\mathit{\lambda }/H\gg \mathrm{1}$) and as outlined in the previous section, the x dependence of the vertical displacement, η, is neglected and η is considered to depend only on the vertical position z and time t. The displacement of the fluid in the vertical direction along the wall, $\mathit{\eta }(x=\mathrm{0},z,t)$, is therefore computed by integrating in time the w component of the velocity field in time. The Eulerian approach discussed previously then gives a good approximation, as a starting point, of the w field (Eq. 22), and η may be derived as follows:

$$\begin{array}{}\text{(34)}& \mathit{\eta }(z,t)=\underset{\mathrm{0}}{\overset{t}{\int }}{w}_{\mathrm{E}}(x=\mathrm{0},z,t)\mathrm{d}t={\displaystyle \frac{Ak}{\mathit{\omega }}}\sin (-\mathit{\omega }t)W\left(z\right).\end{array}$$

Using the continuity equation, Eq. (21), u_E-L is derived accordingly. Specific attention needs to be paid to the z variation in the pycnocline's displacement, η. In this regard, using the chain rule in differentiating W(z−η), one obtains

$$\begin{array}{}\text{(35)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial u_{\text{E-L}}}{\partial x}}& {\displaystyle }=-{\displaystyle \frac{\partial w_{\text{E-L}}}{\partial z}}\\ {\displaystyle }& {\displaystyle }=Ak\cos (kx-\mathit{\omega }t)\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}\right)W^{\prime }(z-\mathit{\eta }),\end{array}\end{array}$$

leading to

$$\begin{array}{}\text{(36)}& u_{\text{E-L}}=A\sin (kx-\mathit{\omega }t)\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}\right)W^{\prime }(z-\mathit{\eta }).\end{array}$$

An extra term depending on the z dependence of η now appears in the expression of u_E-L to account for the movement of the pycnocline in contrast with Eq. (23). Note also that the prime denotes a derivative with respect to the argument of W.

As a result, the deep-water boundary forcing functions for the Euler–Lagrange approach are defined as

$$\begin{array}{}\text{(37)}& {\displaystyle }{f}_u^{\text{E-L}}(z,t)=A\sin (-\mathit{\omega }t)\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}\right)W^{\prime }(z-\mathit{\eta }),\text{(38)}& {\displaystyle }{f}_w^{\text{E-L}}(z,t)=-Ak\cos (-\mathit{\omega }t)W(z-\mathit{\eta }),\text{(39)}& {\displaystyle }{f}_{\mathit{\rho }}^{\text{E-L}}(z,t)=\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z-\mathit{\eta }\right)-\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z\right).\end{array}$$

By comparing the spatial structure of the two approaches (Fig. 3), the main difference resides in the structure of the velocity fields: in the Euler–Lagrange approach, the velocity field's deformation tracks that of the pycnocline (plotted in black in Fig. 3c and d) in contrast to the purely Eulerian case where the velocity field treats the pycnocline's position as constant (Fig. 3a and b). Panel (e) shows how the vertical eigenfunctions computed for the wave-displaced pycnocline at the wave peak and trough are offset from the corresponding eigenfunction computed for the initial undisturbed stratification. In the Euler–Lagrange approach, this vertical wave-induced pycnocline displacement is indeed accounted for. For a more detailed discussion of modified vertical eigenfunctions due to wave-induced displacements of the pycnocline and a more thorough mathematical elaboration thereof, the interested reader is referred to the paper by Lamb (1998).

3.3.2 Waves of finite wavelength

So far, the dependence of η on the along-wave position x has been neglected in a first approximation since the considered waves are assumed to be long. Practically, numerical simulations are often required to consider waves of finite wavelength. In this case, to accurately satisfy the continuity equation, Eq. (21), the x dependence must be accounted for and the displacement, η, is expressed as follows:

$$\begin{array}{}\text{(40)}& {\displaystyle }\begin{array}{rl}{\displaystyle }\mathit{\eta }(x,z,t)& {\displaystyle }={\int }_{\mathrm{0}}^t{w}_{\text{Eul}}(x,z,t)\mathrm{d}t\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{Ak}{\mathit{\omega }}}\sin (kx-\mathit{\omega }t)W\left(z\right).\end{array}\end{array}$$

Integrating Eq. (35) now needs to take into account the x dependence of η, leading to a new expression for the horizontal velocity perturbation:

$$\begin{array}{}\text{(41)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{u}_{\text{E-L}}=& {\displaystyle }-{\displaystyle \frac{\mathit{\omega }}{kW\left(z\right)}}\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}\right)W(z-\mathit{\eta })\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathit{\omega }{W}^{\prime }\left(z\right)}{k{W}^{\mathrm{2}}\left(z\right)}}\mathrm{\Phi }(z-\mathit{\eta })+B,\end{array}\end{array}$$

where Φ is the antiderivative of W and B is an integration constant. The second term on the right-hand side of Eq. (41) appears to be, by evaluation of terms across the height of the domain, orders of magnitude smaller than the first one for a typically used wave and is therefore dropped:

$$\begin{array}{}\text{(42)}& u_{\text{E-L}}=-{\displaystyle \frac{\mathit{\omega }}{kW\left(z\right)}}\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}\right)W(z-\mathit{\eta })+B.\end{array}$$

B can be derived from the fact that the horizontal velocity is zero at the depth of the pycnocline, leading to

$$\begin{array}{}\text{(43)}& {\displaystyle }\begin{array}{rl}{\displaystyle }B=& {\displaystyle \frac{\mathit{\omega }}{kW\left(z_{\mathrm{p}}\right)}}\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}(\mathrm{0},z_{\mathrm{p}},t=\mathrm{0})\right)\\ {\displaystyle }& {\displaystyle }\times W(z_{\mathrm{p}}-\mathit{\eta }(\mathrm{0},z_{\mathrm{p}},t=\mathrm{0}\left)\right).\end{array}\end{array}$$

Since the right-hand side of Eq. (42) is not defined on the boundaries for z=0 and $z=-H$, it can be extended using a continuous linear extension, finally resulting in

$$\begin{array}{}\text{(44)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }{u}_{\text{E-L}}(z,t)=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{\displaystyle \frac{\mathit{\omega }}{kW\left(z\right)}}\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}\right)W(z-\mathit{\eta })+B& \text{for }\mathrm{0}>z>-H\\ {\displaystyle \frac{\mathit{\omega }}{k}}\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}(\mathrm{0},\mathrm{0},t)\right)+B& \text{for }z=\mathrm{0}\\ {\displaystyle \frac{\mathit{\omega }}{k}}\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}(\mathrm{0},-H,t)\right)+B& \text{for }z=-H\end{array}\right..\end{array}\end{array}$$

This approach is hereafter referred to as “optimized Euler–Lagrange”.

Note that, when used to implement deep-water wave-generating boundary conditions, in all final expressions for u_E-L, w_E-L and ${\mathit{\rho }}_{\text{E-L}}^{\prime }$ derived in this section or Sect. 3.3.1, the value of x is set to zero without any loss of generality in a manner similar to the Eulerian approach as discussed in Sect. 3.2.

As a result, the deep-water boundary forcing functions for the optimized Euler–Lagrange approach are defined as follows:

$$\begin{array}{}\text{(45)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }{f}_u^{\text{optimized E-L}}(z,t)=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{\displaystyle \frac{\mathit{\omega }}{kW\left(z\right)}}\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}\right)W(z-\mathit{\eta })+B& \text{for }\mathrm{0}>z>-H\\ {\displaystyle \frac{\mathit{\omega }}{k}}\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}(\mathrm{0},\mathrm{0},t)\right)+B& \text{for }z=\mathrm{0}\\ {\displaystyle \frac{\mathit{\omega }}{k}}\left(\mathrm{1}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}(\mathrm{0},-H,t)\right)+B& \text{for }z=-H\end{array}\right.,\end{array}\text{(46)}& {\displaystyle }{f}_w^{\text{optimized E-L}}(z,t)=-Ak\cos (-\mathit{\omega }t)W(z-\mathit{\eta }),\text{(47)}& {\displaystyle }{f}_{\mathit{\rho }}^{\text{optimized E-L}}(z,t)=\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z-\mathit{\eta }\right)-\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z\right).\end{array}$$

The associated adjustments introduced in the pressure boundary condition due to the presence of a time-dependent boundary-normal velocity field at the deep-water boundary are discussed in Sect. 2.4.

Finally, per the previous discussion of Eq. (30), we emphasize that, by construction, neither of the outlined variants of the Euler–Lagrange approach are designed as exact solutions to the linearized Euler equations: perfectly shaped monochromatic waves should not be expected. As illustrated by the results in Sect. 3.4.3, the waves generated through the most suitable of the two Euler–Lagrange approaches (dictated by the aspect ratio $\mathit{\lambda }/H$ at hand) are far more robust than those produced by the purely Eulerian approach. As a result, higher-amplitude longer waves can be used as forcing of fully non-linear simulations (see Sect. 3.4.4).

3.4 Simulations of periodic internal waves in uniform-depth water

3.4.1 Numerical setup

Across all numerical simulations conducted in this study, the wave-based Reynolds number is held constant at ${\mathit{Re}}_{\mathrm{w}}=\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{5}}$. Such a value of Re_w is representative of the laboratory scale yet is sufficiently high to avoid any attenuation in wave amplitude in the propagation zone. The Schmidt number, $\mathit{Sc}=\mathit{\nu }/\mathit{\kappa }$, is fixed at unity. A wave-based Froude number of Fr=0.2 is linked to generated waves that may confidently be characterized as finite-amplitude and can support the development of sufficiently strong non-linear effects as they propagate away from the forcing boundary.

The quasi two-layer continuous stratification profile for $\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z\right)$, given by Eq. (3), is kept the same across all runs. A relatively thin pycnocline with $\mathit{\delta }/H=\mathrm{0.09}$ is used, with a non-dimensional density jump across the pycnocline of $\mathrm{\Delta }\mathit{\rho }/{\mathit{\rho }}_{\mathrm{0}}=\mathrm{2}{N}_{\mathrm{0}}^{\mathrm{2}}\mathit{\delta }/g=\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{3}}$ located at relative position $z_{\mathrm{p}}/H=-\mathrm{0.4}$. The particular value of $\mathit{\delta }/H$ is chosen to mimic the thinner pycnocline of the early fall stratification profile in a long deep lake (Schweitzer, 2010). The deep-water-generated wave used in these simulations is chosen to have an aspect ratio of ${\mathit{\lambda }}_x/H=\mathrm{10.12}$. Such a value of ${\mathit{\lambda }}_x/H$ qualifies the wave as finite-length, albeit not short. Finally, for the particular thin-pycnocline stratification profile and choice of λ_x, the amplitude coefficient A leading to a value of Fr=0.2 corresponds to a value of $A/\left(\mathit{\delta }c\right)=\mathrm{5}$. Euler–Lagrange effects will clearly be present. The choice of ${\mathit{\lambda }}_x/H$ further motivates the question as to whether the fully optimized Euler–Lagrange approach is needed.

All simulations are performed in a uniform-depth tank of depth H and length L=10λ_x. The domain is chosen to be sufficiently long to allow for the development of non-linear effects within the generated waves. Uniformly sized rectangular spectral elements with 224 points per wavelength λ_x are employed in the horizontal direction, whereas 161 points span the entire water column in the vertical direction. The resolutions are given in Table 1, and the elements are uniformly spaced in both length and height. The internal grid point distribution in each element consists of non-uniformly distributed two-dimensional Gauss–Lobatto–Legendre (GLL) integration points (Canuto et al., 2007).

Table 1 Grid point count and resolution for the two-dimensional simulations in a uniform-depth tank.

3.4.2 Limitations of the Eulerian approach

The limitations of the Eulerian approach for the wave forcing are visible in a linear INSE solver (not shown here) but are more readily demonstrated in the framework of a non-linear solver of this type. As shown in Fig. 4a, the Eulerian-generated u velocity field – namely, the shear layer between the upper and lower lobes of opposite velocity tracks horizontally along the location of the undisturbed pycnocline (similar to the top-left panel of Fig. 3) and does not follow the actual pycnocline location (see bottom-left panel of Fig. 3). Immediately visible non-physical numerical features emerge near the forcing boundary at near grid scale, as evidenced by the lobes of alternating sign in the vertical velocity in that region. These spurious vertical velocities are a factor of 2 larger than the theoretically prescribed ones within the target wave, leading to a commensurate reduction in the time step by virtue of the Courant–Friedrichs–Lewy (CFL) condition. Finally, non-negligible regions with density values that exceed the bounds of the background stratification by a factor of 2 to 2.5 are observed (see the blanked-out regions in Fig. 4c). These spurious numerical effects intensify as more waves are generated for the value Fr=0.2. Further intensification of these effects is observed at Fr=0.5 (not shown). In this case, the non-physical non-linear interactions are strong enough to further amplify the near-source spurious vertical velocities and cause an aggressive and prohibitive reduction in the time step.

Figure 4 Fully non-linear simulation using Eulerian wave-generation boundary conditions at $t/T=\mathrm{5.8}$ for Fr=0.2 and $\mathit{\lambda }/H=\mathrm{10}$. The initial undisturbed pycnocline location is represented as a black line. Velocities (a, b) are normalized with their maximum values and the adjusted density (c) by the density jump at the pycnocline. Near-source near-grid-scale lobes of vertical velocity are a factor of 2 to 2.5 larger than that of the prescribed wave. White regions in the density contours correspond to values exceeding the color bar limits which are set by the undisturbed background density profile.

Numerical experimentation indicates that for cases with a well-defined pycnocline, the Eulerian approach produces robust waves and is effectively only valid for small-amplitude (Fr≲0.05, ${\mathit{\eta }}_{\text{max}}/H\lesssim \mathrm{0.02}$) and short-wavelength waves ($\mathit{\lambda }/H\lesssim \mathrm{5}$). It is important to note that the ad hoc linearizing that leads to the Eulerian approach is effectively a linear truncation of a Taylor series. In the amplitude tending to the zero limit, the Eulerian and Euler–Lagrange approaches match. However, even at moderate amplitudes, a significant mismatch is observed (see dotted black lines in Fig. 4). Our approach retains the notation of the Eulerian approach, which is easier to implement in a software setting but effectively introduces the higher-order terms in the Taylor series (at the cost of some algebra).

3.4.3 Linearized Navier–Stokes simulations

To differentiate the features strictly resulting from the differences in the wave-generation approach from the ones resulting from the non-linear effects downstream of the source, the first set of simulations are restricted to solving the linearized incompressible Navier–Stokes equations under the Boussinesq approximation. The non-linear terms have been dropped in Eqs. (4), (5) and (6) analogously to what has been done in the linearized Euler equations, i.e., Eqs. (19), (20) and (21). Per the discussion of the previous section on the limitations of the Euler approach, only the Euler–Lagrange approach is considered here.

The sensitivity to including the x dependence of the isopycnal displacement, η, in the Euler–Lagrange approach is assessed in Fig. 5 by examining the velocity and density fields that are produced by the approaches outlined in Sect. 3.3.1 and 3.3.2 per the corresponding expressions for $f_u^{\text{E-L}}$ and $f_u^{\text{optimized E-L}}$ in Eqs. (37) and (45). The forcing functions f_w and f_ρ have the same structure in both cases and are given by Eqs. (38) and (39), noting any adjustments for the x dependence of η in the optimized Euler–Lagrange formulation. Results are shown after approximately seven wave periods since the initiation of deep-water boundary wave forcing.

Figure 5 Comparison of the velocity and density structures generated by the two different Euler–Lagrange wave-generation approaches in a fully linear simulation at $t/T=\mathrm{7.4}$. The left panels, panels (a) and (c), use the Euler–Lagrange approach and the right panels, panels (b) and (d), use the optimized one. The bottom panels, panels (c) and (d), correspond to an enhanced view of the respective top panels, panels (a) and (b). Velocities are normalized with their maximum values and the adjusted density by the density jump at the pycnocline.

Close to the wave source, periodic shorter-wavelength features are observed for both approaches. These smaller-scale oscillations result from neither of the Euler–Lagrange approaches being an exact solution of the linearized Euler equations as discussed in Sect. 3.3. The amplitude and downstream persistence, however, of these shorter-wavelength effects is markedly weaker in the optimized Euler–Lagrange approach (right panels in Fig. 5), because of the finite wavelength of the generated wave (see Sect. 3.3). To this end, the optimized Euler–Lagrange approach is the method of choice in the fully non-linear simulations given that our baseline wave has an aspect ratio of $\mathit{\lambda }/H$: using it minimizes any possible non-linear interactions between the above parasitic smaller-scale waves and the main target wave which otherwise pose non-trivial challenges for the robustness of the latter wave.

3.4.4 Fully non-linear simulations

The resulting velocity and density fields, obtained by solving the fully non-linear Navier–Stokes equations under the Boussinesq approximation (Eqs. 4, 5 and 6), forced by the optimized Euler–Lagrange approach are shown in Fig. 6 after 10 wave periods. The spurious numerical features close to the boundary have been found to be significantly weaker (not shown here) as compared to what is shown in Fig. 4, whereas the shear layer of the horizontal velocity field tracks the oscillating pycnocline according to Fig. 3 and does not affect the wave generation. Additionally, no spurious mass generation is observed, with density values restricted within the limits dictated by the background stratification. Finally, near-grid-scale vertical velocity near the source remains very small in magnitude and such that the wave-induced vertical velocity is the only factor controlling the time step, as expected. The non-linear response of the generated wave may now be examined along the propagation path without contamination by spurious non-linear interactions due to small-scale near-source transients.

Figure 6 Fully non-linear simulation using optimized Euler–Lagrange wave-generation boundary conditions at $t/T=\mathrm{10}$. Velocities (a, b) are normalized with their maximum values and the adjusted density (c) by the density jump at the pycnocline.

The structure of the generated waves is indeed visibly modified by non-linearity as they propagate away from their source, with different waveform geometries becoming immediately identifiable as a function of distance from the source. Figure 7 attempts to offer such a waveform classification across three different sub-windows along the propagation path. Figure 7a shows waves of depression that develop close to the source. Since the particular waves have large flat plateaus and narrow troughs, a clear similarity with cnoidal waves (Boyd, 2015) is suggested. During a transitional phase, shown in Fig. 7c, the wave troughs broaden. Further downstream, the waves tend to assume a near-sinusoidal shape with peaks and troughs of comparable width (Fig. 7e).

Figure 7 Exploded view of full density contours at different downstream locations for the waves shown in Fig. 6, illustrating the development of waveforms (a, c, e). Respective streamwise Fourier spectra $\widehat{\mathit{\rho }}$ of ρ(x), computed at $z/H=z_{\mathrm{p}}=-\mathrm{0.4}$ and $t/T=\mathrm{10}$ and normalized by the maximum peak depending on the wavenumber normalized by the prescribed wave number k_x (b, d, f).

A more quantitative description of the different types of observed waveforms is enabled by examining the corresponding along-wave spectral content. One-dimensional spatial fast Fourier transforms (FFTs) of the density $\widehat{\mathit{\rho }}\left(k\right)$ are computed streamwise for each sub-window focused at the depth of the undisturbed pycnocline. Special attention needs to be paid when computing the FFTs for the simulations at hand, since the internal grid point distribution within each spectral is non-uniform (Sect. 3.4.1). A non-uniform FFT algorithm is therefore used (Dutt and Rokhlin, 1993; Potter et al., 2017), as it is well tested and readily available.

Closer examination of the right column of panels in Fig. 7, suggests that the along-wave spectral content has power spectral density in regions not specified by the forcing. In particular, a strong second harmonic persists at a downstream distance as large as 60H. Further downstream of the forcing boundary, the amplitude of this harmonic significantly attenuates, resulting in a wave that is closer to being monochromatic (as confirmed by the visualization of Fig. 7e). In the context of a fully non-linear simulation with a sloping boundary, adjusting the length of the section of the computational domain over which the waves propagate prior to reaching the slope allows one to decide how much the waves are allowed to naturally adjust due to their finite amplitude and dispersion. Equivalent simulations, which separate the slope from the source by only a fraction of the horizontal wavelength (Masunaga et al., 2015, 2016), are not expected to support a non-linearly adjusted (and potentially steepened) waveform as the incident wave reaches a slope.

The vertically integrated kinetic energy (KE) at any down-stream position x is

$$\begin{array}{}\text{(48)}& \text{KE}=\underset{-H}{\overset{\mathrm{0}}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{0}}(u^{\mathrm{2}}+w^{\mathrm{2}})\mathrm{d}z.\end{array}$$

KE is shown as a function of time and position for both the fully linear case (Fig. 8a) and non-linear case (Fig. 8b). The theoretically prescribed characteristic of energy transport given by the wave speed, $c={\mathit{\lambda }}_x/T$, is also plotted. Figure 9 presents the interpolation of the kinetic energy along the prescribed characteristic of energy transport in the linear case, shown in white in Fig. 8a. The slope of the KE contours, a measure of the group velocity, appears to match well with the theoretical value. Viscous decay can be considered negligible since the wave-based Reynolds number is chosen to be ${\mathit{Re}}_{\mathrm{w}}=({\mathit{\lambda }}_x^{\mathrm{2}}/\mathit{\nu })/T\gg \mathrm{1}$. The deviation along the characteristic in Fig. 9 can therefore be attributed to the dispersive aspect of the continuous two-layer stratification. In the non-linear case, characteristics also appear to be parallel to the theoretical solution even if the KE is not constant along it due to non-linearities generating extra wavelengths.

Figure 8 Vertically integrated kinetic energy KE normalized by its maximum value KE^* as a function of the downstream position and time for both the linear and non-linear case. The white line corresponds the theoretical energy transport characteristic, $x=c\times t$.

Figure 9 Kinetic energy along the prescribed characteristic of energy transport in the linear case shown in Fig. 8a.

4 Discussion

The generation process of finite-amplitude periodic waves through a time-dependent deep-water Dirichlet boundary condition has been examined for the case of a quasi two-layer continuous stratification (Fig. 1). In the case of a linear stratification, N² is constant in time and space and therefore Eq. (16) has an analytical solution, W_n. Such a solution leads to an exact expression of the perturbation fields for the linear Euler equations under the Boussinesq approximation with an analytical vertical structure and explicit time dependence (see Appendix B). The amplitude of the generated waves can non-trivially exceed the limit prescribed by linear theory without any impact on wave robustness as evidenced by the deep-water waves used in the linearly stratified simulations of Winters (2015).

As described in Sect. 3, the quasi two-layer stratification studied here appears to be more complex. In this context, we do not have an analytical expression for W(z,t), leading to the different approximations introduced in this study. These extra layers of approximations will therefore tighten the amplitude limitations of the wave that can be generated. Nevertheless, depending on the chosen stratification, numerical experiments (not shown here) have demonstrated that robust waves of up to Fr=0.5 are achievable using the optimized Euler–Lagrange approach.

Another important feature that has been demonstrated by Fig. 7 is the fact that in the quasi two-layer continuous stratification the forcing produces wave-trains that are non-monochromatic. Equation (44) reveals that the optimized Euler–Lagrange approach results from the multiplication of temporarily oscillating terms, leading to the appearance of multiple harmonic wavelengths in the generated wave. Sufficiently downstream of the source, the strength of these harmonics seems to diminish for the Fr=0.2 waves shown in Fig. 7e, with wave-induced perturbations that may be regarded as assuming a near-sinusoidal waveform. Nonetheless, the wave trains at the same downstream location in our experiments, with Fr=0.5 (not shown here), are found to remain remarkably non-linear, with extremely steep fronts therein.

Per the literature review in the introduction, the only other computational study considering the generation of long finite-amplitude waves in a two-layer stratification with sufficient distance for the waves to develop downstream of the source that the authors are aware of is that of Dauhajre et al. (2021). We suspect that no issues were reported with regard to the deep-water generated waves for two reasons. First, motivated by apparently different objectives than this study, the use of the wave aspect ratio of $\mathit{\lambda }/H=\mathrm{200}$ is remarkably long and, most likely, restrictive if a turbulence-resolving capability (and not a turbulence parameterization) is preferred. Additionally, the inferred normalized pycnocline thickness value of $\mathit{\delta }/H=\mathrm{0.35}$ is more representative of that in the oceanic continental shelf and not of a deep and long seasonally stratified lake, the primary motivator of this study. Most importantly, noting that Ak=W_max, where W_max is the maximum wave-induced vertical velocity, Dauhajre et al. (2021) work with a typical value of $A/\left(\mathit{\delta }c\right)\approx \mathrm{1.8}$. These non-dimensional parameter values along with any reduction in the curvature at the base of the pycnocline through the insertion of a weakly yet non-trivially stratified lower layer may diminish the intensity of any Euler–Lagrange effects per Eq. (30) and the associated discussion. Finally, it is unclear how the numerical dissipation built into the K-profile parameterization (Large et al., 1994) actively used by Dauhajre et al. (2021) may have damped out any near-source short-wavelength initialization transients such as those reported in Sect. 3.4.3.

5 Conclusions

This study has examined the formulation of robust finite-amplitude periodic internal-wave-generating boundary conditions for a non-linear stratification, highlighting extra levels of subtlety compared to the linear stratification case while relying on a higher-order-of-accuracy spectral element method to discretize the governing equations. The commonly used Eulerian approach, which relies on a fixed reference frame, is found to develop non-trivial errors when implemented in simulations with a sharp quasi two-layer continuous stratification and higher-amplitude internal waves with a horizontal-wavelength-to-depth (wave aspect) ratio that is finite albeit not excessively large. This results in errors because the prescribed wave forcing assumes a fixed/unperturbed pycnocline and does not account for the upstream and vertical wave-induced displacement of the pycnocline. This mismatch between fixed wave-forcing and moving pycnocline is shown to scale with the parameter $A/\left(\mathit{\delta }c\right)$, where A is a measure of wave amplitude, δ is the pycnocline thickness and c is the wave propagation speed. Simulations with values of $A/\left(\mathit{\delta }c\right)=\mathrm{5}$ show spurious mass generation near the wave-generating source, with accompanying unphysical near-grid-scale vertical velocities that can detrimentally reduce the computational time step and even prohibitively restrict it for a long enough time and higher values of the wave-induced Froude number, Fr. The minimal numerical dissipation of the spectral element method cannot damp these spurious numerical features.

For values of $A/\left(\mathit{\delta }c\right)$ sufficiently larger than unity, an Euler–Lagrange approach needs to be used instead in the wave generation, which does account for the above pycnocline displacement. Although an exact solution of the linearized Euler equations under the Boussinesq approximation is not actually attained through this approach, the resulting waves are sufficiently robust: they can propagate away from the source; non-linear adjustments of their waveform are possible through leveraging the higher-order-accuracy spectral element scheme.

The findings of this study will serve as a platform to enable a detailed numerical study of internal swash zones (ISZs), which are zones driven by the interaction of long periodic non-linear internal waves with a sloping boundary. Such simulations will aim to investigate the parameter space in two dimensions, which would include the wave Froude number; pycnocline thickness and depth; wave–aspect ratio; slope value; and the role of no-slip vs. free-slip boundary conditions, particularly on the slope. Select two-dimensional studies will operate as the springboard for full-scale three-dimensional turbulence-resolving simulations These larger simulations may invariably be restricted by existing computational resources to wave aspect ratios in the range [20,40]. As our interests are motivated by internal swash zones in seasonally stratified deep lakes, we will use a two-layer continuous stratification with thinner pycnoclines typical of such environments (Schweitzer, 2015). As such, the parameter $A/\left(\mathit{\delta }c\right)$ will be non-trivially larger than unity. To address this region of parameter space, an Euler–Lagrange approach is needed to account for the wave-induced displacement of the isopycnal field in both vertical and horizontal directions. The optimized Euler–Lagrange approach will be used to generate robust high-amplitude deep-water internal waves at values of the Froude number of up to Fr=0.2.

A parallel avenue of future investigation into the findings of this paper may be their translation to experimental internal-wave generators. Horizontally oscillating paddles are reported as limited to significantly short waves with an aspect ratio of around 10 (Ghassemi et al., 2022). The vertically stacked plate/eccentric camshaft structure of Mercier, Gostiaux and co-workers (Mercier et al., 2010; Gostiaux et al., 2006), through its ability to reproduce a baroclinic structure in the vertical, may be the most amenable experimental technique to adopt aspects of the optimized Euler–Lagrange approach presented here.

Appendix A: Spectral element method

Following the work of Diamantopoulos et al. (2022), the Poisson problem (see Eq. 14) can be rewritten as follows:

$$\begin{array}{}\text{(A1)}& -{\mathrm{\nabla }}_{xz}^{\mathrm{2}}p(x,z)=f\end{array}$$

The discretization of the Laplacian in the non-homogeneous directions x and z is presented. Let 𝒱^N⊂H¹(Ω) be a finite subspace where $p,v\in {\mathcal{V}}^N$ are a part of the solution. Accordingly, the weak form of Eq. (A1) under the Galerkin approximation becomes

$$\begin{array}{}\text{(A2)}& \underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }v\cdot \mathrm{\nabla }p\mathrm{d}\mathrm{\Omega }=\underset{\mathrm{\Omega }}{\int }vf\mathrm{d}\mathrm{\Omega }+\underset{\partial \mathrm{\Omega }}{\oint }v\mathrm{\nabla }p\cdot \mathbf{n}\mathrm{d}S,\end{array}$$

where $\mathrm{\nabla }p\cdot \mathbf{n}=\partial p/\partial \mathbf{n}$ is the natural boundary condition (Deville et al., 2002). By defining 𝒱^N as the finite subspace spanned by two-dimensional Lagrangian basis functions up to the order of N, i.e., ${\mathcal{V}}^N=\text{span}\mathit{\left\{}{h}_{\mathrm{1}}\right(x,z),\mathrm{\dots },h_n(x,z\left)\mathit{\right\}}$, p and v are approximated as $p={\sum }_k{p}_k{h}_k$ and $v={\sum }_k{v}_k{h}_k$, where $k\in \mathit{\{}\mathrm{1},\mathrm{\dots },n\mathit{\}}$ is the corresponding index set, and n is the total number of degrees of freedom on the x–z plane. Note that each element has the same polynomial order. Thus, the discretized Eq. (A2) is written in matrix form:

$$\begin{array}{}\text{(A3)}& K\mathbf{p}=M\mathbf{f}\Rightarrow K\mathbf{p}=\mathbf{g},\end{array}$$

where $K_{ij}={\int }_{\mathrm{\Omega }}\mathrm{\nabla }{h}_i\cdot \mathrm{\nabla }{h}_j\mathrm{d}\mathrm{\Omega }$ and $M_{ij}={\int }_{\mathrm{\Omega }}{h}_i{h}_j\mathrm{d}\mathrm{\Omega }$ are the respective entries of the assembled stiffness and mass matrices (Deville et al., 2002), where $i,j\in \mathit{\{}\mathrm{1},\mathrm{\dots },n\mathit{\}}$. Note that the viscous/diffusive equations for the velocity and density field follow the same weak-form-based formulation and discretization.

A non-overlapping domain decomposition (DD) method with iterative substructuring/static condensation is used when solving for the pressure (Karniadakis and Sherwin, 2013). In tandem with a logically Cartesian topology, Eq. (A3) is broken down into a hierarchy of smaller problems with homogeneous Dirichlet boundary conditions for the two levels of the condensation (Karniadakis and Sherwin, 2013; Huismann et al., 2017; Deville et al., 2002). Once the second and last stage of DD is reached, a Schur complement problem on the vertical interfaces, Γ^v, of the subdomains is iteratively solved. In the context of the hierarchy of problems, a subsequent backward sweep ensures the solution on the global computational domain.

A simple strategy is adopted for the numerical solution of the viscous and diffusive parts of the solver. It is during this step of the solver where boundary conditions for the velocity field and the density perturbation are enforced (Diamantopoulos et al., 2021). Following the discretization presented above, the respective Helmholtz matrix, H^u, is given by

$$\begin{array}{}\text{(A4)}& {\mathbf{H}}^u=\mathit{\alpha }\mathbf{K}+\mathbf{M},\end{array}$$

where the time-step coefficient α scales linearly with νΔt or κΔt.

Appendix B: Details on the derivation of the vertical structure in the linear stratification case

In the case of a linear stratification, we have by definition

$$\begin{array}{}\text{(B1)}& N\left(z\right)=N_{\mathrm{0}}.\end{array}$$

Equation (16) becomes a classic second-order linear differential equation, analogous to that of a simple harmonic oscillator. The solution is oscillatory in z:

$$\begin{array}{}\text{(B2)}& W_n\left(z\right)=\sin \left({\displaystyle \frac{n\mathit{\pi }z}{H}}\right)\text{with}{k}_n={\displaystyle \frac{n\mathit{\pi }}{H}}{\left({\displaystyle \frac{{\mathit{\omega }}^{\mathrm{2}}}{N^{\mathrm{2}}-{\mathit{\omega }}^{\mathrm{2}}}}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

Appendix C: Subtleties of the implementation of the time-dependent deep-water boundary conditions

To implement the different time-dependent deep-water boundary condition approaches described in Sect. 3.2 and 3.3, the eigenfunction W(z) and corresponding eigenvalue k need to be calculated from Eqs. (16) and (17). A high-order spectral element method (Diamantopoulos et al., 2022) is used to for this purpose.

The values of the eigenfunction and its vertical derivative on locations offset from the actual grid points, W(z−η) and $W^{\prime }(z-\mathit{\eta })$, are required for either of the Euler–Lagrange approaches (see Sect. 3.3). These values are obtained at each time step through a cubic spline interpolation in the vertical.

Additionally, to reduce transient-driven contamination of the generated deep-water waves and force both velocity components and density perturbation to be zero at the deep-water boundary at time t=0, the amplitude of the boundary forcing is ramped up in time through application of an exponential envelope. The three forcing expressions for u, w and ρ^′ (see Sect. 3.2 and 3.3) are multiplied by an envelope function, f(t), defined by

$$\begin{array}{}\text{(C1)}& f\left(t\right)=\left(\mathrm{1}-\exp \left({\displaystyle \frac{t}{\mathit{\tau }}}\right)\right),\end{array}$$

where τ is a characteristic timescale of the ramp-up constrained by τ≪T and set to $\mathit{\tau }=\frac{\mathrm{5}}{\mathrm{100}}T$ in this study.