Variational Water Wave Modelling: from Continuum to Experiment

Variational methods are investigated asymptotically and numerically to model water waves in tanks with wave generators. As a validation, our modelling results using (dis)continuous Galerkin finite element methods will be compared to a soliton splash event resulting after a sluice gate is removed during a finite time in a long water channel with a contraction at its end.


Introduction
A popular approach in the modelling of nonlinear water waves is to make the approximations that the three-dimensional fluid velocity u is irrotational and divergent-free, such that u = ∇φ and ∇ · u = ∇ 2 φ = 0, and that the dynamics is inviscid, such that the dynamics is governed by variational and Hamiltonian dynamics [10,13]. At least symbolically one can invert this Laplace equation for the interior potential φ and reduce the dynamics to the free surface, expressed in terms of the potential φ s at the free surface and the position of this free surface. For non-overturning waves, this free surface dynamics can be expressed in terms of the water depth h = h(x, y, t) and φ s (x, y, t) = φ(x, y, z = b + h, t) with horizontal coordinates x and y as well as time t. Here the fixed topography is denoted by b = b(x, y). The free surface thus lies at the vertical level z = b(x, y) + h(x, y, t), parametrised by x and y.
One then often considers the initial value problem governed by autonomous Hamiltonian dynamics for h and φ s with initial conditions h(x, y, 0) and φ s (x, y, 0) without any forcing or dissipation. In practical situations, however, waves are generated continuously by wave-makers or temporarily by opening a sluice gate, both involving time dependent internal or boundary conditions. This implies that the dynamics is non-autonomous, including explicit dependence of the equations on time. Sometimes, these non-autonomous aspects can be included in the variational principles governing the wave dynamics.
We will therefore start to formulate finite-dimensional variational dynamics in which the variational principle indeed depends explicitly on time. The forced-dissipative nonlinear pendulum with the harmonic oscillator as linearisation is a first example of such a non-autonomous variational principle. It is expressed in terms of the coordinate or angle q = q(t) and its conjugate momentum p = p(t). We will show that a fairly generic non-autonomous variational principle can be transformed into an autonomous one by lifting it to a larger phase space.
The variational principles for water waves can be derived from the incompressible Euler equations with a dynamic free surface, using both the velocity potential and the pressure as Lagrange multipliers to impose two constraints. In this approach [10,4], the pressure variable eventually disappears from the formulation. Both driven wave makers and the removal of a sluice gate can be included in the variational principle for water waves. It adds explicit time-dependence to the principle. In Bokhove et al. [3] and Thornton et al. [17], this was investigated for water wave dynamics in a vertical Hele-Shaw cell with a moving interface between water and air, concerning damped motion with a wave pump and with damped wave sloshing, respectively. The Hele-Shaw cell is a narrow wave tank with a wave-maker, e.g., it has a width of 2 mm, length of 0.6 to 1 m, and depth of circa 0.1 m, between two closely spaced glass plates. The bulk damping due to the dominant friction against these glass plates leads to an exponential time dependent factor e γt in the variational principle. Its inverse expresses the exponential decay due to friction with γ > 0 proportional to both the viscosity and the square inverse of the gap width between the glass plates (e.g., Bokhove et al. [3]).
Such explicit time dependence makes the numerical discretisation of the water wave problem more complicated and prone to numerical instabilities. In Gagarina et al. [6], this challenge was overcome, verified and validated for driven potential flow water wave dynamics. As a second example and to minimise overlap with that research, we illustrate the numerical modelling of variational water wave dynamics with a reduced water wave model that includes a simple, new model of a sluice gate. By modifying work of Pego & Quintero [15], we derive a model similar to Benney and Luke's [1] but remain entirely within the variational framework from the onset. The velocity potential at the bottom topography Φ and free surface deviation η = h − H 0 herein play a similar role as q and p for the nonlinear pendulum, with z = H 0 the water level at rest. The above two examples allow us to illustrate our time discretisation techniques for non-autonomous water wave dynamics. Finally, we validate the new model discretisation against a soliton splash phenomenon in a water wave channel with a removable sluice gate and a V-shaped channel contraction at its end.
While some techniques have been published before [3,6,7,17], the derivation of the modified Benney-Luke equations, the systematic way to lift non-autonomous to autonomous variational principles in time, and the use of this lifting in finding extended time discretisations are novel. Also new is the soliton splash simulation.

