Limiting amplitudes of fully nonlinear interfacial tides and solitons

Nonlin. Processes Geophys., 23, 285–305, 2016 www.nonlin-processes-geophys.net/23/285/2016/ doi:10.5194/npg-23-285-2016 © Author(s) 2016. CC Attribution 3.0 License. Limiting amplitudes of fully nonlinear interfacial tides and solitons Borja Aguiar-González1,2 and Theo Gerkema3 1Departamento de Física, Facultad de Ciencias del Mar, Universidad de Las Palmas de Gran Canaria, 35017 Las Palmas, Spain 2NIOZ Royal Netherlands Institute for Sea Research, Department of Ocean Systems Sciences, and Utrecht University, P.O. Box 59, 1790 AB Den Burg, Texel, the Netherlands 3NIOZ Royal Netherlands Institute for Sea Research, Department of Estuarine and Delta Systems, and Utrecht University, P.O. Box 140, 4400 AC Yerseke, the Netherlands Correspondence to: Borja Aguiar-González ([email protected]) Received: 4 January 2016 – Published in Nonlin. Processes Geophys. Discuss.: 19 January 2016 Revised: 22 June 2016 – Accepted: 7 July 2016 – Published: 18 August 2016 Abstract. A new two-fluid layer model consisting of forced rotation-modified Boussinesq equations is derived for studying tidally generated fully nonlinear, weakly nonhydrostatic dispersive interfacial waves. This set is a generalization of the Choi–Camassa equations, extended here with forcing terms and Coriolis effects. The forcing is represented by a horizontally oscillating sill, mimicking a barotropic tidal flow over topography. Solitons are generated by a disintegration of the interfacial tide. Because of strong nonlinearity, solitons may attain a limiting table-shaped form, in accordance with soliton theory. In addition, we use a quasi-linear version of the model (i.e. including barotropic advection but linear in the baroclinic fields) to investigate the role of the initial stages of the internal tide prior to its nonlinear disintegration. Numerical solutions reveal that the internal tide then reaches a limiting amplitude under increasing barotropic forcing. In the fully nonlinear regime, numerical experiments suggest that this limiting amplitude in the underlying internal tide extends to the nonlinear case in that internal solitons formed by a disintegration of the internal tide may not reach their table-shaped form with increased forcing, but appear limited well below that state. 1 Introduction Tidally generated internal solitons are a widespread phenomenon in the oceans and have been observed and studied for decades (see e.g. Apel et al., 2006). They are intrinsically linked to the internal tide, which itself is generated by barotropic tidal flow over topography. As the internal tide steepens, it may split up into groups of internal solitons, which therefore appear at the tidal period. For internal solitons as such, an archetypal model has been the Korteweg–de Vries (KdV) equation, which is based on the assumption of weak nonlinearity and weak nonhydrostaticity. The equation gives prediction for the relation between amplitude, width and phase speed of the solitons, as well as the shape itself. In the KdV equation there is, mathematically speaking, no limit to the amplitude that solitons may reach (although, of course, at some point the underlying assumption of weak nonlinearity would be violated). This behaviour changes fundamentally if a higher-order (i.e. cubic) nonlinear term is included, leading to the so-called extended KdV (eKdV) equation, as discussed in e.g. Helfrich and Melville (2006). This extended version produces qualitatively different solitons: their amplitude is limited (for a given configuration of layers) and they broaden as they reach their maximum amplitude, the so-called “table-top” solitons. This behaviour is confirmed by fully nonlinear soliton models, as derived by Choi and Camassa (1999) and Miyata (1985, 1988) (denoted as the MCC equations for brevity). In this paper, we focus on another limiting factor, which comes into play even before solitons arise, namely in the internal tide itself. In a purely linear system, the amplitude of the internal tide increases linearly with the strength of the barotropic tidal flow. Here we study how this changes if one includes quasi-linear terms, i.e. retaining products of barotropic and baroclinic fields in the advective terms while still ignoring interactions of the baroclinic field with itself. Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union. 286 B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons We demonstrate that a saturation in the amplitude of the internal tide occurs, and increasing the barotropic flow further does not produce a larger internal tide. As a consequence, when one includes the genuinely nonlinear effects, i.e. products of baroclinic terms, resulting solitons may stay well below their formal limiting amplitude, no matter how strong the forcing. To study these effects we derived a set of fully nonlinear, weakly nonhydrostatic model equations, by extending the MCC equations with a barotropic tidal forcing over topography and with Coriolis effects, which have previously been shown to play a key role in soliton generation from internal tides (Gerkema and Zimmerman, 1995). To avoid having to deal with nonlinearities in the barotropic tide itself (which cannot be formally neglected in a fully nonlinear model), we mimic the interfacial wave generation, replacing a barotropic tidal flow over topography with a horizontally oscillating topography. (There is no exact equivalence between the two, but we demonstrate that, for the parameters used here, the difference remains small.) An alternative approach will also be discussed later. The presence of a topography greatly complicates the subsequent handling of the equations, but we demonstrate that the set of equations can be obtained and can be cast in a form amenable to numerical solving. An extension of the MCC theory with Coriolis effects (MCC-f) was already derived by Helfrich (2007), who investigated the decay and return of internal solitary waves with rotation. We focus on the novel aspect of studying the wave evolution and limiting amplitudes of fully nonlinear, weakly nonhydrostatic internal tides and solitons when forcing and rotational effects are added. We denote our extension of the MCC theory as forced-MCC-f(or forced-MCC in the absence of rotation), for brevity. The paper is organized as follows. We derive a new twofluid layer model consisting of a set of forced rotationmodified Boussinesq equations in Sect. 2. We start with the basic equations and assumptions. Then, we scale equations (Sect. 2.1) and vertically integrate them over the layers (Sect. 2.2). Up to this point, the resulting equations are exact but do not form a closed set. The set is closed by making an expansion in a small parameter measuring the strength of nonhydrostaticity (Sect. 2.3). The resulting model turns out to be equivalent to the Choi–Camassa equations plus additional terms which represent the forcing and rotation effects. Prior to discussing the numerical experiments, we address in Sect. 3 some aspects related to the oscillating topography, the governing nondimensional parameters and the parameter values used for the runs. In Sect. 4 we investigate the factors limiting the growth of tidally generated solitons by first examining the generation of quasi-linear internal tides within the parameter space of this study. Next, in Sect. 5 we solve the full set of forced-MCC-fand explore the conditions by which tide-generated fully nonlinear solitons may actually attain a limiting amplitude. The main findings and conclusions are presented in Sect. 6. The numerical methods and schemes are described in Appendix A. The full set of model equations as solved in the code is presented in Appendix B together with its (quasi)- linearized form. In Appendix C we compare, within the parameter space of this study, the case of an oscillating topography with the case of a tidal flow over a topography at rest. 2 Derivation of the forced-MCC-fmodel We start from the continuity and Euler equations and consider a two-fluid layer system (Fig. 1) with a jump in density across the interface and in which each layer is composed of a homogeneous, inviscid, and incompressible fluid; we apply the Boussinesq approximation. We also assume uniformity in one of the horizontal directions, taking ∂/∂y = 0. Hence, the continuity and momentum equations read as ui,x +wi,z =0,(1) ρui,t +uiui,x +wiui,z −f vi= −pi,x ,(2) vi,t +uivi,x +wivi,z +f ui=0,(3) ρwi,t +uiwi,x +wiwi,z= −pi,z −ρig, (4) where ρiis density, (ui,vi,wi) are the velocity components in Cartesian coordinates, piis pressure, gthe gravitational acceleration, fthe Coriolis parameter (f=2sinφ, at latitude φ) and ρthe mean density. The subscript i=1 (i=2) refers to the upper (lower) layer, and a stable stratification, ρ1< ρ2, is assumed. Boundaries are defined at the surface, taken to be a rigid lid, which is located at z=H1, and at the bottom, located at z= −H2+H(x,t). The time dependence of the bottom will later be specified as a horizontal oscillation, mimicking a barotropic tidal flow over topography. The kinematic boundary conditions at the surface and bottom read as w1=0 at z=H1,(5) w2=Ht+Hxu2at z= −H2+H(x,t). (6) At the interface, z=Z(x,t), the boundary conditions are given by the continuity of normal velocity and pressure: wi=Zt+uiZxand p1=p2at z=Z. (7) For later convenience, we write pressure as the sum of hydrostatic and dynamic parts, the latter being denoted by primes: pi=ρ1gH1−ρigz +p0 i(t,x,z). In the horizontal momentum equation, this amounts to replacing pi,x with p0 i,x, whereas the vertical momentum equation (4) gives ρwi,t +uiwi,x +wiwi,z= −p0 i,z. Nonlin. Processes Geophys., 23, 285–305, 2016 www.nonlin-processes-geophys.net/23/285/2016/ B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons 287 0 Upper surface (rigid-lid) z=H1 z=0 Interface: z=Z(x,t) Bottom z=-H2+H(x,t) Upper layer: H1,ρ1 Lower layer: H2,ρ2 x→ Figure 1. The two-fluid layer system for which the forced-MCC-fequations are derived. The horizontal dashed grey line indicates the level z=0, the level at which the interface resides at rest. The second equation in (7), expressing continuity of pressure at the interface, now becomes (p0 1−p0 2)|z=Z=(ρ1−ρ2)gZ. 2.1 Scaling The next step is to bring the equations into an appropriate dimensionless form for which we introduce the following scales. The scale for the undisturbed water depth is taken to be D, and the typical wavelength L. Crucially, we will assume waves to be long, i.e. nonhydrostatic effects to be weak. This will be expressed by the small parameter1,δ= D L21. Since we allow waves to have large amplitudes (i.e. to be strongly nonlinear), we may take horizontal current velocities to scale with c0=(g0D)1/2, where g0is reduced gravity, g0=g (ρ2−ρ1)/ρ. (Notice that the exact linear long-wave phase speed for interfacial waves, cp, is similar to c0, but has H1H2/D instead of Din the square root.) Thus, uand vwill be scaled with c0. As the interfacial displacement is allowed to be large, an appropriate scale of Zis D. The typical scale of wnow follows from the continuity equation as Dc0/L. Finally, the scale of pressure follows from assuming a primary balance between the acceleration terms ρ utand pxin the horizontal momentum equation. In summary, then, we can introduce the following dimensionless variables, indicated by asterisks, x=L x∗, z =D z∗, t =(L/c0) t∗, p0 i=(ρ c2 0) p0∗ i, ui=c0u∗ i,vi=c0v∗ i, wi=(D/L) c0w∗ i.(8) With these scales, the dimensionless continuity and Euler equations yield (for convenience, we drop the asterisks right away) ui,x +wi,z =0,(9) ui,t +uiui,x +wiui,z −µ vi= −p0 i,x,(10) vi,t +uivi,x +wivi,z +µ ui=0,(11) 1In Choi and Camassa (1999) a small parameter was used instead, which relates to ours as δ=2. δwi,t +uiwi,x +wiwi,z= −p0 i,z.(12) Here µis the scaled Coriolis parameter, µ=f L/c0. Furthermore, we introduce the dimensionless quantities ζ,hi, and h via (Z,H1,H2,H ) =D(ζ,h1,h2,h), so that the scaled form of the boundary conditions is w1=0 at z=h1,(13) wi=ζt+uiζxat z=ζ(x,t), (14) p0 2−p0 1=ζat z=ζ(x,t), (15) w2=ht+u2hxat z= −h2+h(x,t). (16) The goal is now to derive a reduced set of equations from Eqs. (9) to (12), in which the boundary conditions (13)–(16) are incorporated by vertical integration, exploiting the smallness of the parameter δ. The procedure is identical to that of Choi and Camassa (1999), but with the additional complications of the Coriolis force, topography, and tidal forcing. 2.2 Vertically integrated equations We vertically integrate the equations over the upper and lower layers and expand them to the orders δ0and δ1to obtain a closed set for the weakly nonhydrostatic equations, following Choi and Camassa (1999). The layer mean f1of a function f1(x,z,t) for the upper layer is defined as f1(x,t) =1 η1 h1 Zζ dzf1(x,z,t), η1=h1−ζ , (17) and for the lower layer as f2(x,t) =1 η2 ζ Z −h2+h dzf2(x,z,t), η2=h2−h+ζ, (18) where ηirepresents the thickness of the layer (depending on the interfacial displacement ζ). Notice that the boundaries of the integral depend on time and space (x) via the interfacial movement ζ(t,x), but also, for the lower layer, via the horizontally oscillating topography2,h(t,x). Before 2For this reason we need to apply the Leibniz integral rule below with respect to xand t. www.nonlin-processes-geophys.net/23/285/2016/ Nonlin. Processes Geophys., 23, 285–305, 2016 288 B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons proceeding, nonlinear terms in the horizontal momentum equations (10) and (11) are rewritten as (u2 i)x+(wiui)zand (uivi)x+(wivi)z, respectively, to facilitate the procedure. After integration of Eqs. (9)–(11) for i=1 and applying the boundary conditions (13)–(15), we obtain the layer-mean equations for the upper layer: η1,t +(η1u1)x=0,(19) (η1u1)t+(η1u1u1)x−µη1v1= −η1p0 1,x,(20) (η1v1)t+(η1u1v1)x+µη1u1=0.(21) For the lower layer one proceeds similarly, except that now both boundaries are variable. Applying the boundary conditions (14)–(16), vertical integration of Eqs. (9)–(11) for i=2 yields η2,t +(η2u2)x=0,(22) (η2u2)t+(η2u2u2)x−µη2v2= −η2p0 2,x,(23) (η2v2)t+(η2u2v2)x+µη2u2=0.(24) 2.3 Expansion in δ The six integrated Eqs. (19)–(24) derived so far are exact but do not form a closed set. The variables η1,η2and ζcount as one unknown, but we also have ui,vi,p0 i,x,uiuiand uivi, giving 11 unknowns for 6 equations. To obtain a closed set, the last two expressions will be cast in terms of uiand viby using the vertical momentum equation, expanded in terms of the small parameter δ. Furthermore, continuity of pressure at the interface is used to connect the pressure in the lower and upper layers (i.e. p0 1,x and p0 2,x). All in all, the six equations are thus modified to contain only six unknowns. With this aim, we make a formal expansion of the unknowns for the lowest (δ0) and next (δ) orders, as, for example, fi=fi (0)+δf i (1)+... At the lowest order (δ0), p0(0), the dynamics is hydrostatic. At the next order (δ), p0(1)brings weakly nonhydrostatic effects into the system. 2.3.1 Lowest order At lowest order, the vertical momentum equation (12) reduces to ∂p0 i (0)/∂z =0 as terms of order δare neglected; therefore, (perturbation) pressure is vertically constant in each layer. For convenience, we introduce P=p0 2 (0), being a function of tand x. It then follows from continuity of pressure at the interface that p0 1 (0)=P−ζ. Thus, to this order of approximation, p0 1,x =Px−ζx+O(δ), (25) and, for the lower layer, p0 2,x =Px+O(δ). (26) Given the zindependence of pressure and returning to the original horizontal momentum equations, it is now natural to assume that the horizontal velocities, too, are independent of zwithin each layer: uiui=u2 i+O(δ), uivi=uivi+O(δ). At lowest order, then, the set of integrated equations is closed; together with the (exact) integrated continuity equations (19) and (22), we have the momentum equations in terms of the six variables ui,vi,ζand P: (η1u1)t+(η1u2 1)x−µη1v1= −η1(Px−ζx)+O(δ), (27) (η2u2)t+(η2u2 2)x−µη2v2= −η2Px+O(δ), (28) (η1v1)t+(η1u1v1)x+µη1u1=O(δ), (29) (η2v2)t+(η2u2v2)x+µη2u2=O(δ). (30) Recall that η1,2can be expressed in terms of ζand thus involve just one unknown. 2.3.2 Next order At order δ, the procedure is to close the set of six vertically integrated equations by deriving expressions for the horizontal pressure gradients p0 i,x as well as for the contributions of uiuiand uiviin the nonlinear terms. The latter problem is particularly simple. At order δ, the products contain one lowest-order field, which is independent of z(e.g. ui(0)=ui(0)); hence, uiui=1 ηiZdzu2 i=1 ηiZdzui(0)2+2δui(0)ui(1)+... =ui(0)2+2δui(0)ui(1)+... =u2 i+O(δ2), so that uiui=u2 i+O(δ2), uivi=uivi+O(δ2). This allows us to write the horizontal momentum equations as (η1u1)t+(η1u2 1)x−µη1v1= −η1(p0 1 (0)+δp0 1 (1))x+O(δ2), (31) (η2u2)t+(η2u2 2)x−µη2v2= −η2(p0 2 (0)+δp0 2 (1))x+O(δ2), (32) (η1v1)t+(η1u1v1)x+µη1u1=O(δ2), (33) (η2v2)t+(η2u2v2)x+µη2u2=O(δ2). (34) The remaining problem is to find an expression for p0 i (1). At order δ, Eq. (12) reads, in terms of the lowest-order vertical velocities, wi(0)t+ui(0)wi(0)x+wi(0)wi(0)z= −p0 i (1) z.(35) By vertically integrating the continuity equation (9), we obtain an expression for wi(0): wi(0)= −zui,x (0)+ci(t,x), Nonlin. Processes Geophys., 23, 285–305, 2016 www.nonlin-processes-geophys.net/23/285/2016/ B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons 289 where ciare “constants” of integration which are determined by using the boundary conditions at the surface (13) and bottom (16). Thus, wi(0)for the upper and lower layers become, respectively, w1(0)=(h1−z) u1,x (0),(36) w2(0)=(h −h2−z) u2,x (0)+D2h, (37) where the operator Diis defined as ∂/∂t +ui(0)∂/∂x. Substituting w1(0)from Eq. (36) and w2(0)from Eq. (37) into Eq. (35), and vertically integrating the result, we get an expression for p0 1 (1)and p0 2 (1). Taking their derivative with respect to xand their mean over each layer, we finally obtain an expression for p0 i (1) xat order δ. Including the lowest-order terms (25) and (26), this allows us to write the horizontal pressure gradient for the upper layer, p0 1,x =p0 1,x (0)+δp0 1,x (1)+O(δ2) =Px−ζx−δ1 3η1 (η3 1G1)x+O(δ2), (38) and, for the lower layer, p0 2,x =p0 2,x (0)+δp0 2,x (1)+O(δ2) =Px−δ1 3η2 (η3 2G2)x+1 2η2G2hx −η2 2(D2 2h)x−ζxD2 2hi+O(δ2), (39) where we introduced for simplicity the term Gi(as in Choi and Camassa, 1999), Gi=ui,xt (0)+ui(0)ui,xx(0)−(ui,x(0))2.(40) With this, the horizontal momentum equations (31) and (32) become (η1u1)t+(η1u2 1)x−µη1v1 = −η1nPx−ζx−δ1 3η1 (η3 1G1)xo+O(δ2)(41) (η2u2)t+(η2u2 2)x−µη2v2 = −η2Px−δ1 3η2 (η3 2G2)x+1 2η2G2hx −η2 2(D2 2h)x−ζxD2 2hio+O(δ2). (42) We have thus obtained a closed set of six dimensionless equations, namely the exact continuity equations (19) and (22), the horizontal momentum equations (41) and (42), as well as (33) and (34); the last four equations involve the weakly nonhydrostatic assumption. The six unknowns are u1,u2,v1,v2,P, and (via η1,2)ζ. Without interfacial forcing and Earth’s rotation, our set of equations correctly reduces to that of Choi and Camassa (1999). We further specify the model by prescribing the oscillating topography, i.e. the forcing to the system, with h=h(X) whereX(x,t) =x−U0cost, (43) and U0is an arbitrary positive constant. We combine the continuity equations (19) and (22) into (η1+η2)t+(η1u1+η2u2)x=0.(44) Given that η1+η2=h1+h2−h, with the two-fluid system depth h1+h2=1, this leads to −ht+(η1u1+η2u2)x=0.(45) If we now substitute the time derivative of the oscillating topography (43), we have (η1u1+η2u2)x=Uhx,(46) with U=U0sin t, (47) which is the velocity of the oscillating topography with amplitude U0, mimicking a barotropic tidal flow. However, the two are not exactly equivalent, since the transformation from one frame of reference to the other involves an acceleration, and is therefore not Galilean. We discuss this aspect further in Appendix C. Equation (46) can be then integrated in x: η1u1+η2u2=Uh +C(t). (48) Far from the sill (i.e. h→0 for x→ ±∞), we impose the flow to be purely baroclinic, so that the left-hand side must be zero, and hence it follows that C(t) =0. Notice that the righthand side is prescribed via the forcing and thus is a known quantity. It allows us to express u2in terms of u1. We can thus combine the horizontal momentum equations (41) and (42), eliminating P, u1,t +u1u1,x +µv1=ζx+1 (1−h) (Uh)t+(η1u2 1+η2u2 2)x −µ(η1v1+η2v2)−η1ζx +δ1−η1 (1−h)"η1G1η1,x +η2 1 3G1,x# +δη2 (1−h) "−η2G2ζx−η2 2 3G2,x +η2G2 2hx+η2 2(D2 2h)x+ζxD2 2h+O(δ2)(49) u2=Uh −η1u1 η2 ,(50) www.nonlin-processes-geophys.net/23/285/2016/ Nonlin. Processes Geophys., 23, 285–305, 2016 290 B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons v1,t +u1v1,x +µu1=O(δ2), (51) v2,t +u2v2,x +µu2=O(δ2), (52) ζt−(h1−ζ)u1,x +u1ζx=0,(53) where the vi-horizontal momentum equations (51) and (52) have been further simplified from Eqs. (33) and (34) by using the continuity equations (19) and (22). Equation (19) has now been expressed in terms of ζfor convenience. The other continuity equation (22) is no longer included since it is already implicitly present via Eq. (50). All in all, we now have five equations for five unknowns (u1,u2,v1,v2and ζ). The numerical methods and schemes used to solve the model are described in Appendix A. The actual form of the model equations as used in the numerical code is presented in Appendix B. In Appendix C we compare, within the parameter space of this study, the case of an oscillating topography with the case of a tidal flow over a topography at rest. Before concluding this section, it is worthwhile noting an alternative approach. Given the assumption of a rigid lid, one could have also taken U=0 in Eq. (48), the topographic motion set to zero, and then prescribed an external barotropic flux via C(t). Imposing a barotropic flux in this manner does not allow for spatial variations of that flux, as would occur with a free surface, for which an additional dynamical equation would be required to solve the barotropic mode. Specification of C(t) is common in fully nonlinear models of this type as, for example, in Lamb (1994) and Vlasenko et al. (2005). However, the choice of an oscillating topography has also proven to be of use in the study of strongly nonlinear interfacial waves. For instance, Grue (2015) recently confirmed findings on the onset of wave train formation observed in experimental measurements by Maxworthy (1979) with a three-dimensional two-layer, fully dispersive and strongly nonlinear interfacial wave model with a time-varying bottom topography. 3 Numerical experiments: preliminary remarks Whilst not designed to represent a specific region of the world oceans, we aim to investigate in a general manner the conditions by which tidally generated solitons may evolve and, eventually, develop limiting amplitudes in ocean-like scenarios. It is then desirable that leading solitons can propagate towards a mature stage before overtaking preceding internal tides; otherwise, although these are form-preserving features, the tracking of their wave properties becomes cumbersome. For this reason, the parameters that we describe in the following were selected to highlight the qualitative features of these nonlinear processes for a broad range of (mimicked) tidal forcing strengths. Although the model is solved and discussed in nondimensional form, we also present the parameter values in dimensional form to put them in an oceanographic context. 3.1 The oscillating topography and the hydraulic state: the Froude number We define the (dimensional) topography analytically following H(X) =HT 1+(x/HL)2,(54) with xbeing the grid positions in space, and HTand HL being the dimensional parameters which set the height and width of a symmetric sill, respectively. In this manner we ensure perfectly smooth second and third derivatives of the dimensionless topography h(X) in the model equations. At this point it is worthwhile recalling that the oscillation of the topography is introduced in dimensionless form as h= h(X) with X(x,t) =x−U0cost, where U0prescribes the strength (velocity amplitude) of the oscillating topography via U=U0sin(t), the mimicked barotropic tidal flow (see Eqs. 43–47). By increasing U0we enhance the forcing via U, which in dimensional form we introduce, respectively, as U0=c0U0and U=c0U. The topographic obstacle (ridge, sill, etc.) is always centered on the xaxis and the length of the xdomain is chosen to be large enough to prevent waves from reaching the boundaries. In all experiments, the topography starts moving to the left at t=0; we start with a system at rest, i.e. U=u1=u2=0 at t=0. The waves are generated near the origin; on the negative (positive) xaxis, waves travel to the left (right). Because the forcing starts asymmetrically, it is expected that wave packets in the front appear rather differently when comparing both sides (negative vs. positive xdomain). These fronts are the transients, which are influenced by the way the experiment is started. A steady solution at both sides of the xaxis is reached after several tidal periods have passed. In this regard, and to avoid transient effects generated at the start of each run, wave properties have been tracked systematically for the third leftward-propagating interfacial wave counting from the front and after nine tidal periods of forcing. To characterize the hydraulic state, we use the Froude number calculated as F r =U0 cp ,(55) where the amplitude of the mimicked tidal flow, U0, is compared to the linear long-wave phase speed for interfacial waves, cp. The strength of U0leads to three different regimes of interfacial wave generation (see e.g. Vlasenko et al., 2005; Da Silva et al., 2015). Accordingly, the hydraulic regime is denoted, hereafter, as subcritical when Fr 1, critical when Fr ≈1, and supercritical when Fr >1. To account for the varying strength of the tidal forcing within a tidal cycle, we introduce the instantaneous Froude number, defined as Fr0=U/cp. Nonlin. Processes Geophys., 23, 285–305, 2016 www.nonlin-processes-geophys.net/23/285/2016/ B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons 291 Importantly, we also use the Froude number in Appendix C to discuss the applicability of our “non-inertial” frame of reference, the oscillating topography, to the ocean case, where the topography is at rest. To this aim we compare the generation of interfacial waves from the (quasi-)linear forced-MCC equations with that from the (quasi-)linear version of the weakly nonlinear model derived in Gerkema (1996), which works with actual tidal motion. Recall that in the quasi-linear case, barotropic advection is retained but baroclinic interactions are neglected. The equations are then still linear with regard to the baroclinic fields, but the coefficients become time-dependent due to barotropic factors (which are prescribed), so that higher harmonics will be generated when the forcing is increased. For clarification, the (quasi-)linearization of the forced-MCC-fequations is presented in Appendix B. Results from this model comparison confirm a near equivalence between both models within the parameter framework of study, which we restrict to 0 < F r < 1.6. This encourages us to discuss our numerical results, henceforth referring to the strength of the topographic oscillation, U0, as the strength of the tidal flow. 3.2 Parameter values We adopt a two-layer system where the total water depth, D, is set to 100 m, with the upper layer always being thinner than the lower layer (H1< H2). The horizontal oscillation of the moving topography is always of semidiurnal frequency. Although the height of the topography varies between runs, its horizontal scale is kept constant and about 20 km (HL=10 km in Eq. 54). Regarding reduced gravity, g0typically ranges from 0.007 m s−2in the Celtic Sea (Gerkema, 1996) to 0.027 m s−2over the Oregon continental shelf (Stanton and Ostrovsky, 1998); we use this range accordingly. In Table 1 the varying parameters are listed. They vary between runs as indicated in bold fonts, one at a time. The theoretical amplitude of the “table-top” soliton predicted from Eq. (3.68) in Choi and Camassa (1999), and beyond which no solitary wave solution exists, is also indicated. In Sect. 4, runs A1, A2 and A3 illustrate the effect of varying stratification via the reduced gravity, g0. Runs A1, B1 and B2 illustrate the effect of varying the topography ratio, ϕT=HT/D, the height of the topography relative to the total water depth. Finally, runs A1, C1 and C2 illustrate the effect of varying the two-fluid layer thickness ratio, γ=H1/H2. Based on the results from the above analyses, we will argue later why in Sect. 5 we focus on a highly stratified regime (g0=0.03 m s−2) for the study of fully nonlinear waves. For convenience, wave properties are scaled as follows. The interfacial displacement, Z, the internal tide amplitude, A, and the soliton amplitude, As, are scaled to the thickness of the upper layer, H1. The soliton phase speed, cs, is scaled to the phase speed of linear long-wave interfacial waves, cp. Horizontal distances along the x-direction and the soliton width, Ls, are scaled to the wavelength of linear longwave interfacial waves, Lp. Finally, we use the scaled Coriolis parameter, µp, which relates to µin Sect. 2.1, following µp=µ/(2π). 4 Numerical experiments: quasi-linear internal tides Tide-generated solitons emerge from nonlinear disintegration of the underlying internal tides and may be, therefore, naturally subjected to the properties of the latter. For this reason, we find it insightful to investigate first the properties of the underlying internal tides, prior to their nonlinear disintegration, within the parameter space of this study. As described in Sect. 3, the quasi-linear case includes advective terms from the interactions between the barotropic and baroclinic flows, while interactions between baroclinic fields, the genuinely nonlinear terms, are still absent. Therefore, higher harmonics are naturally generated when the forcing is increased. The linear case, where all advective terms are absent, is included here to assess potential departures from the quasi-linear case. Accordingly, Fig. 2 presents the internal tide response to the strength of the tidal forcing for runs A1 to C2 (see Table 1). The minimum forcing strength for all cases is U0=5 cm s−1. In subsequent data points, the increase in U0 is 10 cm s−1from U0=10 cm s−1and onwards, reaching up to Fr∼1.5. In the purely linear experiments, the amplitude of the internal tide increases linearly with the barotropic tidal flow strength. However, the quasi-linear internal tide exhibits a limiting amplitude in all runs, when the tidal forcing increases well above F r =1, a feature that seems to have passed unnoticed in earlier studies. For weak forcing (F r  1), the amplitude of the quasi-linear internal tides approaches the linear case; the advective terms then become very small. This pattern indicates that the decisive factor in the amplitude saturation of quasi-linear internal tides lies in the barotropic advection, which is absent in the linear case. Regarding the comparison between runs with different parameters, we find the following. In Fig. 2a, the increase in stratification causes an earlier deviation between the amplitude growth of the quasi-linear and linear cases, hence occurring at a lower Froude number for runs with a higher stratification (cf. runs A1, A2 and A3). The same effect is observed in Fig. 2b when the height of the topography is increased. The higher the topography, the earlier a deviation from the linear case appears in the Froude number space (cf. runs A1, B1 and B2). Finally, no significant differences emerge regarding the rise of the quasi-linear departure in Fig. 2c, where the thickness of the upper layer varies (cf. runs A1, C1 and C2). These results indicate that the deviation from the linear case arises for lower Froude numbers if either the strength of the stratification or the height of the topography increases. www.nonlin-processes-geophys.net/23/285/2016/ Nonlin. Processes Geophys., 23, 285–305, 2016 292 B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.1 0.2 0.3 0.4 0.5 0.6 (a) Fr −A/H1 A1–L A1–QL A2–L A2–QL A3–L A3–QL 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.1 0.2 0.3 0.4 0.5 0.6 (b) Fr A1–L A1–QL B1–L B1–QL B2–L B2–QL 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.1 0.2 0.3 0.4 0.5 0.6 (c) Fr A1–L A1–QL C1–L C1–QL C2–L C2–QL Figure 2. Amplitude of the linear, L, and quasi-linear, QL, internal tide scaled to the thickness of the upper layer vs. the Froude number. Varying parameters between panels are (a) the strength of stratification, g0(runs A1, A2 and A3); (b) the topography ratio, ϕT(runs A1, B1 and B2); and (c) the two-fluid layer thickness ratio, γ(runs A1, C1 and C2). The run time is nine tidal periods. See Table 1 for further details. Table 1. Summary of runs. Varying parameters are the reduced gravity, g0(m s−2), the topography ratio,ϕT, and the two-fluid layer thickness ratio,γ. The theoretical maximum amplitude, Am, as predicted from Eq. (3.68) in Choi and Camassa (1999), is also indicated. Run g0ϕT=HT/D γ =H1/H2−Am/H1H1,H2ρ1,ρ2 (in m s−2) (in m) (kg m−3) A1 0.03 0.4 0.43 0.67 30, 70 1022, 1025.15 A2 0.02 0.4 0.43 0.67 30, 70 1023.05, 1025.15 A3 0.01 0.4 0.43 0.67 30, 70 1024.1, 1025.15 B1 0.03 0.35 0.43 0.67 30, 70 1022, 1025.15 B2 0.03 0.3 0.43 0.67 30, 70 1022, 1025.15 C1 0.03 0.4 0.33 1 25, 75 1022, 1025.15 C2 0.03 0.4 0.25 1.5 20, 80 1022, 1025.15 Although not shown, it is worth mentioning that the wavelength of the quasi-linear tides does not deviate from the linear case in any of the settings of study and is independent of the strength of the tidal forcing (and hence of the Froude number) and of the height of the topography. However, as predicted from linear theory for interfacial waves, an increase in g0or H1(with H1< H2and Dbeing constant) generates longer internal tides. The amplitude saturation described above is further illustrated in Fig. 3 for run A1, where snapshots of leftwardpropagating quasi-linear internal tides are shown for various forcing strengths (see legend). This spatial view shows how the increase in the forcing transforms the wave from a sinusoidal to an asymmetric shape, indicative of the presence of higher harmonics, while the amplitude becomes saturated. These findings raise the question as to whether solitons emerging from a disintegration of the initially quasi-linear internal tides may be subjected to saturation before they reach a limiting “table-top” shape. We examine this question in the next section by focussing on runs A1, B1 and C1, varying the height of the topography and the thickness of the upper layer while preserving a high stratification. The latter allows us to investigate the broadest range of interfacial wave amplitudes, as suggested by Fig. 2a. 5 Numerical experiments: fully nonlinear internal tides and solitons In this section we investigate the conditions by which tidally generated fully nonlinear solitons may attain a limiting amplitude. Special attention is devoted to factors conditioning the growth of fully nonlinear waves as “table-top” solitons. The main question to address is whether the amplitudes of tidally generated solitons may be subjected to limiting amplitudes of the underlying quasi-linear internal tides, as we hypothesized in the previous section, thus qualifying predictions from classical eKdV and MCC theories. 5.1 Tide-generated “table-top” solitons: run A1 In Fig. 4a a spatial overview of leftward-propagating internal tides and solitons is shown after nine tidal periods of run time. The tidal forcing is fairly strong and this leads to the generation of “table-top” solitons in a supercritical regime (Fr =1.13, U0=90 cm s−1). In subsequent panels, a set of snapshots zooms in on the spatial domain of Fig. 4a to highlight the different stages of the nonlinear disintegration of the internal tides. Nonlin. Processes Geophys., 23, 285–305, 2016 www.nonlin-processes-geophys.net/23/285/2016/ B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons 293 −7 −6.5 −6 −5.5 −5 −4.5 −4 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Z/H1 x/Lp Fr =1.51, U0= 120 cm s−1 Fr =1.39, U0= 110 cm s−1 Fr =1.26, U0= 100 cm s−1 Fr =1.13, U0=90cms −1 Fr =1.01, U0=80cms −1 Fr =0.88, U0=70cms −1 Fr =0.76, U0=60cms −1 Fr =0.63, U0=50cms −1 Fr =0.50, U0=40cms −1 Fr =0.38, U0=30cms −1 Fr =0.25, U0=20cms −1 Fr =0.13, U0=10cms −1 Fr =0.06, U0=5cms −1 Figure 3. Snapshots of the interfacial displacement of leftward-propagating quasi-linear internal tides for run A1 (H1=30 m; Lp=35.49 km). The amplitude saturation is evident as the tidal forcing is increased and the flow becomes supercritical (see legend). The run time is nine tidal periods. At a first stage, Fig. 4b, the internal tide splits up into two different groups of rank-ordered solitons: a train of depressions on the leading edge, and a train of elevations, after the former packet, with initially smaller amplitudes. At a later stage, Fig. 4c, the largest elevations have reached the smaller depressions in the train, and three leading solitons at the front present almost equal amplitudes. Previous solitary wave packets, already propagating away from the generation area, are shown in Fig. 4d and e and correspond to preceding disintegrated internal tides. The “table-top” soliton observed at the leading edge of every preceding internal tide emerged in all cases from the first of the three solitons described previously in Fig. 4c. As the leading soliton evolves and reaches its maximum amplitude, it also broadens, as predicted by soliton wave theory (Helfrich and Melville, 2006), in comparison with subsequent solitons of smaller amplitude (Fig. 4d, e). The observed increase in the distance between the “table-top” soliton and subsequent (smaller) solitons also indicates that, as expected from theory, the leading soliton moves (phase speed) faster than solitons in the tail. Because tidally generated solitons are part of the evolving internal tides, z=0 cannot be used as a reference level to compute the amplitude down to the trough of the soliton (see Figs. 1 and 4). Similarly, the soliton width cannot be measured taking z=0 as a reference level. A criterion is required to adopt a suitable reference level for calculating the soliton amplitude, As, and width, Ls. Here we introduce the reference level Za, which for every leftward-propagating soliton indicates where the first spatial derivative of the interfacial displacement, Z, becomes zero. Accordingly, the soliton amplitude, As, is defined as the vertical distance between Zaand the trough of the leading soliton, located at Zb(see e.g. in Fig. 4c–e). The soliton width, Ls, is defined as the horizontal distance between Zcand Zd, located halfway of the vertical distance spanning As(see also e.g. in Fig. 4c–e). Finally, the soliton phase speed, cs, is computed by subtracting the velocity of the (mimicked) tidal flow, U, from the velocity of the soliton embedded within the internal tide. Using the above criteria, Fig. 5 presents the wave evolution of leading solitons under different forcing strengths (see legend) towards a fully developed stage. Contrary to what one might expect, the amplitudes of the leading solitons decrease during their evolution (Fig. 5a). This can be ascribed to their tide-generated nature. At an early stage, the disintegration of the internal tide leads at its front to a large depression, and this subsequently evolves into a mature leading soliton propagating through the tail of the preceding internal tide (see Fig. 4c–e). The soliton reaches its maximum amplitude slightly before the flow becomes critical (Fr =0.88) and attains the “tabletop” form in the supercritical regime when forced with a stronger tidal flow (Fr =1.13). Unexpectedly, when the tidal forcing is increased even further, the soliton width starts to decrease while keeping its maximum amplitude (cf. Fig. 5a and b). This feature is unlike classical eKdV and MCC theories, suggesting that limiting factors related to the forcing may be acting. Generally speaking, we distinguish between two types of solitons regarding their timescales of growth (see Fig. 5a and www.nonlin-processes-geophys.net/23/285/2016/ Nonlin. Processes Geophys., 23, 285–305, 2016 300 B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons tions (Miyata, 1985, 1988; Choi and Camassa, 1999), extended here with forcing terms and Coriolis effects (forcedMCC-f). The focus is on the effects of the forcing, which represents a novelty in the existing literature and provides a closer view to an ocean-like scenario. The mechanism for internal tide generation is represented by a horizontally oscillating sill, mimicking a barotropic tidal flow over topography. Solitons are generated by a disintegration of the internal tide. The application of an oscillating topography is not completely equivalent to the oceanic case of a tidal flow over a topography at rest. For this reason we have restricted our analyses to a parameter space where a semi-equivalence between both forcing systems was demonstrated (Appendix C). This agreement encourages us to conclude that our findings are not an artifact caused by the use of a mimicked barotropic tidal flow. Of course the findings presented here cannot describe the whole variety of the specific oceanic conditions. However, we believe that this study improves our understanding of the generation and evolution of tide-generated solitons. Numerical solutions show that strongly nonlinear tidegenerated solitons attain in some cases a limiting tableshaped form, in agreement with classical soliton theory. However, results also suggest that tide-generated solitons may alternatively be limited by saturation of the underlying quasi-linear internal tide. In the purely linear system the amplitude of the internal tide increases linearly with the strength of the barotropic tidal flow. But in the quasi-linear case, as the forcing becomes stronger, advective terms become stronger too and cannot be neglected. (Again, in the quasi-linear case, barotropic advection is retained, but interactions of the baroclinic field with itself are neglected). As a result, a saturation in the amplitude of the internal tide occurs; a further increase in the tidal flow does not produce a larger internal tide. This effect seems to have passed unnoticed in previous studies, but might be a key factor in the subsequent disintegration of the internal tide into solitons. It implies that when one includes the genuinely nonlinear effects, i.e. products of baroclinic terms, resulting solitons may stay well below their formal limiting amplitude, no matter how strong the forcing. Interestingly, an increase in the tidal forcing above the value that generates table-shaped solitons, or above the value that simply generates solitons attaining an earlier limitation in growth, causes a narrowing and, subsequently, a decrease in amplitude. The upshot is that increasing the tidal forcing above a certain strength does not lead to larger solitons, but, counterintuitively, to smaller ones. Motivated by the above finding, we performed analogous runs using the full set of weakly nonlinear equations derived in Gerkema (1996). Because these equations are built around the framework of the classical KdV theory and Klein– Gordon equations, one would not expect that an amplitude saturation of solitons could occur. Nevertheless, results (not shown here) demonstrate that both the quasi-linear internal tides and weakly nonlinear tide-generated solitons also exhibit a limiting amplitude. Noting that this model works with an actual tidal flow over a topography at rest, it seems reasonable to argue that the limiting factor is inherent to the tidal forcing. This supports the idea that the forced-MCC-f equations represent an insightful tool for the fully nonlinear framework, where tidally generated solitons may attain limiting amplitudes with or without reaching a table-shaped form. Another departure from classical theories is that strongly nonlinear tide-generated solitons may exhibit larger maximum amplitudes than predicted from eKdV and MCC solutions, while soliton phase speeds are always smaller. We attribute these differences to the fact that tide-generated solitons ride on internal tides and, hence, their wave properties are not simply the response to a two-fluid layer system as such, as in eKdV and MCC solitons, but are also subjected to the forcing of the system, to a variable background flow, and to interfacial displacements of the internal tide itself. In this context, numerical results also show that solitons propagate freely from the source only when the tidal flow is small (subcritical flow), while an increase in the tidal forcing (critical and supercritical flow) generates accelerating and decelerating phases of the soliton speed. In relation to the rotational cases, and in agreement with previous studies (Gerkema and Zimmerman, 1995; Gerkema, 1996), numerical results from the forced-MCC-fequations show that when rotation becomes stronger, the dispersive effect of the Coriolis force becomes stronger too and overcomes nonlinearities, thus preventing the internal tide from disintegrating into solitons. Before concluding we must note, reiterating arguments by Ostrovsky and Grue (2003), that fully nonlinear, weakly nonhydrostatic models entail a paradox to the effect that strongly nonlinear solitons appear from a set of equations that have strong nonlinearity but weak dispersion, while the very existence of solitons presumes a balance between the two. In our case, the MCC-type model is used, involving only the lowestorder nonhydrostatic dispersive terms. Despite the small parameter featuring in the nonhydrostatic terms, they may actually become large in practice (i.e. in the numerical runs) if internal wave profiles are steepening, hence contradicting the original assumption. Indeed, there is no guarantee that the higher-order dispersive terms, which were dropped from these equations, would always remain small. A suggestion for future work is, therefore, to check our results against a numerical computation with a fully nonlinear nonhydrostatic set of equations. 7 Data availability The modeling data used by this study are freely available but not otherwise published in any publicly accessible database. The data can nonetheless be provided on request via e-mail to the first author. Nonlin. Processes Geophys., 23, 285–305, 2016 www.nonlin-processes-geophys.net/23/285/2016/ B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons 301 Appendix A: Numerical strategy We define a grid in time and space for discretization of the various derivatives of the system. Then, tn=n1t and xj=j1x are introduced for integer values of n(time step) and j(spatial step), where 1t and 1x are the magnitude of the steps. Timeand spatial-dependent variables are described as e.g. y(tn,xj), at any time and position. Thus, yn jis the value of the variable yat the current time and spatial step, nand j, respectively. Consequently, n+1 represents the “next time step”, and so n−1 is the “previous time step”, and analogously for jin the spatial grid. The various derivatives in the model are discretized with centered difference approximations (Durran, 1999) as follows: yt(tn,xj)b= yn+1 j−yn j 1t ,(A1) yx(tn,xj)b= yn j+1−yn j 1t ,(A2) yxx(tn,xj)b= yn j+1−2yn j+yn j−1 (1x)2,(A3) yxt (tn,xj)b= yn+1 j+1−yn j+1−(yn+1 j−1−yn j−1) 21x1t ,(A4) yxxt(tn,xj)b= yn+1 j+1−yn j+1−2(yn+1 j−yn j)+(yn+1 j−1−yn j−1) (1x)21t .(A5) Initially the system is at rest with horizontal velocities, ui and vi, and displacement of the interface, ζ, being all zero at the first two time levels (n−1, n). The thickness of the upper, h1, and lower, h2, layers, together with the topography, h(X), describes the two-layer system. At the next time step (n+1), we start to move the topography to the right, creating the effect of a tidal motion flowing to the left. For given U, i.e. scaled velocity of moving topography (Eq. 47), and time step, the excursion of the topography is a known quantity which is used to shift (first, second and third) spatial derivatives of h(X) at every new time step. The time derivatives of the vimomentum and continuity equations (51), (52) and (53) are solved numerically using the third-order Adams–Bashforth approximation (Durran, 1999), for which v1,v2and ζat the next time step (n+1), and at all jpositions, are determined in terms of the known quantities at the previous two time steps (n−1, n). However, solving numerically u1from Eq. (49) is not straightforward, as we deal with three different time derivatives of u1accompanied by space–time-dependent coefficients. Thus, after collecting the various time derivatives involving u1on the one side, and all remaining terms on the other side, the horizontal momentum equation of u1takes an expression of the form a u1,t +b u1,xt +c u1,xxt =Y (tn,xj), (A6) where a,band crepresent spatial derivatives of space–timedependent variables (ζ(x,t) and h(x,t)), and Y (tn,xj)represents a collection of known quantities whose values may be dependent on time and/or space. In the remainder, we describe the numerical method to solve this set of partial differential equations. If we treat the time derivative as a collective term on the left-hand side, we can write (au1+b u1,x +c u1,xx)t=Y (tn,xj) +(atu1+btu1,x +ctu1,xx), (A7) which leads us to the introduction of a new variable, U1, which groups coefficients a,b, and cand time derivatives of u1and turns our problem into a numerically solvable expression of the form U1,t =Y (tn,xj)+(atu1+btu1,x +ctu1,xx). (A8) It is important to recall here that Y (tn,xj)and the spatial derivatives of u1are both evaluated at the current time step (n); the time derivatives of a,band c, which involve values of ζat the current (n) and new time step (n+1), have been previously evaluated with Eq. (53) via Adams–Bashforth approximation. This allows us to rewrite the above expression as U1,t =R(tn,xj)(A9) by grouping all known quantities on the right-hand side under the variable R(tn,xj). Next we need to discretize the time derivative of U1, but before doing that, we discretize its spatial derivatives using Eqs. (A2) and (A3), resulting in U1=aj−2cj 21x u1j+−bj 21x +cj (1x)2u1j−1 +bj 21x −cj (1x)2u1j+1, which we rewrite by introducing factors d,eand fas follows: U1j=dju1j+eju1j−1+fju1j+1.(A10) If we now discretize the time derivative of U1and apply Adams–Bashforth, we obtain a numerically solvable expression for U1at the next time step, which reads Un+1 1j=Un 1j+1t 12 23Rn j−16Rn−1 j+5Rn−2 j,(A11) where U1 n+1 jactually includes Un+1 1j=dn+1 jun+1 1j+en+1 jun+1 1j−1+fn+1 jun+1 1j+1.(A12) To close our system, we still need to obtain un+1 1jfor all j terms. To that end, the equation above is more complicated to solve and gives rise to implicit equations, as we not only www.nonlin-processes-geophys.net/23/285/2016/ Nonlin. Processes Geophys., 23, 285–305, 2016 302 B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons have the unknown un+1 1j, but also un+1 1j−1and un+1 1j+1, which come from the mixed second and third derivatives of u1in Eq. (A6). However, this is a well-known problem that can be solved using the tridiagonal matrix algorithm (TDMA), also known as the Thomas algorithm (Logan, 1987). The choice of the space–time steps 1t and 1x is based on two main requirements. Firstly, the resolution in x(1x) must be sufficiently fine to resolve third-derivative terms and ensure that any short, solitary-like waves are properly resolved. Nevertheless, Kelvin–Helmholtz instabilities may not be filtered out. In this type of model, Jo and Choi (2002) found that solitary waves of sufficient amplitude could be unstable at high wavenumbers due to Kelvin–Helmholtz instability. Thus, if the grid resolution is too fine, unstable short waves will emerge near the wave crest and ultimately explode numerically (Jo and Choi, 2002; Helfrich and Melville, 2006; Helfrich and Grimshaw, 2008). In some cases, the instability can be controlled by filtering out wavenumbers above a threshold (W. Choi, personal communication, 2007; cited in Helfrich and Grimshaw, 2008). For our numerical experiments we consider a 1x course enough to prevent the problem. A second condition follows from the requirement of stability. Thus, for a given spatial step one may take the Courant–Friedrichs–Lewy condition for the linearized equations as an indication of the required time step. The criterion implies that 1x/1t should be larger than the phase speed of the wave, taking special care where the advection by the barotropic tidal flow (here mimicked with the moving topography) should be added to the phase speed to apply the criterion properly (Gerkema, 1994). For the simulations we present, it was not needed to filter out wavenumbers above a threshold to control Kelvin– Helmholtz instabilities, as we designed the space–time grid to avoid this problem. However, in some cases, especially in the simulations where the forcing was fairly strong, an additional trick was needed to retain stability around the generation area (Gerkema, 1994). In those cases averages were taken in the vicinity of the top of the sill (around the steepest part of the topography), where the instabilities arose. At one particular point (xj,tn) in space–time, new values of ui,vi and ζwere calculated by taking the average of the old values at xj−1,xjand xj+1, and subsequently in time between tnand tn−1. The disturbance provoked by this procedure was tested and found to be a minor effect only, as it was only applied over the closest region to the top of the topography. Appendix B B1 Forced-MCC-fmodel equations In Appendix A, the numerical scheme used to solve the model is explained using a generic expression (A6) for the uihorizontal momentum equation (49). Here we present the full set of nondimensional forced equations actually used for the numerical solving of the model. The procedure is as follows. Firstly, all terms of the uihorizontal momentum equation (49) are worked out and grouped according to their physical effects (i.e. linear, nonlinear and dispersive effects from the upper and lower layers, and from topography), leaving unknown quantities involving time derivatives of u1on the left-hand side. The resulting expression (B1) resembles (A6), where coefficients a,band cinvolve derivatives of space– time-dependent variables and Y (tn,xj)is represented here by the sum of all terms on the right-hand side, a u1,t +b u1,xt +c u1,xxt = linear +nonlinear +dispersive1+dispersive2 +dispersivetopo +δη2 (1−h) "(η2hx−η2ζx)φx−η2 2 3φxx +φ(η2 2hxx +ζxhx)i,(B1) u2=Uh −η1u1 η2 , (50) v1,t = −µu1−u1v1,x +O(δ2), (51) v2,t = −µu2−u2v2,x +O(δ2), (52) ζt=(h1−ζ)u1,x −u1ζx, (53) with φ=1 η2hhUt+U2hx+(u1−u2)(η1u1,x −u1ζx)+u2Uhxi,(B2) a(ζ,h) =1+δη2 1−h(η2hx−η2ζx)(η1/η2)x −η2 2 3(η1/η2)xx +η1 η2η2 2hxx +ζxhx#,(B3) b(ζ,h) =δ1−η1 1−hη1ζx +δη2 1−hη1 η2 (η2hx−η2ζx) −2η2 2 3(η1/η2)x#,(B4) c(ζ,h) = −δ1−η1 1−hη2 1 3−δη2 (1−h) η1η2 3,(B5) linear =µv1+ζx+1 1−hhhUt+U2hx+u2hti,(B6) Nonlin. Processes Geophys., 23, 285–305, 2016 www.nonlin-processes-geophys.net/23/285/2016/ B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons 303 nonlinear = −u1u1,x +1 1−h(u1−u2)ζt+u1η1u1,x +u2η2u2,x −µ(η1v1+η2v2)−η1ζx,(B7) dispersive1=δ1−η1 1−h h−η1ζxu1u1,xx −(u1,x)2 +η2 1 3(u1u1,xxx −u1,xu1,xx )#,(B8) dispersive2=δη2 1−hh−η2ζxu2u2,xx −(u2,x)2 −η2 2 3(u2u2,xxx −u2,xu2,xx )#,(B9) dispersivetopo =δη2 (1−h) u2hxη2u2,xx +ζxu2,x +η2 2Uthxx +U2hxxx +2Uu2,x hxx +2u2Uhxxx +3u2u2,xhxx +u2 2hxxx +ζxUthx+U2hxx +2u2Uhxx +u2 2hxxi.(B10) B2 Linear and quasi-linear forced-MCC-fmodel equations The quasi-linear forced-MCC-fmodel follows from neglecting the purely nonlinear terms and weakly nonhydrostatic dispersive terms in Eqs. (B1) and (50)–(53). The equations are linear with regard to the baroclinic fields, but the coefficients become time-dependent due to barotropic advection (which is prescribed) and, therefore, higher harmonics will be generated when the forcing is increased. The quasilinear version of the forced-MCC-fequations reads as u1,t =µv1+ζx +1 1−hhhUt+U2hx+u2ht −µ(h1v1+(h2−h)v2)−h1ζx],(B11) u2=Uh −h1u1 h2−h,(B12) v1,t = −µu1,(B13) v2,t = −µu2,(B14) ζt=h1u1,x.(B15) We notice that the linear runs were actually done somewhat indirectly by taking the quasi-linear version of forcedMCC-fequations above, Eqs. (B11)–(B15), and reducing the forcing by a factor of 100 since the quasi-linear terms cannot be removed explicitly in this model setting without also removing the forcing. Afterwards we enhance the amplitude in the plots accordingly. By reducing the forcing, we effectively enter the linear regime. Appendix C: Oscillating topography vs. tidal flow A Galilean transformation involves two frames of reference which move with constant and rectilinear speed with respect to each other. Hence, observations made in one frame can be converted to another, as physical laws are identical. However, our oscillating topography is not an inertial frame since it is accelerated with respect to a situation where the topography is at rest (as in the ocean). It is, therefore, not evident that the results from the two frames are equivalent. We use the generation model of weakly nonlinear, weakly nonhydrostatic interfacial waves derived in Gerkema (1996), which works with tidal motion over a fixed topography, as a benchmark for testing the impact of our “non-inertial” frame of reference. If we compare interfacial waves generated from the nonlinear version of both models, differences are expected to arise from the fact that forced-MCC equations are fully nonlinear. For this reason we restrict the comparison to the linear and quasi-nonlinear model versions. If the results between the forcing systems turn out to be similar, it seems reasonable to assume an equivalence in the nonlinear case as well. In Appendix Fig. C1, interfacial waves generated from both models are presented for various numerical experiments under a fairly strong forcing, i.e. when both models may be expected to deviate more noticeably from each other. Our interests focus then on the upper limit of the supercritical regime (Fr>1) that we can reach while still preserving a good agreement between both generation mechanisms. The different settings in Fig. C1 differ in the strength of stratification from top to bottom panels, while the thicknesses of the upper and lower layers (H1=30 m, H2=70 m) and the height and width of the sill are kept fixed (HT=40 m and HL=10 km in Eq. 54). Results from Fig. C1 indicate that in all cases a close correspondence exists between numerical solutions from Gerkema (1996) (grey line) and the forced-MCC equations (black line), suggesting only a minor impact of the noninertial nature of our frame of reference when reaching up to Fr∼1.5. These results encourage us to interpret the speed of the oscillating topography as the “strength of the tidal flow” within the parameter space of this study. www.nonlin-processes-geophys.net/23/285/2016/ Nonlin. Processes Geophys., 23, 285–305, 2016 304 B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons Figure C1. Linear (left panels) and quasi-linear (right panels) interfacial waves generated via a tidal flow over a sill from the model equations derived in Gerkema (1996) (grey line) and via a horizontally oscillating sill from the model equations derived in this study (black line). The Froude number and corresponding strength of the (mimicked) tidal flow are indicated in the upper-right corner of each panel. For scaling purposes one must note that the wavelength of the linear long-wave interfacial wave, Lp, varies from top to bottom panels as Lp=35.5 km (g0=0.03 m s−1) in (a) and (b);Lp=29 km (g0=0.02 m s−1) in (c) and (d); and Lp=20.5 km (g0=0.01 m s−1) in (e) and (f). The run time is nine tidal periods. Nonlin. Processes Geophys., 23, 285–305, 2016 www.nonlin-processes-geophys.net/23/285/2016/ B. Aguiar-González and T. Gerkema: Limiting amplitudes of fully nonlinear interfacial tides and solitons 305 Acknowledgements. Financial support was provided by the Spanish government (Ministerio de Ciencia e Innovación) through a PhD grant (FPU) awarded to the first author (AP2007-02307), and through the Netherlands Organization for Scientific Research (NWO), section Earth and Life Sciences (ALW), via ZKO grant no. 839.08.431 awarded to the “INdian-ATlantic EXchange in present and past climate” (INATEX) programme. The authors also wish to thank the reviewers for their helpful remarks and suggestions. Edited by: V. Shrira Reviewed by: V. I. Vlasenko and two anonymous referees References Apel, J. R., Ostrovsky, L. A., Stepanyants, Y. A., and Lynch, J. F.: Internal solitons in the ocean, Tech. Rep., Woods Hole Oceanographic Institution, 2006. Choi, W. and Camassa, R.: Fully nonlinear internal waves in a twofluid system, J. Fluid Mech., 396, 1–36, 1999. Da Silva, J. C. B., Buijsman, M. C., and Magalhaes, J. M.: Internal waves on the upstream side of a large sill of the Mascarene Ridge: a comprehensive view of their generation mechanisms and evolution, Deep-Sea Res. Part I, 99, 87–104, 2015. Durran, D.: Numerical methods for wave equations in geophysical fluid dynamics, vol. 32, Springer Verlag, 1999. Gerkema, T.: Nonlinear dispersive internal tides: generation models for a rotating ocean, Ph.D. thesis, Utrecht University, 1994. Gerkema, T.: A unified model for the generation and fission of internal tides in a rotating ocean, J. Mar. Res., 54, 421–450, 1996. Gerkema, T. and Zimmerman, J. T. F.: Generation of nonlinear internal tides and solitary waves, J. Phys. Oceanogr., 25, 1081–1094, 1995. Gerkema, T. and Zimmerman, J. T. F.: An introduction to internal waves: Lecture Notes, Royal NIOZ, Texel, 2008. Grue, J.: Nonlinear interfacial wave formation in three dimensions, J. Fluid Mech., 767, 735–762, 2015. Helfrich, K. and Grimshaw, R.: Nonlinear disintegration of the internal tide, J. Phys. Oceanogr., 38, 686–701, 2008. Helfrich, K. R.: Decay and return of internal solitary waves with rotation, Phys. Fluids., 19, 026601, doi:10.1063/1.2472509, 2007. Helfrich, K. R. and Melville, W. K.: Long nonlinear internal waves, Annu. Rev. Fluid Mech., 38, 395–425, 2006. Jo, T.-C. and Choi, W.: Dynamics of strongly nonlinear internal solitary waves in shallow water, Stud. Appl. Math., 109, 205–227, 2002. Kakutani, T. and Yamasaki, N.: Solitary waves on a two-layer fluid, J. Phys. Soc. Japan, 45, 674–679, 1978. Lamb, K. G.: Numerical experiments of internal wave generation by strong tidal flow across a finite amplitude bank edge, J. Geophys. Res., 99, 843–864, 1994. Logan, J. D.: Applied mathematics: a contemporary approach, Wiley, 1987. Maxworthy, T.: A note on the internal solitary waves produced by tidal flow over a three-dimensional ridge, J. Geophys. Res., 84, 338–346, 1979. Miyata, M.: An internal solitary wave of large amplitude, La mer, 23, 43–48, 1985. Miyata, M.: Long internal waves of large amplitude, in: Proceedings of the IUTAM Symposium on Nonlinear Water Waves, 399–400, 1988. Ostrovsky, L. A. and Grue, J.: Evolution equations for strongly nonlinear internal waves, Phys. Fluids, 15, 2934–2948, 2003. Ostrovsky, L. A. and Stepanyants, Y. A.: Do internal solitions exist in the ocean?, Rev. Geophys., 27, 293–310, 1989. Stanton, T. and Ostrovsky, L.: Observations of highly nonlinear internal solitons over the continental shelf, Geophys. Res. Lett., 25, 2695–2698, 1998. Vlasenko, V., Stashchuk, N., and Hutter, K.: Baroclinic tides: theoretical modeling and observational evidence, Cambridge University Press, 2005. www.nonlin-processes-geophys.net/23/285/2016/ Nonlin. Processes Geophys., 23, 285–305, 2016 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.