Variational modelling of extreme waves through oblique interaction of solitary waves: application to Mach reﬂection

. In this work, we model extreme waves that oc-cur due to Mach reﬂection through the intersection of two obliquely incident solitary waves. For a given range of incident angles and amplitudes, the Mach stem wave grows linearly in length and amplitude, reaching up to 4 times the amplitude of the incident waves. A variational approach is used to derive the bidirectional Benney–Luke equations, an asymptotic equivalent of the three-dimensional potential-ﬂow equations modelling water waves. This nonlinear and weakly dispersive model has the advantage of allowing wave propagation in two horizontal directions, which is not the case with the unidirectional Kadomtsev–Petviashvili (KP) equation used in most previous studies. A variational Galerkin ﬁnite-element method is applied to solve the system numerically in Firedrake with a second-order Störmer– Verlet temporal integration scheme, in order to obtain stable simulations that conserve the overall mass and energy of the system. Using this approach, we are able to get close to the 4-fold amplitude ampliﬁcation predicted by Miles.


Introduction
Offshore structures such as wind turbines, ships and platforms are designed to resist loads and stresses applied by winds, currents and water waves.These three factors can cause damage or destroy these structures when their effect is underestimated.Designers and engineers must take into account the effect of not only each of these phenomena separately, but also their interaction, which can increase their adverse effects.In this work, we focus on the impact of extreme waves created from the propagation of an obliquely incident solitary wave along the side of a ship (a wave-structure inter-action) or its impact with another identical obliquely incident wave (a wave-wave interaction).These two cases are mathematically equivalent since reflection at a rigid wall (represented here by the ship's side) is modelled through the boundary condition of no normal flow at the wall, which is equivalent to the intersection of two identical waves travelling in opposite directions, in which case a virtual wall is formed.The study of extreme, freak or rogue waves resulting from reflection at a wall or interaction of waves has spawned different theories in the last 50 years, some of which are subsequently reviewed.
The objective of the present work is to apply a theory first introduced by Miles (1977a, b) and based on experiments from Perroud (1957), where he described analytically the behaviour of an incident solitary wave interacting with a wall.For a specific range of angle of incidence ϕ i and scaled amplitude a i of the wave, the reflection of the soliton may result in three wave fronts: the incident and reflected waves (of respective amplitudes a i and a r ), as well as a Mach stem wave (of amplitude a w ) propagating along the wall with an increasing length (see Fig. 1).
This theory holds in the case of small-but-finite wave amplitude, shallow-but-finite water depth, and weak nonlinearity, that is, and is based on an interaction parameter, first defined as which enables the prediction of the amplitude and direction of propagation of each wave front.The most important observation is the transition at κ = 1 from a regular reflection Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
Figure 1.Top: top view of a channel containing an incident solitary wave propagating in the x direction with amplitude a i .The side wall is oblique and makes an angle ϕ i with the x direction.Bottom: top view of the reflection pattern when the incident wave impacts upon the wall.The pattern is composed of three waves: (1) the incident wave, (2) a reflected wave of amplitude a r that forms an angle ϕ r with the angle perpendicular to the wall, and (3) a Mach stem wave propagating along the wall with amplitude a w and an angle ϕ w with the wall.
(κ ≥ 1) to a Mach reflection (κ < 1), which has led to the following definition of the stem-wave amplification, , for κ ≥ 1, so that α w = a w /a i is the quotient of the stem-wave and incident-wave amplitudes.Equation (3) shows that at the transition point where κ = 1 the stem wave may grow up to 4 times the amplitude of the incident wave, leading to extreme loading on offshore structures.The aim of the present study is to develop a (numerical) model that can accurately simulate the evolution of the stem wave so that the distance and direction of propagation required to reach the 4-fold amplitude can be estimated.A challenging aspect is that it takes a long time and large distance of propagation before the stem wave reaches its maximum amplitude, which was a limiting factor in previous experimental and numerical studies.Kodama et al. (2009) extended Miles' theory to the Kadomtsev-Petviashvili (KP) limit, in which the assumptions are where H 0 , a 0 and λ 0 are the water depth, the wave amplitude and wavelength, respectively.While the KP limit still considers shallow-but-finite depth and small-but-finite amplitudes, the main difference with Miles' theory concerns the condition on the angle ϕ i .Yeh et al. (2010) explained that, in contrast to Miles' theory, wherein the soliton propagates in one direction only (the Korteweg-De Vries -KdV -limit), the KP limit assumes a quasi-two-dimensional approximation, and therefore the condition tan 2 ϕ i = O( ) cannot be simplified to ϕ 2 i = O( ) as in Miles' assumptions.The quasi-twodimensional KP soliton is not a solution of the KdV equation, but it can be transformed to an asymptotic KdV soliton via some manipulations detailed in Yeh et al. (2010).However, the width of the obtained KdV soliton is proportional to with a KP the scaled amplitude of the initial KP soliton, and it therefore depends on the angle ϕ i .This is physically unrealistic since the KdV soliton should have the same shape whatever its direction of propagation.For this reason, Yeh et al. ( 2010) brought a "high-order correction" to the solution, setting the amplitude of the KdV soliton to be so that its width depends on its amplitude a KdV , but not on any angle.Taking this into account, they slightly modified the definition (Eq.2) of the interaction parameter κ to where a i = a KdV /H 0 is the scaled amplitude of the incident wave, leading to what we will hereafter identify as the "modified Miles' theory" for the expected stem-wave amplification: Using this modified interaction parameter (7) in Eq. (3), they found much better agreement between previous numerical simulations (Funakoshi, 1980;Tanaka, 1993) and modified Miles' theory.Moreover, Kodama et al. (2009) showed that the stem wave resulting from the interaction of two solitary waves with small incident angles is an exact solution of the KP equation.Solving this KP equation, they could describe the exact solution depending on the angle of incidence and the amplitude of the initial waves, and validate their theory with numerical simulations (Kodama et al., 2009;Li et al., 2011).Both the amplitude and length of the stem wave indeed followed their predictions in the case of regular and Mach reflection.The numerical scheme could not simulate the highest amplitudes that Miles predicts for κ ≈ 1.Recently, Ablowitz and Curtis (2013) studied Mach reflection for the Benney-Luke approximation, showing that, in that case, modified Miles' theory applies asymptotically, leading to amplifications of up to 3.9.The purpose of the present work is to derive and apply a stable numerical scheme able to estimate the solution over a long distance of propagation, in order to model highamplitude waves and to confirm the transition from regular to Mach reflection happening for κ ≈ 1.We develop a model similar to the one of Benney and Luke (1964), which is an asymptotic approximation of the potential-flow equations for small-amplitude and long waves.Whilst it has the advantage of conserving both the nonlinear and dispersive properties of the waves (essential to the modelling of a freak wave, for instance), it does not require a mesh moving vertically with the free surface since the model is reduced to the horizontal plane.Pego and Quintero (1999) derived these modified Benney-Luke equations and Bokhove and Kalogirou (2016) recently used them to simulate a soliton splash resulting from a wave running in a restricted channel.Their simulations were in reasonably good agreement with experiments, which confirms that the Benney-Luke approximation is an accurate model of water waves.The present approaches are necessary to determine how, in future work, we can impose the line solitons on the wave makers to generate a 4-fold amplified wave in the middle of a wave basin and measure its impact on offshore structures.The variational technique used in the present approach enables us to express the equations as a Hamiltonian system to which robust time integrators can be applied (Hairer et al., 2006;Gagarina et al., 2016).The space and time Galerkin finite-element method used to discretize the present model ensures the overall conservation of mass, energy and momentum, which are essential in the highamplitude and long-distance propagating waves studied here.
The remainder of this paper is organized as follows: the modified Benney-Luke-type model is derived from the variational principle for an inviscid and incompressible fluid (Luke, 1967) in the potential-flow approximation, using the small-amplitude and small-dispersion scaling of Pego and Quintero (1999).In order to apply modified Miles' theory and verify our numerical results against Kodama's exact solution, the KP limit is obtained from the Benney-Luke approximation, leading to a new variational principle for KP.A careful scaling is then defined to obtain an asymptotic soliton solution of our present model, based on the exact solution of the KP equation from Kodama et al. (2009).The corresponding interaction parameter is consequently derived, leading to another version of modified Miles' theory (Eqs. 3 and 7), later used to compare our numerical simulations with respect to Miles' expectations.The finite-element method is then used to discretize the equations in space together with the second-order Störmer-Verlet temporal scheme that ensures stable simulations.Results are finally discussed and compared to the expectations. .Three-dimensional water-wave domain with rest depth H 0 , velocity potential φ(x, y, z, t), total depth h(x, y, t) and free surface deviation η(x, y, t).
2 Water-wave model

Introduction
Our water-wave model is derived using a variational approach that ensures conservation of mass, momentum and energy.In a basic sea state with extreme waves, these conservation properties are essential given the different length scales involved.Starting from Luke's variational principle for an inviscid fluid with a free surface (Luke, 1967), a model similar to the one derived by Benney and Luke (1964) for small-amplitude and long waves is obtained.The (numerical) method developed by Bokhove and Kalogirou (2016) is used to derive the relevant variational principle for our Benney-Luke model.This asymptotic model conserves the nonlinear and dispersive properties of the sea waves, which enables comparison with the Kadomtsev-Petviashvili (KP) model for which the modified Miles' theory as expressed in Eqs.(3) and (7) applies.

From Luke's variational principle to the Benney-Luke set of equations
Water-wave equations are often adequately described by the potential-flow approximation.In the absence of vorticity, the fluid velocity u = (u x , u y , u z ) can be expressed as the gradient of the so-called velocity potential φ(x, y, z), such that u = ∇φ.The deviation from the surface at rest H 0 is defined by η(x, y, t) so that the total depth h(x, y, t) can be expressed as h(x, y, t) = H 0 +η(x, y, t) (cf.Fig. 2).We consider a flat sea bed lying at z = 0, with vertical walls at ∂ b , where b is the horizontal plane of the bed coordinates b = {0 ≤ x ≤ L x , 0 ≤ y ≤ L y }.Luke (1967) described an inviscid and incompressible fluid with a free surface in the potential-flow approximation through the following variational principle: www.nonlin-processes-geophys.net/24/43/2017/ Nonlin.Processes Geophys., 24, 43-60, 2017 where g is the acceleration of gravity.The gradient ∇ is defined on b only, such that ∇ = (∂ x , ∂ y ) T is the horizontal gradient.The velocities at the walls and sea bed are assumed to be zero, that is, n • ∇φ = 0 on ∂ b , with n the outward horizontal normal and ∂ z φ = 0 at z = 0.The boundary conditions at the free surface z = h and the equations of motion in the domain are obtained from Eq. ( 8) as The amplitude parameter = a/H 0 1, with a the amplitude of the waves, and the small dispersion parameter µ = (H 0 /λ 0 ) 2 1, with λ 0 the horizontal wavelength, have been introduced by Milewski and Keller (1996) and Pego and Quintero (1999) to scale Eq. ( 8).The scaled variational principle is where This scaling focusses on small-amplitude long waves.From now on, the hats on the variables introduced in Eq. ( 11) are omitted.
To derive the Benney-Luke model, the velocity potential φ is expanded in terms of the sea-bed potential φ(x, y, 0, t) = (x, y, t) and the dispersion parameter µ as in Bokhove and Kalogirou (2016): Combining this expansion with the system of Eq. ( 9) and retaining terms up to second order, Eq. ( 12) becomes (see Bokhove and Kalogirou, 2016, for details) Substituting Eq. ( 13) into the variational principle Eq. ( 10) yields the variational principle under the Benney-Luke ap- Arbitrary variations in both and η, together with boundary conditions n • ∇ = 0 and n • ∇ = 0 at ∂ b , lead to the Benney-Luke equations Equation ( 15) will be solved numerically as explained in Sect. 4.However, to test our Benney-Luke model on modified Miles' theory (Eqs.3 and 7), it must first be compared to the KP theory for which Kodama et al. (2009) have shown that modified Miles' theory holds.

From the Benney-Luke set of equations to the Kadomtsev-Petviashvili equation
The Kadomtsev-Petviashvili equation for small-amplitude solitons can be derived from the Benney-Luke variational principle Eq. ( 14) and Eq. ( 15) through the transformations Substituting scalings Eq. ( 16) into Eq.( 15a), η can be expressed from as Substituting Eq. ( 16) into the transformed variational principle Eq. ( 14) yields Nonlin.Processes Geophys., 24, 43-60, 2017 www.nonlin-processes-geophys.net/24/43/2017/ Subsequent elimination of η using Eq. ( 17) and truncation to O( 2 ) gives the variational principle for KP in terms of η ≈ X : Note that we consider an infinite plane, with vanishing at the boundaries |X, Y | → ∞, such that the boundary terms arising from the integration by parts vanish in Eq. ( 19b).Since δ is arbitrary, the variational principle Eq. ( 19) yields the following equation for the leading-order scaled potential : From Eq. ( 17), at leading order in O( ), η can be expressed as η = X and, therefore, taking the partial derivative of Eq. ( 20) with respect to X leads to the KP equation for η: A solution of the KP Eq. ( 21) is found by substituting the following soliton solution ansatz, the form inspired by Yeh et al. (2010) Eq. ( 9), into Eq.( 21): where ϕ is the angle of incidence, A is the amplitude of the soliton, and B and C are coefficients to be determined via direct substitution.The KP soliton is then found to be with 3A/4 and A the prescribed amplitude.Using Eq. ( 17) at leading order, i.e. η = X , the solution for thus becomes 3 Comparison with modified Miles' theory and Kodama's exact solution 3.1 Introduction to Kodama's exact solution Kodama et al. (2009) have studied the reflection pattern for "symmetric V -shape initial waves consisting of two semi- infinite line solitons with the same amplitude", in a system of coordinates ( X, Y , τ ) related to our system of coordinates Eq. ( 16) (X, Y, τ ) via They solved the KP equation for which the surface deviation solution η is given by where A is the amplitude of the soliton, ϕ is the angle of incidence at the wall, and C is a constant defined as C ≡ 1 2 A+ 3 4 tan 2 ϕ.They showed that in this specific case, the transition from regular to Mach reflection occurs when tan ϕ = 2 A. (28) www.nonlin-processes-geophys.net/24/43/2017/ Nonlin.Processes Geophys., 24, 43-60, 2017 F. Gidel et al.: Variational model of obliquely interacting solitary waves Moreover, Kodama et al. (2009) defined exactly the incident, reflected and stem solitons resulting from the interaction as a O-type soliton in the case where tan ϕ > √ 2 A, and a (3142)type soliton in the case where tan ϕ < √ 2 A. The O-type soliton consists of two line solitons travelling in the x direction, each having a specific amplitude and angle with respect to the y axis (see Fig. 3).The (3142)-type soliton consists of two other line solitons, also travelling in the x direction with their own amplitudes and angles with respect to the y axis, but this soliton also has the property of being non-stationary, i.e. that while it propagates along the x axis, a new line soliton is progressively created and grows parallel to the y axis at the intersection of the two initial line solitons.In the case of both O-type and (3142)-type solitons, one of the line solitons can be associated with the incident solitary wave presented in the introduction, the second line solitons with the reflected wave (with a different amplitude and angle), and the intersection of the two line solitons as the stem wave, growing in length only when the angle of the incident wave is smaller than the critical angle (Eq.28).These two solitons are represented in Fig. 3, obtained from Kodama et al. (2009).A comparison between these theoretical solitons and those obtained numerically from the V-shape initial soliton showed very good agreement, confirming that the incident, reflected and stem waves described by Miles are indeed asymptotically equivalent to the O-type and (3412)-type solitons, depending on the initial angles.In the case of a symmetric initial pattern, that is, for two initial line solitons of equal amplitude and angle of incidence, Kodama et al. (2009) gave the expression of the maximal amplitude of the intersection wave as Since the condition tan ϕ = √ 2 A is equivalent to Miles' condition κ = 1, we can define the interaction parameter corresponding to the KP Eq. ( 26) as Substitution of the interaction parameter (Eq.30) into the amplification expectations (Eq.29) indeed yields Miles' predictions (Eq. 3) for α w = a max / A.

Application to the present Benney-Luke model
In Sect.2.3, the Benney-Luke model was reduced to the KP Eq. ( 21).This equation for the surface deviation η is slightly different from the one used by Kodama et al. (2009) and introduced in Eq. ( 26).In order to compare our numerical solutions to Kodama et al. (2009)'s result Eqs. ( 29)-( 30), our KP Eq. ( 21) is (re)scaled using the coefficients introduced in Eq. ( 25), which yields Eq. ( 26) used by Kodama et al. (2009).
Using the same transformations Eq. ( 25) in the KP soliton solution Eq. ( 27), we can obtain a solution for our KP Eq. ( 26) in terms of the original variables (X, Y, τ, η) introduced in Eq. ( 16), given by The connection between the above solution (Eq.31) and the previously presented solution (Eq.23) can be established by applying the following transformations in Eq. ( 31): with tan 2 ϕ, yielding the solution (Eq.23) derived in Sect.2.3.Therefore, applying scaling (Eq.32) to the critical condition (Eq.28) yields the critical condition for Eq. ( 21) and solution Eq. ( 23), given by We then apply scaling Eq. ( 16) to transform solution Eq. ( 23) for η back to the original Benney-Luke approximation Eq. ( 15) used in our simulations, in which case the asymptotic solutions for η and become where the soliton has been localized around the position (x 0 , y 0 ) at time t = t 0 .Finally, by setting the solutions (34) of the Benney-Luke equations can be rewritten as η(x, y, t) This solution is used as an initial condition at time t = 0 in the simulations.Condition Eq. ( 33) defines the following relation between ϕ i , a i and in our Benney-Luke scaling, for Eq. ( 15): This condition is equivalent to Miles' condition κ = 1 and therefore we can define our Benney-Luke interaction parameter as Note, however, that taking into account the remark from Kodama (2010) about the quasi two-dimensionality of the KP limit, as explained in the introduction, the interaction parameter defined in Eq. ( 38) must be corrected to in order to satisfy Miles' prediction Eq. ( 3).As shown in the potential-flow Eq. ( 9) for the Benney-Luke approximation, the small-amplitude parameter is defined as = a/H 0 .Therefore, in the specific case where a i = 1 and = a KdV /H 0 , the interaction parameter Eq. ( 7) is recovered.The diagram in Fig. 4 summarizes the equations and solutions derived thus far, in each scaling.In the next section, we explain how the Benney-Luke system of equations is discretized in both space and time in order to be solved numerically.

Numerical implementation
As a first step in the computational solution, the Benney-Luke model needs to be discretized in space and time, on a meshed domain.This section explains the methods used to discretize the domain and the equations.

Space discretization: finite-element method (FEM)
A continuous Galerkin finite-element method is used to discretize the solutions in space.The variables η and are approximated by the finite-element expansions where the subscript h denotes the discretized form of the solutions with basis functions ϕ j (x, y), and i, j ∈ [1, N ] with 2N unknowns.The Einstein notation for the implicit summation of repeated indices is used.To avoid the secondorder derivative in the fourth term of the variational principle Eq. ( 14), the auxiliary variable is introduced, so that, in the variational principle Eq. ( 14), the term µ 3 ( ) 2 can be written as which leads to the variational principle In keeping with Eq. ( 40), the second-order Galerkin expansion for q is now expressed as Substituting expansions Eqs. ( 40) and ( 44) into the variational principle Eq. ( 43) yields the space-discrete variational principle with ˙ i the time derivative of i .Its variation with temporal end-point conditions δ i (0 A matrix form of Eq. ( 46) can be found in Bokhove and Kalogirou (2016).Rather than using this matrix form directly, we accommodate only the spatial discretization using Firedrake (Rathgeber et al., 2016;Balay et al., 2016Balay et al., , 1997;;Dalcin et al., 2011;Hendrickson and Leland, 1995), "an automated system for the portable solution of partial differential equations using the finite element method (FEM)".This automated system uses the finite-element method to solve partial differential equations, and requires specification of the following: the domain in which the equations are solved, and the kind of mesh to use (e.g.quadrilaterals, the spatial dimension); the order and type of polynomials used; the type of expansion for the unknowns (e.g.continuous Galerkin, Lagrange polynomials); the function space of the unknowns and test functions; and the weak formulations discretized in time.
In the present case the domain is defined as a horizontal channel ending in an oblique wall, and quadrilaterals are used for its discretization (see details in Sect.5.1.2).Here, we chose to use quadratic polynomials to expand , q and η.The resulting weak formulations implemented in Firedrake in terms of h , q h and η h are the following: The forms given in Eq. ( 47) are convenient since they highlight the unknowns h , q h and η h as well as the test functions δ h , δq h and δη h .The final step is to discretize the equations in time, with a second-order Störmer-Verlet scheme, as explained in the next section.

Time discretization: second-order Störmer-Verlet scheme
After incorporating the FEM expansions, the space-discrete form of the variational principle Eq. ( 43) can be written in the Hamiltonian form where M ij and A ij are the mass and stiffness matrices, respectively defined as and the Hamiltonian is where Note that dH /dt = 0 due to skew symmetry.Gagarina et al. (2016) have shown that, for a generic Hamiltonian system in the form here with P = { i } and Q = {(M ij + µ/2A ij )η j }, robust variational time integrators conserving the overall mass and energy can be formulated.To derive these time schemes, P and Q are discretized on each time interval [t n , t n+1 ] as the approximated momentum P τ and coordinate Q τ and expanded with coefficients P m and Q m and linear continuous basis functions ψ m and ψ m : The linear basis functions ψ m and ψ m are continuous within each time interval, but admit discontinuities at the interface between two time slots.Therefore, to discretize Eq. ( 52), the notion of jumps [[.]] and averages {{.}} β α for a time-dependent function d(t) must be introduced (Gagarina et al., 2016): The coefficients α and β are real numbers defined such that α + β = 1 and α, β ≥ 0. The notation d n,± denotes the left and right traces of d(t) at time t n , that is, Discretization of the variational principle Eq. ( 52) then yields (Gagarina et al., 2016) where N is the number of finite time intervals [t n , t n+1 ] that divide the time domain [0, T ].Gagarina et al. (2016) showed that to obtain a second-order Störmer-Verlet scheme, P and Q must be discretized with mid-point and trapezoidal rules, respectively, that is, Substituting Eqs. ( 57)-( 58) into the discretized variational principle Eq. ( 56) yields The variations of Eq. ( 59), when augmented by endpoint conditions δP 0,− = δ(2P −1/2 − P −1,+ ) = 0, δQ 0,− = 0, δP N,+ = 0 and δQ N,+ = 0, yield the following scheme: Setting P n = αP n,+ + βP n,− with α = 0 and β = 1 ensures stability of the numerical scheme (Gagarina et al., 2016 Substituting these conditions into Eq.( 60) yields the continuity condition [[Q]] t n = 0 for Q in Eq. (60a), and the secondorder Störmer-Verlet scheme is recovered, with the stability condition where ω is the (maximum) frequency of the discrete waves.
Setting the vectors P = { i } and Q = M ij + µ 2 A ij η j , the variational principle Eq. ( 48) for Benney-Luke equations can therefore be discretized as in Eq. ( 61), leading to Eq. (A1) in Appendix Sect. A. Since the space discretization is performed internally within Firedrake, the weak formulations Eq. ( A1) can be implemented with the full form of the variables h , q h and η h and test functions δ h , δq h and δη h yielding the time discretization of Eq. ( 47), namely Time-step Eqs.(63a), ( 63b) and (63c) are implicit, while Eq. ( 63d) is explicit.Although the equations are nonlinear, the step Eqs. ( 63b), ( 63c) and ( 63d) are linear with respect to the unknowns, q n+1/2 h , η n+1 h and n+1 h , respectively.Therefore, linear solvers are used to solve the three weak formulations (Eq.63b, c, d), which reduces the computational cost by assembling the Jacobian matrix only once instead of computing it at each time step.The implementation of such linear and nonlinear solvers is straightforward in Firedrake, since functions that solve weak formulations for specific unknown and test functions already exist (Rathgeber et al., 2016;Balay et al., 1997Balay et al., , 2016;;Hendrickson and Leland, 1995;Dalcin et al., 2011).

Numerical results
In this section, the domain is specified and discretized in order to evaluate and η numerically.The numerical evolution of the stem-wave amplitude is compared to the predictions from our modified Miles' theory Eqs. ( 3) and ( 39).Finally, the angles of propagation of the reflected and stem waves are measured and compared to the values predicted by theory.

Orientation of the channel
The interaction of two solitary waves can be modelled using either two obliquely intersecting channels, with incident solitons propagating along each channel (see scheme a in Fig. 5), or from the reflection of a soliton at a wall with the no-normal flow condition at the wall (see scheme b in Fig. 5).While the first case a is more relevant to the theme of this paper, we choose to model case b to reduce the size of the domain by half and thus to reduce the computational cost.Since cases a and b are mathematically equivalent, the results and conclusions obtained with half of the domain will also be valid for the intersection of two oblique channels.
The domain is described by the length of the wall L w , the length of the channel L c , and the angle of incidence ϕ i .The channel should be long enough, compared to the wavelength of the incident wave, in order for the boundaries to be far enough from the initial soliton to be considered infinitely distant.From Eq. ( 34), the width of the initial soliton depends on √ 3 /4µ, and since µ is set to 0.02 for every simulation, the width of the soliton varies with from 2.5 (when = 0.20) to 4 (when = 0.12).We set L c = 5 to leave enough space between the extremities of the soliton and the boundary of the channel for every case, for the extremities to be considered infinitely distant from the soliton boundaries.To allow the stem wave to grow and reach its maximal amplitude, the wall also needs to be long compared to the wavelength.This constraint was a limit in previous numerical and experimental studies (Tanaka, 1993;Li et al., 2011) since it requires robust and stable numerical schemes and large wave basins.We set the wall length to 200 ≤ L w ≤ 600 depending on the value of , that is, more than 100 times the incident-wave width.When considering half of the domain represented in Fig. 5b, we chose to set the wall in the x direction, in which case the initial soliton must propagate in an oblique direction and is therefore equivalent to a KP soliton, as defined in  Eq. ( 36).Alternatively, we can let the initial soliton propagate in the x direction, in which case the wall is oblique and the expression of the KP-type soliton (Eq.36) can be simplified to a KdV-type soliton propagating in the x direction, as (Drazin and Johnson, 1989) η(x, y, t) The behaviour of the incident and stem waves in the cases of an oblique incident soliton (Eq.36) and a soliton propagating in the x direction only (Eq.64) are compared in Fig. 6.The initial solitons have amplitude a i = 1.0, smallamplitude parameter = 0.14 and small-dispersion parameter µ = 0.02.The angle between the direction of propagation of the solitons and the wall is ϕ i = π/6 in both cases.The dashed lines represent the evolution of the interpolated amplitude of incident solitons with time.While the initial amplitude was a i = 1.0 in both cases, we observe that both amplitudes first increase before decreasing to an asymptotic value slightly smaller than 1.0 (a i = 0.93).This behaviour is not expected for solitons, which should keep a permanent shape.However, we solve here the Benney-Luke equations for which the KP soliton is only an asymptotic (and hence not exact) solution because we recall that the transformation (Eq.16) from the Benney-Luke model to the KP theory is not exact since it requires a truncation to O( 2 ).In the numerical simulations represented in Fig. 6, = 0.14, so the condition O(1) is respected only asymptotically; this is a possible explanation of the observed variation in amplitude.Figure 6 shows that the incident KP and KdV-type solitons Eqs. ( 36) and ( 64) converge, and that both do so to the same surface deviation, a i = 0.93.This same limit shows that the approximation error from Benney-Luke to the KP soliton is asymptotically the same as from Benney-Luke to KdV.The stem-wave amplitudes (solid lines in Fig. 6) resulting from the interaction of the KP-type (Eq.36) and KdV-type (Eq.64) initial solitons with the wall are both amplified at the same speed and with the same amplification factor, which confirms that the KP-type and KdV-type initial solitons Eqs. ( 36) and ( 64) give the same results.The small variations in the curves are due to the mesh resolution which is not fine enough to resolve a regular amplitude.However, the computed approximation is sufficiently accurate to provide an estimate of the asymptotic amplitude of the stem wave.Since we have demonstrated that the two types of initial solitons Eqs. ( 36) and ( 64) evolve similarly to give the same results, subsequent simula-Figure 7. Domain discretization using quadrilaterals in Gmsh.In order to reduce computational requirements, mesh refinement is restricted to only the region adjacent to the wall.
tions will be conducted using only a unidirectional soliton, as defined by Eq. ( 64), which is a solution of both the KP and KdV equations.

Mesh
In order to evaluate and η at any position in the channel, the domain is discretized using quadrilaterals.This is done using the Gmsh mesh generator (Geuzaine and Remacle, 2009).Since the domain is large, we define a heterogeneous mesh with areas of higher refinement along the wall, where the solution needs to be more accurate.Moreover, the end of the domain is truncated with a blunt wall instead of the sharp angle, to avoid boundary quadrilaterals having internal angles that are too acute.The final domain comprising different mesh refinements is presented in Fig. 7, in which the insets show the aforementioned refined mesh and right-hand boundary quadrilateral elements.

Amplification of the stem wave
The numerical amplification of the stem wave is compared with the predictions of modified Miles' theory applied to our Benney-Luke model Eqs.( 3) and (39), namely The interaction parameter defined in Eq. (39) depends on three parameters: the scaled amplitude of the incident soliton a i , its angle of incidence ϕ i , and the small-amplitude parameter .From Miles' theory, a change in these parameters will modify the behaviour of the reflected and stem waves.
Figure 8 shows a comparison between predictions Eqs.(3) plitude.In performing the computations required for Fig. 8, we defined the maximal amplification as follows: when the stem wave reaches its maximal amplitude a w max , we measure the amplitude of the incident wave a i at the same x position.The new incident amplitude a i is used to adjust the interaction parameter and to compute the amplification of the stem wave α w = a w max /a i .The grid refinement is 0.25×0.25 in the finest area (e.g. at the wall) and 0.4 × 1.5 elsewhere.
The numerics follow the theoretical curve, but a slight difference between the present results and those expected from modified Miles is noticeable.As alluded to beforehand, we assume that this is due to the fact that the soliton used as an incident wave is an asymptotic but not an exact solution of the Benney-Luke equations.The scaling from Benney-Luke to KP is not exact but asymptotic, with a truncation at second order, which leads to a slight difference in the final wave amplification.This observation agrees with the conclusions of Ablowitz and Curtis (2013) on the asymptotic amplification of the stem wave in the case of the Benney-Luke model.The shift is probably also increased by the mesh resolution, which could be optimized to get a better estimate of the incident-wave amplitude in order to limit the error caused by its approximation.New simulations with higher mesh resolution are expected to verify the current results.However, the present Benney-Luke model still predicts the evolution of the stem-wave amplitude very well, enabling it to reach up to 3.6 times the initial amplitude.The stem-wave maximal amplification is reached for κ = 0.9733, marginally smaller than the κ = 1.0 predicted by Miles.While the model from Kodama et al. (2009) is expected to predict the evolution of the stem wave based on the KP equation perfectly, they were unable to reach more than 3.2 times the initial amplitude in their numerical simulations.

Angle of the stem and reflected waves
Miles' theory also predicts different directions of propagation of the stem and reflected waves in the cases of regular and Mach reflections.While in the first case, characterized by κ ≥ 1, the angle of the reflected wave ϕ r is expected to be equal to that of the incident soliton ϕ i , it should become larger than ϕ i in the case of Mach reflection, i.e. when κ < 1: Moreover, in the case of regular reflection, the stem wave is expected to propagate along the wall with a constant length, while for Mach reflection, its length should increase linearly, making a non-zero angle ϕ w with the wall: Predictions Eqs. ( 65) and (66) will be checked numerically next.

Regular reflection
Figure 9 shows numerical results and predictions for the case where κ = 1.12 ≥ 1.The wall makes an angle of 30 • with the direction of propagation of the initial solitary wave, hence ϕ i = 30 • .In the bottom-right plot of Fig. 9, there is an angle of 60 • between the reflected and stem waves, which means that the angle ϕ r between the reflected wave and the line perpendicular to the wall is equal to 30 • , that is, equal to ϕ i .This observation holds at any time, and therefore the expectations (Eq.65) for the reflected waves are satisfied in the case of regular reflection.The stem wave propagates along the wall without increasing in length, and therefore no angle can be measured between the stem wave and the wall, i.e. ϕ w = 0, as predicted in Eq. ( 66) for regular reflection.These results, together with Fig. 8 for the amplification of the stem wave, confirm modified Miles' theory in the case κ ≥ 1 for both the reflected and stem waves.

Mach reflection
Figure 10 shows numerical results and schematic expectations for the propagation of the reflected and stem wave for κ = 0.58 < 1.In the bottom-right sub-figure, the angle between the incident and reflected waves can be measured, as represented in the top-right sub-figure, in order to check that ϕ r is larger than ϕ i .The total angle ϕ r +ϕ i measures 70 • , with the initial incident angle set to ϕ i = 30 • .Therefore, ϕ r is 40   gle ϕ w with the wall.In the bottom-right figure, a top view of the numerical results at different times from t = 0.28 to t = 1.12 highlights the increase in the stem wave's length as it propagates along the wall.The dashed orange line that connects the solutions confirms that the wavelength increases linearly.

Conclusions and discussions
The present model Eq. ( 15), together with the new scaled interaction parameter Eq. ( 39), shows good agreement with the predictions of Miles regarding the amplification of the stem wave and the angles of the reflected and stem waves.Two different regimes can be observed in the numerical results, with different behaviours of the waves in the case of Mach and regular reflections, which confirms the conclusions obtained by Ablowitz and Curtis (2013) regarding the ability of the Benney-Luke model to predict reflection of obliquely incident solitary waves.Our simulations do not allow determination of the exact value of the interaction parameter at the transition from Mach to regular reflection, but currently the maximal amplification is reached at κ = 0.9733, which is very close to the predicted maximal amplification at κ = 1.0.The maximal amplification obtained herein is α w = 3.6, which is higher than the amplifications obtained with most previous models and experiments (Kodama et al., 2009;Li et al., 2011;Tanaka, 1993;Funakoshi, 1980), but still slightly lower than the expected 3.9 amplification from Ablowitz and Curtis (2013).This agrees with the conclusion of Ablowitz and Curtis (2013) concerning the impact of on the amplification near κ = 1.While they obtained the maximal amplification α w = 3.9 for = 0.10, our amplification α w = 3.6 is obtained for = 0.17, which is larger than 0.1 and thus leads to a larger difference with Miles' prediction of α w ≈ 4.Moreover, thanks to the robust scheme, which ensures stable simulations over the large domain despite the different length scales involved, used to derive and discretize our equations, our model is the first model able to describe numerically the dynamic development of the stem wave up to such high amplitudes.Previous studies (Kodama et al., 2009;Li et al., 2011;Tanaka, 1993;Funakoshi, 1980) were not able to attain such high amplifications because of numerical limitations such as insufficient computational resources.Ablowitz and Curtis (2013) obtained the highest numerical amplification α w = 3.9 by considering the final state, initialized asymptotically using the KP two-line solution.This last approach gives an accurate understanding of the asymptotic maximal amplification of the stem wave with the Benney-Luke model, but does not describe the development of the stem wave along the wall.The description and understanding of the wave propagation along the wall is however fundamental for application purposes.The present results, although currently limited by computational resources, allow us to consider the relevance of obliquely interacting solitary waves in maritime engineering.More advanced simulations should enable determination of the value of κ at the transition from Mach to regular reflection and to reach higher amplification of the stem wave.
There are some limits to the current model.As already concluded in previous studies, the wave needs to propagate over a long distance (relative to its wavelength) in order to reach its maximal amplitude.Consequently, the numerical domain needs to be large and the mesh fine enough to estimate the wave crests accurately.This numerical requirement increases the computational time.A compromise between the accuracy of the simulations and the running time is therefore needed.This constraint is important because near the transition from Mach to regular reflection a slight change in the incident wave amplitude modifies dramatically the interaction parameter and consequently the predictions of the stem and reflected waves.Therefore, a careful analysis of the numerical results must be made.For the same reason, simulations for κ ≈ 1 and large amplifications α w ≈ 4 are extremely difficult to obtain, since a slight change in the initial settings (a i , and ϕ i ) modifies completely the behaviour of the resulting waves.Indeed, Li et al. (2011) conjectured that the transition between Mach and regular reflection in the neighbourhood of κ = 1 might appear gradually and not as abruptly as expected from Miles' predictions Eq. (66).
One may wonder how likely it is that solitary waves would undergo reflection in an open ocean.Interaction of obliquely incident waves on the sides of a ship leads to an increasing wave amplitude, sometimes reaching the deck.This phenomenon is called "green water" and has been studied experimentally and numerically by the Maritime Research Institute Netherlands (MARIN) to limit the damage caused by waves on ships (Buchner et al., 2014).When the incident wave interacts with a ship moving downwind, the effective ship length increases, leaving more time for the stem wave to develop to its maximum amplitude.Peterson et al. (2003) also studied the formation of extreme waves in shallow water and explained under which conditions they are likely to occur and threaten ships.Kalogirou and Bokhove (2016) developed numerical models of waves impacting buoys and ships.An extension of our oblique-wave interaction simulations to www.nonlin-processes-geophys.net/24/43/2017/ Nonlin.Processes Geophys., 24, 43-60, 2017 F. Gidel et al.: Variational model of obliquely interacting solitary waves wave interactions with ships will be an interesting extension of our present work.The present model can also be used to predict the impact of extreme (i.e.freak or rogue) waves on structures.Indeed, when the stem wave reaches more than twice the amplitude of the incident wave, it can be viewed as a freak wave since it has similar properties in terms of nonlinearity, dispersivity and high amplitude.Table 1 shows the distance required by the stem wave to reach more than twice the incident-wave amplitude in several cases parameterized by different values of .For each value of the small-amplitude parameter , the numerical (dimensionless) distance L n required to reach at least twice the amplitude of the initial wave has been measured from the simulations.Then, the definition of the smallamplitude parameter = a 0 /H 0 and the choice of a sea state with characteristic wave height a 0 = 3 m enables computation of the corresponding water depth H 0 .The physical distance L r required by the wave to propagate up to twice the characteristic wave height can then be obtained from scaling Eq. ( 11), using L r = L n × H 0 / √ µ.The value of the smalldispersion parameter µ is set to 0.02 as in the results section.Finally, the wavelength λ 0 can be obtained from the definition of the small-dispersion parameter µ = (H 0 /λ 0 ) 2 .In a wave tank where waves can be generated in different directions, the angle of propagation and initial profile of two solitary waves can be defined from the asymptotically exact solution Eq. (36) of our model Eq. ( 15), so that their interaction will lead to a stem wave.The evolution of the stem wave can be predicted from the present model, so an offshore structure such as a scaled ship or a wind turbine can be placed at a position where the stem wave will reach more than twice the initial amplitude of the solitary waves.For instance, a scaling of 1/10 between the values of H 0 , L r and λ 0 in Table 1 and experiments leads to achievable incidentwave amplitudes and distances of propagation in MARIN's shallow-water basin.From the amplitude of the stem wave at a given position, the impact of the wave on structures can be estimated and the predictions yielded by the model tests can be validated.The model can help the maritime industry to design safer offshore structures that can resist extreme-wave impacts.
Finally, the present work can also be used as a starting point for the modelling of the interaction of three obliquely incident line solitons, which should lead to a 9-fold-amplified resulting wave that can also be generated in wave tanks.1

Data availability
The implementation of our discretization of the Benney-Luke equations is an example in Firedrake, www.firedrake.org (Bokhove and Kalogirou, 2016).In addition, the expanded program used to perform the numerical simulations is archived using Zenodo: https://zenodo.org/badge/latestdoi/79556994.

Figure 3 .
Figure 3. O-type and (3142)-type solitons as represented by Kodama et al. (2009).Top: evolution (from left to right) of the O-type soliton, consisting of two line solitons with different amplitudes and angles with respect to the y axis.As it propagates, the shape of this soliton remains unchanged.Bottom: evolution (from left to right) of the (3142)-type soliton, consisting of two line solitons travelling in the x direction with different angles and amplitudes.As the soliton propagates, a new line soliton is created at the intersection of the two initial line solitons, leading to a stem wave.Figure obtained from Kodama et al. (2009).© IOP Publishing.Reproduced with permission.All rights reserved.

Figure 4 .
Figure 4. Schematic plan showing the link between the scaling of the three systems of equations involved in the derivation of the exact solution and critical condition for which Miles' and Kodama's predictions hold in the Benney-Luke approximation.

Figure 5 .
Figure5.Definition of the domains in the two cases described in the text: (a) intersection of two channels, with two obliquely incident solitons interacting at a virtual wall, and (b) half of the domain with a soliton propagating in one channel and colliding with an oblique wall.This wall is in the x direction (in which case the soliton has a two-dimensional propagation of direction) or oblique, in which case the incident soliton propagates in the one-dimensional direction (x).

Figure 6 .
Figure6.Soliton surface deviations obtained for an initial amplitude a i = 1.0 and angle ϕ i = π/6 rad.Blue: behaviour of the incident (dashed line) and stem (full line) waves when the incident soliton propagates in an oblique direction; red: behaviour of the incident (dashed line) and stem (full line) waves when the incident soliton propagates in one direction.The dashes lines essentially coincide after t > 30.

Figure 9 .
Figure9.Numerical results and predictions for the reflected and stem waves in the case of regular reflection, i.e. κ > 1. Top left: top view of the numerical evolution of the incident, reflected and stem waves.Top right: schematic plan view of the expected evolution of the stem and reflected waves at two different times t 1 and t 2 with t 1 < t 2 .The stem wave should propagate along the wall with constant length.The angle ϕ r of the reflected wave is expected to be constant and equal to the incident-wave angle ϕ i .Bottom centre: side view of the time evolution of the incident, reflected and stem waves, highlighting the amplification of the stem-wave amplitude compared to the initial solitary-wave height.

Figure 10 .
Figure10.Numerical results and predictions for the reflected and stem waves in the case of Mach reflection, i.e. κ < 1. Top left: schematic plan view of the numerical evolution of the incident, reflected and stem waves.Top right: top-view scheme of the predicted evolution of the stem and reflected waves at two different times t 1 and t 2 with t 1 < t 2 .The stem wave should grow linearly in length, leading to an angle ϕ w > 0 with the wall.The angle ϕ r of the reflected wave is expected to be constant and larger than the incident-wave angle ϕ i .Bottom centre: side view of the time evolution of the incident, reflected and stem waves, highlighting the amplification of the stem-wave amplitude compared to the initial solitary-wave height.

Table 1 .
Prediction of the minimal distance needed by the stem wave to reach at least twice its initial amplitude in a sea state with characteristic wave height a 0 = 3 m.The dispersion parameter µ is set to 0.02, while the small-amplitude parameter varies from 0.12 to 0.20, leading to different wave evolutions.The numerical distance needed to reach more than twice the incident-wave amplitude is measured from the numerical simulations.The corresponding water depth, real distance of propagation and wavelength are computed from the definition of , µ and scaling (11).These values are approximate.