Variational mechanics continuous in time 2.1 General formulation
We will consider non-autonomous Hamiltonian dynamics cast in the variational principle with time t in the time interval t ∈ [0, T ] for some T > 0, generalised coordinate (vector) Q = Q(t) and conjugate momentum (vector) P = P (t), and Hamiltonian H = H(P , Q, t). Note that the Lagrangian density L = L(Q, P , t) depends explicitly on time, in order to model forcing as well as (simplified) damping within the variational principle. The variations in (1) are defined as Hence, the variation in (1) becomes in which end-point contributions cancel because of the conditions δQ(0) = δQ(T ) = 0. The latter follow since Q(0) = Q 0 is a given value and by using time reversibility Q(T ) = Q T can, likewise, be considered as given. Their variations are thus zero. Since variations δQ(t) and δP (t) are arbitrary and pointwise in time, Hamilton's equations follow from (4) as Finally, energy is generally not conserved when H also depends explicitly on time, since

Kamiltonian formulation
The non-autonomous variational principle (1) can be lifted to an autonomous variational principle (cf. Goldstein [8]) by considering time as an auxiliary variable with a new time coordinate. First, we rename time as τ in (1) to obtain Subsequently, we take τ = τ (t) with a new time coordinate t, such that dτ /dt = 1 and τ (0) = 0 for convenience, and introduce a conjugate variable p for τ . A so-called "Kamiltonian" is then defined as K(Q, P , τ, p) ≡ H(Q, P , τ ) + p, and the corresponding variational principle and variations are where we used δQ(0) = δQ(T ) = 0 as well as δτ (0) = δτ (T ) = 0. Hence, the extended Hamilton's or "Kamilton's" equations become We note that this last equation is decoupled from the rest of the equations. Moreover, the Kamiltonian K is by construction conserved in time: dK/dt = 0. The above reformulation into a "Kamiltonian" system indicates that a variational discretisation of an autonomous variational principle can thus, in principle, be extended systematically to one for a non-autonomous system.

Linearisation
Hamilton's equations (5) with ∂H/∂P = P for general (nonlinear) dynamics can be linearised and diagonalised in a series of uncoupled oscillators. This will be used later to interpret the numerical verification.
Autonomous case: Hamilton's equations (5) have fixed points (Q,P ) obeying ∂H(Q,P )/ ∂Q = 0 andP = 0. When we substitute the decomposition Q =Q+Q and P =P +P into the variational principle (1) and keep terms quadratic into the (small-amplitude) perturbation variables Q and P , we obtain a series of coupled linear oscillators. After an eigen-analysis, diagonalisation and under certain assumptions, a system of uncoupled harmonic oscillators remains. The one with the largest eigenfrequency ω = ω 1 with q 1 = Q 1 and p 1 = P 1 will therefore satisfy the following principle which variation yields the classic harmonic oscillator upon using the usual end point conditions. That is dq 1 /dt = p 1 and dp 1 /dt = −ω 2 q 1 .
non-autonomous case: In contrast to the autonomous case, {∂H(Q,P )/∂Q = 0,P = 0, p = 0} is not a fixed point of the Kamiltonian system (9). However, the Kamiltonian form can be Taylor expanded around τ = 0 and the fixed points of the autonomous case, focussing on one oscillator as before. Motivated by our examples, this heuristically yields the following "archetypical case" for a = ∂ 2 V (q 1 , τ )/∂q 1 ∂τ | q1=Q1,τ =0 and H = p 2 1 /2 + V (q 1 , τ ) and by using the transformation ωq 1 = ωq 1 + aτ /ω,p 1 = p 1 in the second line, cf. [7]. Consequently, the Kamiltonian system decouples with a new autonomous Hamiltonian 1 2p as an invariant in time t. This regularises the Kamiltonian to one with at most quadratic terms, including a term linear in p 1 .

Forced-dissipative nonlinear oscillator
As first, low-order example, consider a forced and damped nonlinear oscillator with generalised coordinate q and momentum p, potential V = V (q), time dependent forcing function F = F (t) and damping coefficient γ ≥ 0 (for a damped oscillator it is given in Olver [14]). We assume that the leading order term in a Taylor expansion of V (q) is quadratic in q. The dynamics become The variational principle (13a) can be put in the form (1) by using P = pe γt/2 and Q = qe γt/2 , in which the leading order term 1 2 aQ 2 (for some a > 0) in a Taylor expansion of V (Qe −γt/2 )e γt is now quadratic in Q, while higher order terms contain decaying exponentials in (nonzero) powers of e −γt/2 . As a consequence, when in the unforced case F (t) = 0 or when lim t→∞ F (t)e γt/2 → 0, the new Hamiltonian lim t→∞H = 1 2 P 2 + 1 2 γP Q + 1 2 aQ 2 will become a quadratic invariant at long times. A Kamiltonian form (8) of the above variational principle, with conjugatep for τ , is Another, second transformation is P = pe γt and Q = q.
It yields the alternative Kamiltonian variational principle The above two transformations are important because they show us how to find time discretisations for the non-autonomous case based on ones for the autonomous variational principle.
A linear forced and damped harmonic oscillator emerges for V (q) = 1 2 ω 2 q 2 with F (t) = a sin(Ωt), forcing frequency Ω and amplitude a, as well as natural frequency ω. The system then becomes dq/dt = p and dp/dt = −ω 2 q − γp + a sin(Ωt). It has an analytical solution, comprised of the sum of homogeneous and quasi-steady state particular solutions, for β = ω 2 − 1 4 γ 2 real, initial conditions q(0) = q 0 and p(0) = p 0 , and with 4 Variational water wavesà la Benney-Luke As a second, high-order example, we derive a reduced water wave model and subsequently discretise it in space using a Galerkin finite element expansion.

Principle for finite, shallow depth
The deviation from a still water surface is denoted by η = η(x, y, t) and the velocity potential as φ = φ(x, y, z, t) with horizontal coordinates x and y, vertical coordinate z and time t. The still water rest level lies at z = H 0 and the flat bottom is at z = 0 with vertical velocity ∂ z φ = ∂φ/∂z = 0. The free surface and bottom potentials are defined as φ s (x, y, t) = φ(x, y, H 0 + η(x, y, t), t) and Φ(x, y, t) = φ(x, y, 0, t), respectively. Horizontal gradients are ∇ = (∂ x , ∂ y ) T with transpose (·) T . The domain has horizontal extent Ω h with vertical walls at ∂Ω h , where the normal gradient n · ∇φ = 0.
The potential flow water wave equations can be derived from Luke's variational principle with constant acceleration of gravity g, see Luke [10]. We added a (gravitational) potential η R = η R (x, t) to model a sluice gate release problem. For flow at rest, φ = 0 and η = η R (x). By adjusting η R (x, t) to become independent of x in a finite time, such a release can be modelled approximately in a relatively straightforward manner.
The following transformations (from Pego & Quintero [15]) are first applied with amplitude parameter = α 0 /H 0 1 (for small wave amplitudes α 0 , see also Fig. 1), which transform the variational principle after dropping the hats into Adopting the second set of scalings (from Pego & Quintero [15]) with small dispersion parameter µ = (H 0 / 0 ) 2 1 (for long waves with wavelength 0 , see Fig. 1), gives after dropping the tildes and a constant The system resulting from (26) is We then expand φ in terms of the small parameter µ about the bottom potential, Substituting this expansion into the Laplace equation (27a) (and using ∆ = ∇ 2 ), results in µ∆Φ + (µ 2 ∆Φ 1 + µ∂ zz Φ 1 ) + (µ 3 ∆Φ 2 + µ 2 ∂ zz Φ 2 ) + · · · = 0. In addition, φ| z=0 = Φ by definition and ∂ z φ| z=0 = (µ∂ z Φ 1 + µ 2 ∂ z Φ 2 + · · · )| z=0 = 0 from equation (27b), therefore Φ 1 , Φ 2 and their vertical gradient vanish at z = 0. At leading order in µ, the equation . Therefore up to second order the expansion becomes We use this expansion in the variational principle (26) and retain terms up to order O( 2 µ, 3 ). After integration in z, the resulting leading-order principle becomes in which we used the (somewhat ad hoc) boundary conditions n·∇Φ = 0 and n·∇∆Φ = 0 at ∂Ω h with outward horizontal normal n, in the integration by parts, as well as suitable end-point conditions at t = 0 and t = T . The variation of (30b) with respect to η and Φ in (30c) yields Hamilton's field equations An alternative formulation lowering the highest order derivatives is which formulation is advantageous in a C 0 finite element formulation, used later. The resulting equations are Upon elimination of the auxiliary variable q, we again find (31).

Spatial Galerkin-Ritz finite element approach
The two variabes η and Φ are expanded in terms of compact (finite element) C 0 -basis functions ϕ l (x, y), yielding approximations in which the Einstein summation convention is used with indices k and l running over the finite degrees of freedom. The domain Ω h is tessellated with quadrilaterals. Direct substitution of (34) into the variational principle (30b) leads to extra jump terms, arising from products of terms with second-order spatial derivatives because the basis functions are only continuous. Instead, we substitute (34), as well as into the variational principle (32a), to obtain This is a so-called Galerkin-Ritz approach, with matrices M kl , S kl and vector R l given by Variation of (36) yields For R l = 0, the Hamiltonian of the system is . Several limits are enclosed in (30b) or (32a) and its spatial discretisation (36): nonlinear, potential flow shallow water equations emerge for µ = 0; linear Benney-Luke type equations for = 0 and µ > 0; and, linear, potential flow shallow water equations for = µ = 0. By defining P k = (M kl + µ 2 S kl )η l and Q k = φ k , (36) obtains the standard form (1), which discretisation in time is considered next.

Defining the time derivative across time nodes
Time interval [0, T ] is partitioned into N time intervals I n = [t n , t n+1 = t n + ∆t n ] for n = 0, . . . , N − 1. The variables in the variational principle (1) or (8) are approximated using a discontinuous Galerkin finite element method in time. These approximations Q δ , P δ are continuous in each time interval but discontinuous across. Across each time node t n = n−1 j=0 ∆t j for n > 0 with t 0 = 0 and t N = T , the contribution of the term P · dQ/dt in the variational principle needs to be determined. It leads to a δ-function at each node in the integrand, but since the variational principle is a weak form, a finite sum of contributions emerges. Since the expansions are discontinuous, limiting values on either side of each time node t n are discontinuous, with Q n,− = lim ↓0 Q δ (t n − ) and Q n,+ = lim ↓0 Q δ (t n + ), and likewise for P n,− and P n,+ . The jump in a quantity across the node will be denoted as, e.g., [Q] = Q n,+ − Q n,− and the weighted average as, e.g., {P } = αP n,+ + (1 − α)P n,− with α ∈ [0, 1].

Limit continuous Galerkin to discontinuous Galerkin:
Consider all time nodes n > 0 and split each time node into a narrow element t ∈ [t n+1 − ∆τ, t n+1 + ∆τ ] of width 2∆τ such that none of these narrow elements overlap another. Hence, the original elements shrink slightly in size. Instead of a discontinuous Galerkin finite element method in time, we have created a continuous finite element in time by using expansions within these narrow elements such that the overall expansion becomes continuous. Consider the mapping t = t n+1 − ∆τ + (2∆τ )τ within this narrow element to a reference coordinateτ ∈ [0, 1]. To enforce continuity, we have values on the left side of this narrow element P L = P (τ = 0) = P n+1,− and Q L = Q(τ = 0) = Q n+1,− and on the right P R = P (τ = 1) = P n+1,+ and Q R = Q(τ = 1) = Q n+1,+ . The goal is now to calculate the contribution of the term P · dQ/dt in the variational principle at a node for a discontinuous Galerkin finite element method, in which the variables are multivalued, as a limit ∆τ → 0 of the continuous Galerkin finite element method defined above. The result will, of course, depend on these continuous finite element expansions. Within these narrow elements, we choose a linear expansion in Q and a quadratic one for P across node t n+1 , but other choices can be explored as well, with α ∈ [0, 1] a free parameter. Hence, using (40) we obtain that The result is the jump in the coordinate Q times a weighted mean of the momentum P of the limit values across the time node at t n+1 , and it is independent of ∆τ . The other terms in the variational principle within these narrow elements contribute to zero in this limit ∆τ → 0 due to the absence of further time derivatives. Consequently, for an arbitrary discontinuous Galerkin finite element expansion within each element, the (intermediate) discrete variational principle is Similarly, the Kamiltonian principle becomes To derive definite time integration schemes, we still need to specify the expansions for P δ , Q δ , p δ and τ δ as well as the quadrature rules to evaluate the integral of the Hamiltonian.

Optimization criteria
The strategy is to optimise the properties of the time discretisation scheme (42) or (43) for various values of α, choices of the approximating polynomials and different quadratures used to approximate the time slab integrals. Formally, a polynomial of order n p will yield a time discretisation scheme that is O(n p +1). So far [7], the following results have been obtained for autonomous Hamiltonian systems: • Using piecewise constant polynomials and exact integration yields first-order accurate stable and symplectic schemes for α = 0, 1/2, 1. The variable in a time element is represented by its mean. Classic symplectic Euler schemes (see, e.g., Hairer et al. [9]) emerge for α = 0 and α = 1 in the autonomous case with H = H(Q, P ).
• For a second order scheme at least piecewise linear polynomials are required. Consequently, two coefficients or degrees of freedom represent each variable in a time element. Classic stable Störmer-Verlet schemes emerge. The four degrees of freedom per time slab in the scheme reduce to three because the variations yield that one of the variables becomes continuous. The scheme is a mixed (dis)continuous Galerkin scheme with in general two implicit stages and one explicit stage.
• Third-order schemes require at least quadratic polynomials. A symmetric quadrature (Simpson's rule) using three quadrature points for each variable leads to a six-stage scheme, only stable for α = 1/2. Dispersion analysis by Gagarina et al. [7] based on the discrete harmonic oscillator reveals that there are two parasitic modes, which influence is neutralised by initialising the scheme continuously, thus removing the initial jumps. The resulting scheme remains discontinuous with in general two coupled implicit and four explicit stages.
All these schemes will be extended to the non-autonomous case next via a Kamiltonian treatment of time in (8). This Kamiltonian can be used to monitor that the discrete Kamiltonian energy shows no drift and has bounded fluctuations.
We note that P n,+ = P n,− ≡ P n and p n,+ = p n,− ≡ p n have become continuous, while Q remains discontinuous. The final time discretisation can thus be viewed as a second-order mixed discontinuous and continuous Galerkin finite element approximation. When we also denote Q n,+ by Q n and replace τ n by t n and τ n+1/2 by t n+1/2 , then the final, compactly written non-autonomous Störmer-Verlet scheme reads The first two steps are, generally, implicit, while the last one is explicit, and the evaluation of time is at t n+1/2 , by construction.
Alternatively, we can reverse the approximation of Q and P (as well as replacing τ and p by q and τ ), such that while taking α = 0. The discrete variational principle then becomes A similar evaluation but now giving Q n+1,+ = Q n+1,− ≡ Q n+1 and q n+1,+ = q n+1,− ≡ q n+1 , then yields the following second Störmer-Verlet scheme for a non-autonomous Hamiltonian system after denoting the discontinuous variables as P n+1,+ = P n+1 , while one can derive that τ n+1,+ = τ n+1,− = τ n+1 = t n+1 is continuous when τ 0,+ = τ 0,− = 0 is continuous initially. Again the first two steps are implicit and the last one is explicit. Depending on the factual Hamiltonian system, one can choose one or the other Störmer-Verlet version. For the spatially discrete Miles' variational principle for potential flow water waves, we used the second version in Gagarina et al. [6]. It makes the free surface elevation continuous in time, while the velocity potential remains discontinuous in time.
In addition, it is easy to see that variations on the above Kamiltonian approaches include the following modified mid-point approximations of the Hamiltonian as well respectively, with corresponding changes in the discrete dynamics.

Variational Benney-Luke water wave model
All simulations of (33) were done within the finite element modelling environment Firedrake [11] on a quadrilateral mesh [12], to third order in space using quadratic Lagrange polynomials, and second or third order in time using symplectic integrators (49) or (59), respectively 1 .
To verify the numerical implementation, we compared numerical solutions with a standing wave solution of the linearised equations (i.e. (33) with = 0). We also compared small-amplitude nonlinear numerical solutions with an asymptotic KdV-soliton (i.e. (33) with both > 0 and µ > 0). A final validation case consists of comparisons between numerical solutions and wave phenomena observed during a soliton splash experiment 2 . Unless otherwise stated, regular rectangular elements are used with time step ∆t = 0.0028, quadratic Lagrange basis functions and Störmer-Verlet in time.

Exact standing wave solutions
An exact standing wave solution of the linearised counterpart of (31) (with = 0 and η R = 0) in a rectangular domain of size L x × L y is given by Φ(x, y, t) =B sin(ωt) cos where A is the wave amplitude, m 1 > 0, m 2 > 0 are integers and The above solution is used in a comparison with the numerical solution of the linearised Benney-Luke system with = 0. The visual comparison after two wave periods can be seen in Fig. 2. Energy fluctuations are small and bounded, as expected.

Small amplitude soliton solutions
Leading order asymptotic solutions in are derived from (31) (for η R = 0) by applying the following transformations in one spatial dimension Substitution of the above scalings into the variational principle (30b) on the real line with lim |ξ|→∞ F, η → 0, subsequent use of (65) to eliminate η and truncation to O( 2 ), Using that η = F ξ at leading order in from (65) into the equation resulting from (66c), gives the following form of the KdV-equation Employing a travelling wave solution ansatz η = η(χ = ξ − cτ ) with c > 0 into (67) then leads to the solutions (cf. Drazin & Johnson [5]) of the Korteweg-de-Vries (KdV) equation, which we used at t = 0 as initial conditions in the numerical computations. A simulation of a soliton reflecting against a solid wall using the nonlinear discrete Benney-Luke equations appears to be good, see lines correspond to the initial wave propagation phase, while dashed lines correspond to the after-reflection phase). The maximum wave amplitude during the reflection is twice the initial amplitude of the soliton, confirming the results of Power and Chwang [16]. The initial propagation before reflection matches the exact solution (68) well. In addition, the energy fluctuations are small and bounded 3 .

Soliton splash experiment
A soliton splash experiment with wave generation, propagation and reflection phenomena is used as a validation of the Benney-Luke model. The wave channel is sketched in Fig. 4. At one end of the channel, there is a lock with a sluice gate and at the other end, there is a linear symmetric V-shaped contraction. At time t = 0 the water is at rest with water level h 1 in the sluice compartment and h 0 < h 1 in the main channel. Several water levels h 0,1 have been considered but we focus on the event with smooth unbroken waves. Specific details about the soliton splash experiment, including the wavetank dimensions and location of sluice gate can be found in Table 2.
We model the sluice gate release in an approximate fashion by adding a horizontal gravitational component, effectively yielding an initial rest level η R (x, 0). This free surface rest level coincides with a field line of a modified gravitational potential, such that (in dimensional units) with 0 < x 1 < x 2 < L x , a sluice gate release time T s > 0, and initial water level difference h 1 −h 0 > 0 at t = 0. Initial conditions are therefore η(x, y, 0) = η R (x, 0) and φ(x, y, 0) = 0. Given the validation case, we took x 1 = s , x 2 = x 1 + 0.1 m (a sluice gate of thickness 10 cm) and T s ≈ h 1 /V g . The average still water level H 0 is calculated as In the numerical simulations, we took a wave amplitude of α 0 = 0.25 m and wavelength of 0 = 2.25 m, both based on observations from the soliton splash experiment.  The energy as function of time in the soliton splash simulations of the nonlinear and linear Benney-Luke system is displayed in Fig. 6. In both cases, the energy remains bounded and shows no drift. Additional figures from the soliton splash simulation can be found in Appendix B.

Conclusion
Mathematical modelling of water waves in wave tanks with a sluice gate has been illustrated by variational methods: by deriving a modified reduced model with a timedependent gravitational potential mimicking a removable "sluice gate" and by discretising this model in space and time using finite element methods, both within a variational framework. As validation, we simulated a soliton splash in a wave channel with a contraction. Future work will explore these methods for wave-energy devices and ships in modest to heavy seas.

Acknowledgments
This work was mostly undertaken during the program "Water Waves" at the Isaac Newton Institute of Mathematical Sciences (Cambridge, 2014). We thank the organisers Tom Bridges, Mark Groves, Paul Milewski and David Nicholls, as well as all participants for creating a stimulating environment. Special regards go to the Firedrake project team, in particular Colin Cotter, David Ham and Lawrence Mitchell, for assistance in using Firedrake. We acknowledge funding from EPSRC grant no. EP/L025388/1 with a link to the Dutch Technology Foundation STW for the project "FastFEM: behaviour of fast ships in waves".  sequence of weak formulations

B Soliton splash simulations
A comparison between simulations of the linear Benney-Luke system (with = 0) in Fig. 7 and the nonlinear Benney-Luke system (with = 0.55) in Fig. 5, shows that the (solitary) waves in the nonlinear simulations have more coherence, as expected, resulting in a good correspondence with the observed soliton splash event. The linear simulations produce a wave that arrives in the channel contraction later and exhibits a smaller amplitude amplification than the nonlinear one. Mid-channel and width averaged profiles of the nonlinear numerical solution for η can be found in Fig. 8. Double resolution runs with 800×40 elements confirm that our simulations have converged.