Numerical Experiments on Extreme Waves Through Oblique–Soliton Interactions

Extreme water-wave motion is investigated analytically and numerically by considering two-soliton and three-soliton interactions on a horizontal plane. We successfully determine numerically that soliton solutions of the unidirectional Kadomtsev–Petviashvili equation (KPE), with equal far-field individual amplitudes, survive reasonably well in the bidirectional and higher-order Benney–Luke equations (BLE). A well-known exact two-soliton solution of the KPE on the infinite horizontal plane is used to seed the BLE at an initial time, and we confirm that the KPE-fourfold amplification approximately persists. More interestingly, a known three-soliton solution of the KPE is analysed further to assess its eight- or ninefold amplification, the latter of which exists in only a special and difficult-to-attain limit. This solution leads to an extreme splash at one point in space and time. Subsequently, we seed the BLE with this three-soliton solution at a suitable initial time to establish the maximum amplification: it is approximately 7.8 for a KPE amplification of 8.4. Herein, the computational domain and solutions are truncated approximately to a fully periodic or half-periodic channel geometry of sufficient size, essentially leading to cnoidal-wave solutions. Moreover, special geometric (finite-element) variational integrators in space and time have been used in order to eradicate artificial numerical damping of, in particular, wave amplitude.

by the aforementioned rogue-wave investigations of short-crested waves of extremely high amplitude. We will focus on a particular type of extreme wave emerging in crossing seas, with waves described analytically as solutions of the Kadomtsev-Petviashvilli equation (KPE) [15] and resolved numerically in higher-order wave models, for which the KPE is an asymptotic approximation.
Water waves are commonly modelled by the potential-flow wave equations (e.g., [24,31]), using incompressible flow with a free surface in the absence of vorticity so that the three-dimensonal velocity field can be expressed as the gradient of a velocity potential, or by asymptotic approximations thereof. Both direct representations, as partial differential equations, or spectral representations are used of these potentialflow equations [22]. Solutions to two approximations to these potential-flow equations will be considered here: the bidirectional Benney-Luke equations (BLE) [4,7,28], and the unidirectional KPE [15,18]. The KPE is also an asymptotic approximation of the BLE. The advantage of using this simplified KPE is that it has exact solutions in the form of web solitons [1,18]. The KPE is a two-dimensional extension of the onedimensional Korteweg-De-Vries (KdV) equation, based on introducing or keeping weak interactions in the other horizontal, i.e. lateral, direction [15]. The famous KdVsech soliton is therefore also a single-line soliton solution of the KPE. Web solutions of the KPE comprise interacting line solitons, each consisting of such a sech-solution with with an orientation differing from that of the others in the far field on an infinite plane. As per the KdV-equation, the KPE also allows periodic wave solutions in the form of cnoidal waves (but these exist for the KPE in two horizontal dimensions), leading to high-amplitude short-crested wave solutions where these "individual" cnoidal waves interact with or cross one another [13]. Such short-crested high-amplitude crossing waves have also been observed, both at sea and in the laboratory [13,21], which links them to some extent to the rogue waves in crossing seas discussed above.
Alternatively, certain two-soliton solutions of the KPE on the infinite plane can be used to describe approximately the interaction of one soliton travelling along a wall encountering, and thus interacting with, a corner or sharp turn in that wall, see Fig. 1a, b). This corresponds to a certain KPE-solution of two-line solitons, with the solid wall seen as the line of symmetry where the normal velocity equals zero. For a certain angle of this corner or turn, the amplification of an incoming soliton of amplitudeÃ can be shown analytically to become fourfold over time [18,25], i.e. to reach an asymptotic amplitude of 4Ã when it has travelled for a while along the straight wall after the turn; for the two interacting line solitons this fourfold amplification is sustained and propagating on an infinite plane. One difference between the infinite-plane and finitedomain solutions is that in the latter case it takes time for the solution interacting with the wall to reach asymptotically the fourfold amplification [18,25]. Note that the KPE is unidirectional so the introduction of actual walls is not allowed, since there cannot be any factual reflection, and the above conclusions are based on comparison with laboratory experiments [23] as well as on numerical simulations of the BLE [11].
As part of a fluid-dynamics' demonstration conducted in 2010 [8], it was shown that a soliton-complex travelling along a channel led to a tenfold amplification into a so-called soliton-splash at the apex of a linear contraction found at the channel end. A simplified set-up instead considers a sech-soliton travelling along the channel which, when interacting with the two turns at the contraction entrance, can lead to two [1,2] [1, 2] [3,4] [ 3,4] [1, 2] [3,4] [ 3,4] [3, 4] t = 0 t = T (a) (b) [1,2] [1, 2] [3,4] [ 5,6] (c) [3,4] [ 5,6] (d) [3,4] [ 3,4] [ 3,4] [ 3,4] [1, 2] [5,6] t = 0 t = T 1 [1,2] [ 5,6] [3, 4] t = T 2 Fig. 1 Sketches of (a, b) two or (c, d) three interacting line solitons. (a, c) portray the infinite-plane cases, (b) shows a soliton turning a corner, and (d) shows a soliton running into a linear contraction. The various far-field line solitons involved, e.g. [1,2], [3,4], [5,6] have been indicated, including their phase shift due to interactions waves growing towards this maximum fourfold amplification in the KPE, and these two colliding amplified waves can in turn form the soliton splash at the apex, see also Fig. 1c, d. Of course, in the closed contraction, reflection occurs at and around the apex, so that the amplification could be larger than in the infinite-plane KPE-analogue of a similar set-up. In the latter case on an infinite plane, three-line solitons will be interacting, with the two interactions of two-times-two-line solitons colliding at one point in space and time into a localised soliton splash, which subsequently decomposes again into two-times-two-line solitons. Several researchers have analysed such three-soliton interactions on an infinite horizontal plane within the KPE framework. Below follows a brief overview. Biondini et al. [5] analysed the interactions of line solitons that lead to large amplitudes. For three interacting solitons, six wavenumbers are needed, which can be ordered as k 1 < k 2 < k 3 < 0 < k 4 < k 5 < k 6 as we will see later. Biondini et al. [5] state an upper bound of the maximum splash amplitude of 1 2 (k max − k min ) 2 with far-field solitons of amplitude 1 2 , which we will analyse in detail for a three-soliton interaction (note that here k max = k 6 and k min = k 1 ). The three-line solitons have amplitudes 1 2 (k 2 − k 1 ) 2 , 1 2 (k 6 − k 5 ) 2 , 1 2 (k 4 − k 3 ) 2 in the far field. However, whether the upper-bound maximum is obtained depends on the details of the soliton configuration. Kodama [19,20] calculated the three-line-soliton interactions using the values (k 1 , k 2 , k 3 , k 4 , k 5 , k 6 ) = (−1.501, −0.501, −0.5, 0.5, 0.501, 1.501), leading to a maximum amplification of circa nine times, in fact 8.6; this is in line with Biondini's estimate given above, leading here to a just-above-ninefold amplification (for a video of such a splash in a miniature-tank set-up see [32]). Using these coefficients, Baker [2] and Gidel [10] checked and reproduced these calculations. Overarching questions are whether such high-amplitude solutions of the KPE will survive as approximate solutions of the potential-flow water wave equations, or higher-order bidirectional approximations thereof, and, furthermore, to what extent such amplifications can be observed in reality?
The fourfold amplification of a soliton against a channel bend as counterpart of two interacting line solitons has been observed in experiments. Li et al. [23] observed a nearly threefold amplification in Mach-stem reflection, i.e. the two-line-soliton analogue of a soliton travelling along a wall and encountering a corner. In Gidel et al. [11], numerical simulations were performed to check whether the fourfold amplification predicted in the KPE is valid for the interaction of a soliton encountering a corner in a wall, while using the more accurate, higher-order (in wave amplitude and dispersion parameters) bidirectional BLE. Such simulations are more challenging than imposing the solution for two-line solitons as an initial-value problem, because the amplitude of the soliton after it turns into a Mach-stem wave at the corner grows only slowly, over time, to its asymptotic value. Such simulations of the Benney-Luke system showed that the maximum amplification obtained reached 3.6, somewhat short of the fourfold amplification in the KPE or the amplification of 3.9 observed for simulations in which the KPE solution was imposed as a pair of initial values in the BLE [1].
The goals of the current paper are "extreme-wave" extensions of these investigations -of two interacting line solitons with equal far-field amplitudes-into the realm of three interacting line solitons of equal far-field amplitudes. The strategy to achieve these goals is as follows: • to derive the conditions for which three-line solitons of equal amplitudeÃ in the far field reach a maximum-amplitude factor of nine, i.e. amplitude 9Ã at one point in space and time, as solution of the KPE; and, • to establish numerically the extent to which this three-line-soliton solution of the asymptotic KPE remains valid in the bidirectional and more accurate higher-order BLE.
The outline of the paper is as follows. The BLE and KPE are introduced in Sect. 2 as well as their respective coordinate systems and the relation between them. The threeline-soliton solution of the KPE is introduced in Sect. 3 in order to allow initialisation within the BLE. Moreover, a proof of the maximum amplification of these solitons is given and detailed further in an Appendix. Numerical simulations of the amplification of interacting line solitons in the BLE, initialised with exact solutions of the KPE, are shown and interpreted in Sect. 4 before we draw a conclusion in Sect. 5.

Mathematical Models
In this section, two mathematical approaches for modelling water waves are considered-the Benney-Luke equations and the Kadomtsev-Petviashvilli equationincluding the scalings used to arrive at a standardised form.

The Benney-Luke Equations
We define a Cartesian coordinate system (x, y, z), in which (x, y) are the horizontal coordinates and z is the vertical coordinate, the latter of which is positive above a flat sea bed at z = 0. In dimensionless form, the water depth is described by h(x, y, t) = 1 + η(x, y, t), where t is time, η(x, y, t) is the free-surface deviation from the rest depth, and is an amplitude parameter. The domain is denoted by Ω, with z ∈ [0, h(x, y, t)] and horizontal extent Ω h . In deriving the BLE, potential-flow theory, which considers the fluid as irrotational, is applied. This allows us to write the fluid velocity in terms of a scalar velocity potential φ through u(x, y, is the three-dimensional Cartesian gradient operator. The corresponding gradient operator in two-dimensional horizontal space will be denoted by ∇ = (∂ x , ∂ y ).
The Benney-Luke approximation is based on the assumption of long and shallow water waves and allows us to reduce the dynamics to two-dimensions by tracking only the free-surface evolution and the velocity potential at a specified height. The BLE were first derived in [4] and in a variational manner also in [7,11,28]. Only these variational versions of the BLE will be considered here. The non-dimensional governing BLE are given by where Φ(x, y, t) = φ(x, y, z = 0, t) is the sea-bed potential. The two boundary conditions (1c), (1d) represent no-flux conditions at the boundary of the domain ∂Ω h , where n is the outward unit normal vector. The Benney-Luke equations include two parameters, and μ, which are both assumed to be small in the current approximation; in particular, the amplitude parameter is defined as the ratio of wave amplitude over the mean water depth, and μ is a dispersion parameter proportional to ratio of wavelength over water depth.
The energy of the Benney-Luke system in (1) is defined by where the first term in the integrand is the potential energy and the remaining terms therein represent the kinetic energy of the system. The total energy defined above is conserved in time, with a proof given in Appendix A.

The Kadomtsev-Petviashvili Equation
The Kadomtsev-Petviashvili equation is a two-dimensional analogue of the KdV equation and describes unidirectional propagation of waves in two dimensions. While the KdV equation is appropriate for wave propagation in the x-direction, the KPE also includes weak diffraction effects in the y-direction (as can be seen below in (3) where the y-scale is smaller than the x-scale). The KPE can be obtained from the BLE by using the following transformation (see also the combined asymptotic and variational approach in [7] and taking μ = 2 , where the domain of validity of the new variables is (X ,Ỹ ,τ ) ∈ R 2 × R + . Subsequently, we introduce the formal perturbation expansions whereΨX =ũ. Substituting these expansions into Eqs. (1a), (1b) and retaining terms up to order 2 , yields the leading-order equation, cf. [11] ∂X 2∂τũ + 3ũ∂Xũ + 1 3 ∂XXXũ + ∂ỸỸũ = 0.
We finally take another transformation in order to cast Eq. (5) in the "standard" KPE form [18,21], as follows applying this on Eq. (5) results in the well-known KPE [15] of the form The coordinates and variables in the KPE (7) are connected to those in the BLE (1) through the rescaling below, obtained by combining scalings (3) and (6), with μ = 2 , and given by which will be used later to set up the KPE-solutions in X , Y , τ -space for the BLE. The function K (X , Y , τ ) emerges in the next section as u(X , Y , τ ) ≡ 2∂ X X ln K (X , Y , τ ) .

Maximum Amplification with Three Equal-Amplitude Line Solitons
After transforming the KPE in the standard form (7), line-soliton solutions of the following form can be constructed using Hirota's transformation [14] u where function K (X , Y , τ ) can be obtained from the Wronskian in which f (n) i denotes the nth partial derivative of f i with respect to X , and the f i are such that they satisfy the linear equations ∂ Y f i = ∂ X X f i and ∂ τ f i = −∂ X X X f i [5,18]. Particular line-soliton solutions are obtained by taking where θ j = k j X + k 2 j Y − k 3 j τ , with coefficients k j being ordered as k 1 < k 2 < · · · < k M , and the a i j 's are the elements of an N × M matrix A (see also [18]).
The solution written in the form ( , the latter explaining why we had used six k j 's, for j = 1, 2, . . . , 6, in the introduction. These three cases are now discussed in turn, so as to set up the analysis and expression for the three-line-soliton case with equal amplitudes in the far field.

Single-Line Soliton
For the single-line soliton, we use N = 1, M = 2 in Eqs. (10) and take a 11 = 1, a 12 = a > 0, resulting in K = f 1 = e θ 1 + a e θ 2 . The solution then becomes which is a line soliton with centreline at X + (k 1 + k 2 )Y = (k 2 1 + k 1 k 2 + k 2 2 )τ + (ln a)/(k 1 − k 2 ), as found by setting the argument of the sech 2 function to zero. Here, the parameter a essentially determines the location of the soliton solution. Note that when k 1 + k 2 = 0 and k 1 < k 2 , this soliton propagates in the positive X -direction and its crest lies parallel to the Y -axis, an observation used later for the case with three-line solitons.

Two Interacting Line or O-Solitons of Equal Amplitude
The two-line soliton has the following set-up with the 2 × 4 matrix such that where a, b are positive constants. Hence, the interim function K becomes and the parameters k 1,2,3,4 satisfy the ordering k 1 < k 2 < 0 < k 3 < k 4 without loss of generality. On the one hand, in the limit Y → ∞ when X > 0, the second term in (14b) dominates since k 4 > k 3 , which recovers the [1,2] single-line soliton u [1,2] On the other hand, taking the limit Y → ∞ when X < 0 causes the second term in (14c) to dominate since k 2 > k 1 , and we obtain the [3,4] single-line soliton u [3,4] The above arguments also hold when in the case where one of the [3,4] or [1,2] single-line soliton propagates in the positive X -direction. We refer to these cases as limits in the three-line soliton case.

Three Interacting Line Solitons of Equal Amplitude
It is possible to find an analytical exact solution of the KPE that describes the interaction of three solitary waves, as shown in [2,18]. Such a solution is called a (3, 3)-soliton and it follows from equation (9) with N = 3, M = 6, and using the matrix where the positive constants a, b, c shift the location of the solution; in what follows we take a = b = c = 1 which corresponds to no shift. The functions f i , i = 1, 2, 3 in (10b) are hence given by and are such that they would each form a single-line soliton when considered alone, as described in Sect. 3.1. The interim function K in (10a) then becomes with the following conditions satisfied between parameters k j , j = 1, . . . , 6: and (x 2 , y 2 ) are the southwest, southeast, northwest, and northeast corners, respectively −k 4 ; more details will be provided later in Sect. 3.4). Note that, without loss of generality, herein the middle line soliton is aligned to propagate in the positive Xdirection, leading to the choice k 3 + k 4 = 0. We also note that the above ordering of parameters k j has been directed by the nature of the far-field expressions of the line solitons, which is explained next. By using the KPE solution (9) in the case of three-line solitons, we can find approximate solutions for each solitary wave far away from the region of interaction. This can be achieved by considering waves along lines X + (k i + k j )Y = c 0 , for some constant c 0 and for |Y | → ∞, where the pair (i, j) takes the values (1, 2), (3,4), or (5,6), and keeping only the dominant exponential terms in (22). The three solitary wave solutions in separation, i.e. far from the interaction region, are then approximately given by the following expressions (see also Sect. 3.4) with The solitary waves indicated by subscripts [5,6], [3,4] or [1,2] are travelling in a northeast, east or southeast direction in the (x, y)-plane when viewed from above, as seen in Fig. 2. The equations in (23) provide information about the amplitude, speed and angle of each solitary wave (angle θ is between wave and positive y-axis, see Fig. 2). In particular, the wave characterised by subscript [i, j] has amplitude [3,4] or [5,6].
Without loss of generality, we take the middle [3,4] soliton to be parallel to the Y -axis and travelling in the positive X -direction, and the other two-line solitons, [1,2] and [5,6], are assumed to have equal and opposite angles ±θ with the Y -axis as that should make the amplitude at the collision of the three solitons maximal. Finally, the outer two solitons are assumed to have equal amplitude, taken to be 1/λ times the [3,4] soliton. We will ultimately take λ = 1 to determine the maximum amplitude of three colliding solitons of equal amplitudeÃ, but here we allow exploration around λ ≈ 1.
Consequently, the following conditions hold for the parameters k 1 , . . . , k 6 , with angle θ > 0. The equations in (24) can be solved to give the solutions In Eqs. (25), we have introduced a new parameter δ defined by which allows us to eliminate tan θ from the expressions for k i . Note that by using definition (26a), we can express the amplitude in terms of θ and δ as follows noting that we do not yet know when and where the maximum will occur. The symmetry of the set-up immediately seems to make clear that the maximum will occur when the three intersection points coincide at a space-time point (X * , Y * , τ * ) based on the far-field expressions.

Proof of Maximum Amplification
We begin by finding the location (X * , Y * , τ * ) in space-time where the KPE solution u given in (9), with K given in (22), attains its maximum. The argument is entirely geometrical: we calculate the two intersection points of the centrelines of the [5,6] and [3,4] line solitons, as well as the the centrelines of the [1,2] and [3,4] line solitons. This yields two Y -positions, Y [1,2] and Y [5,6] , the mean of which is Y * . The X -position X * and time τ * when the maximum is attained then follow by satisfying Y [1,2] = Y [5,6] = Y * . The solution (22) is rewritten as follows We note that all expressions in (27) are equivalent, but each one of them will be used to obtain the dominant term in the limits Y → ±∞ and for different values of X . For X > 0, Y → ∞, we find that e θ 4 e θ 3 and e θ 6 e θ 5 in (27b), making the term underlined therein dominant, thereby leading to the [1, 2] soliton (23a). However, for X < 0, Y → −∞ we obtain the double-underlined term as the dominant one (in particular, substitute X + (k 1 + k 2 )Y = 0 to determine the dominant term), which is the one we need for comparing the pair of intersection points of the (perturbed) [5,6] and [1,2] line solitons with the (perturbed) [3,4] line soliton. For X ≈ 0, Y → ∞, we find that e θ 6 e θ 5 and e θ 1 e θ 2 , so the underlined term in (27c) dominates over the other terms and leads to the [3,4] soliton (23c). For X < 0, Y → ∞, we find that e θ 1 e θ 2 and e θ 3 e θ 4 in (27d), making the underlined term therein dominant, thereby leading to the [5,6] soliton (23b).
Recall that we have set a = b = c = 1. Therefore in the far-field, we obtain, K [5,6] ≈ A 135 e θ 1 +θ 3 e θ 5 + A 136 A 135 e θ 6 = A 135 e θ 1 +θ 3 e θ 5 + G [5,6] with G [1,2] = A 235 /A 135 , G [5,6] = A 136 /A 135 . From (22) we find that in which we have used k 1 = −k 6 , k 2 = −k 5 , k 3 = −k 4 and k 4 = k 5 , the latter of which is true when δ = 0 (see (25)). The limit δ → 0 provides the optimal angle: it will be used later to establish the maximum amplification. We note that, in the above calculations, it is important to take either the branches for X < 0 or those for X > 0 in determining the average Y = Y * , on which plane the maximum will have to occur. The corresponding three-line-soliton solutions for X < 0 and Y → ±∞ thus become, using (28), u [5,6] u [3,4] The centrelines of these line solitons coincide with the arguments of the sech 2 functions in (30) being zero, i.e. when while again using k 1 = −k 6 , k 2 = −k 5 , k 3 = −k 4 . Hence, the centreline of the [3,4] soliton resides at Substitution of this relation (32) respectively, the average of which yields The Y -coordinates in relations (33) equal Y * when τ * = 0 for which X * = 0 from (32). Hence, the space-time point with the maximum amplitude is where we have assumed that this is where the three-soliton centrelines coalesce. The next step is to assess the maximum possible amplification After defining θ lmn = (k 2 l + k 2 m + k 2 n )Y * , with l, m, n = 1, 2, . . . , 6 and l = m = n, we observe and determine the following relations Using these relations together with those in (29), the following relations emerge when taking the derivatives in ∂ X X K as well as the definition (coming from (34)) we find that where equal pairs have been indicated (the four pairs are in principle different prior to fixing the k i 's, except for using the properties k 1 = −k 6 , k 2 = −k 5 , k 3 = −k 4 ). Hence, the amplification at the peak relative to the largest soliton with amplitudeÃ and with λ ≥ 1 becomes as follows in which the definitions of k 1 , . . . , k 6 found in (25), with their dependence on δ, λ,Ã, have been substituted as well as (29), the latter also in the limit δ → 0. Quantitative illustrations of the dependence upon δ of both the absolute error in (40b) (computed relative to (40a)), and the maximum peak-amplification factor predicted by (40b) are shown in Fig. 3 for the specific value λ = 1.
Since both λ > 0 and δ > 0, the series form (40b) implies that the peakamplification factor, μ p say, is maximised when δ = 0, at which its value is given by (40c). This follows from the positive-definiteness of the coefficient of √ δ in (40), wherein neglected terms are of order O(δ 5/2 ). Immediately obvious from (40b) is that ∂ ∂δ μ p is of order O(δ −1/2 ) as δ → 0, whereas ∂ ∂λ μ p is of order O(δ 1/2 ) as δ → 0. Thus the maximum value of μ p is respectively cusped and smooth with respect to variations in δ and λ.
That the maximum peak-amplification factor (40c) as δ → 0 is the unique global maximum (with respect to both δ and λ) can be proven explicitly using the algebraic manipulator Maple. Differentiating the exact expression (40a) with respect to δ yields a rational fraction in which both numerator and denominator are cumbersome transcendental functions of λ, δ > 0; their presentation is facilitated using the pictorial form in Fig. 4, in which the negative sign initiating the numerator and positive-definiteness of all other terms confirms that ∂ ∂δ μ < 0 for all λ, δ > 0. The conclusion ∂ ∂λ μ < 0 for all λ, δ > 0 can be similarly reached; the details are omitted. Hence, since (40a) is positive-definite and partial derivatives with respect to both its independent variables are globally negative-definite in the domain of definition, the maximum value (40c) as δ → 0 is unique. Fig. 3 (a) Log-log plot of the absolute error, as a function of δ, between exact maximum peak-amplification factor (40a) and series approximation (40b) for the "base state" λ = 1. The dashed line indicates 23 δ 5/2 /(2 11/4 3 3/2 ), confirming the error in (40b) with λ = 1, departure from which at the smallest value of δ signifies entry of the calculations into the rounding plateau of the double-precision arithmetic used. (b) Semi-log plot of the maximum peak-amplification factor (40a) as a function of δ, for λ = 1: the nine-fold amplification as δ → 0 is evident Fig. 4 Maple-generated explicit forms of the numerator dPn and denominator dPn of the derivative with respect to δ of (40a). Clearly the numerator and denominator are respectively negative-and positive-definite for all λ, δ > 0. An analogous inference follows when the derivative with respect to λ is considered, though the corresponding formulae are not presented Note that Y * in (34) diverges when k 4 → k 5 in the limit k 5 − k 4 = δ Ã → 0, as follows and also that the term A 146 F → 0 when δ → 0. Moreover, in the above calculation we have used that A 245 → 0 and A 236 → 0 when δ → 0, as shown in (29). However, it may be useful to adjust the A-matrix in (20) by using a, b, c = 1 such that Y * resides by construction at Y = 0, if possible. We have therefore established that the maximum amplification occurs for δ = 0, shown both asymptotically from the Taylor expansion in (40) and graphically, see Fig. 5. The reason that the above result comprises a valid rather than heuristic proof, despite the fact that it is based on the far-field expressions, is as follows. Prior to reaching the maximum at (X * , Y * , τ * ), there is a phase shift between the branches of [1,2], [5,6] for Y < Y * and Y > Y * , seemingly meaning that the far-field expressions cannot be used to find the interaction point geometrically. However, the phase shift swaps directions when it passes through the time at which the maximum occurs. That is, the symmetry of the set-up implies that the phase shift becomes zero at the maximum point (X * , Y * , τ * ) in space-time. Given that we have calculated and therefore know (X * , Y * , τ * ), an a posteriori proof that it is indeed the maximum is provided in Appendix B.

Use of Geometric/Conservative Finite-Element Method
The BLE are a Hamiltonian pair of wave equations that conserve mass, energy and phase-space volume. Since we are interested in assessing numerically the maximum wave amplification of interacting solitons, it is important to choose a numerical method that preserves wave amplitude and, preferably, the entire geometric or Hamiltonian structure of the underlying partial differential equations. Without such conservation In all simulations, we take = 0.05, μ = 2 . The expression taken for tan θ = 2 9 1/6 1 4 √ is such that to make the tangent of the angle between η 0 [1,2] and the y-axis for BLE equal to 1/4. SP2/SP3 concern simulations in a domain periodic in the x-direction Table 2 Values of parameters k i , obtained from (25) with the parameters in Table 1 Simulation properties, the maximum amplification obtained numerically could be too low due to artificial dissipation; additionally, (computationally expensive) high resolution may be required to minimise such dissipation. We therefore solve the BLE by employing a continuous Galerkin (CG) finite-element method (FEM) based on a geometric or variational space-time discretisation. In particular, we employ a second-order Störmer-Verlet scheme in time, and Lagrange interpolation polynomials of order one, two, three and four to discretise space, denoted by CG1/CG2/CG3/CG4 in the rest of this section. The reader is directed to references [7,8,11] for details on this geometric-discretisation approach. In addition, the method is implemented and available in Firedrake [29], which is "an automated system for the solution of partial differential equations using the finite element method".

Numerical Strategies
Two solutions of the KPE will be used as initial conditions for simulations of the BLE, in order to test whether the calculated maximum amplification is reached and/or maintained. We therefore run two sets of simulations, each using the parameter values shown in Tables 1 and 2, and described as follows: • simulations SP2 use initial conditions formed by the two-line soliton solution (9), with K (X , Y , τ ) given in (14) and k 1 < k 2 = 0 = k 3 < k 4 ; and, • simulations SP3 use initial conditions formed by the three-line soliton solution (9), with K (X , Y , τ ) given in (22) and the k i 's given in (25).
However, these KPE solutions hold on an infinite horizontal plane, and the region of interaction, where the far-field line solitons cross, propagates in the x-direction. Consequently, a numerical simulation ideally needs to occur within a domain that is sufficiently large to admit influence or decay, to nearly zero, of the effect of the inter-actions near the domain boundaries. The remaining issue is then how to deal with the moving single-line solitons at these domain edges. To achieve energy conservation, such far-field single-line solitons should be normal to the domain walls: clearly, this would necessitate awkwardly shaped domains, which concomitantly places challenging computational demands on the numerical methodology deployed. Accordingly, the following strategy has been devised for the simulations. The solutions can be set to become approximately periodic when the domain is sufficiently large. Consider Fig. 2 as one cell in a domain either periodic in X -or periodic in both the X -and Y -directions, in which we have placed the X -shaped [1,2], [5,6] line solitons exactly at the corners of a periodic domain. Such domains can be patched together to form crossing seas. Hence, a domain periodic in only the X -direction or in both X -and Y -directions is chosen. This approach has the additional advantage that energy will be approximately conserved, with any energy oscillations remaining bounded and decaying to zero as O(Δt 2 ) when the time step Δt → 0, given the second-order geometric time integrators and geometric spatial FEM used. However, while the reconstructed η is periodic, the velocity potential Φ is not and requires to be split into Φ = (U 0 x + c 0 ) +Φ for a suitable choice of U 0 and c 0 . The rewritten BLE for the periodic η,Φ are given below in Sect. 4.3.
Furthermore, in the numerical simulations, we halve the X -and Y -directional periodic domain with respect to Y = Y * . After that, whilst retaining X -directional periodicity, we impose Neumann boundary conditions on the top and the bottom Y = Y * of the half domain using the symmetry of K . These Neumann boundary conditions are approximate -with ∂ y Φ 0 = 0, ∂ y η 0 = 0 -and required for the proof, because ∂ y Φ 0 and ∂ y η 0 are not exactly zero on the top boundary whereon there will be "kinks" due to the periodicity procedure imposed. Note that K is symmetric about K is symmetric with respect to Y * + p or Y * − p. Thus, the symmetry allows ∂ Y K = 0 on Y = Y * + pk for an integer k. Therefore, we can halve the X -, Y -directional periodic domain about Y = Y * , or Y = Y * + p while imposing Neumann boundary condition on Y = Y * and Y = Y * + p. To assess the approximations made, we will shortly compare two simulations, one of which uses the full periodic domain and the other the half-periodic domain.
The resolution for each simulation is provided in Table 3. Meshes are regular and, given Δx, Δy, the time step is determined by trial and error such that the energy oscillations are stable, which we will confirm per simulation when discussing numerical results. By using a posteriori the linear dispersion relationship of the BLE and the stability criterion for the Störmer-Verlet time-stepping scheme, as well as our choice of Δx, Δy, we verified that our chosen time step was (linearly) stable.

Initial Conditions and Set-Up in (Half) Periodic Domain
As described in the previous section, the solutions can be made approximately periodic when the domain is sufficiently large. Assuming a doubly periodic domain with and are both satisfied on Ω h = {(x, y) ∈ [x 1 , x 2 ] × [y * , y 2 ]} for appropriately chosen in which we use transformations (8) with t = t 0 . For the given initial conditions, U 0 can be calculated using the periodicity ofΦ, and c 0 is found by settingΦ(x 1 , y) = Φ(x 2 , y) = 0, as given below in which c 0 (y) is defined to make U 0 x 1 + c 0 = Φ 0 (x 1 , y). Hence, when defining the line y) such that c 0 = F(0).

Simulation SP2 with Two-Soliton Interaction
The initial condition for the case with two interacting solitons in an x-periodic domain is set up as follows. Given the KPE solution u(X , Y , τ ) and K (X , Y , τ ) in (14) for the two interacting web solitons, we use (43) to find the initial condition for the BLE at τ = 0 and corresponding t = 0. Furthermore, the four values k 1 , k 2 = k 3 = 0, k 4 = −k 1 used are defined in Table 2. We introduce the initial condition using the expressions of the far-field solitons [3,4] , for a [3,4] from the arguments in (19). We choose the ends of the domains x 1 , x 2 , y 2 such that to make the [1,2] line soliton (45a) exactly pass through the corner (x 2 , y 2 ), and the [3,4] line soliton (45b) pass through the corner (x 1 , y 2 ). Firstly, we choose y 2 to be sufficiently larger than y * , here, y 2 = y * + 40 = −20 + 40 = 20. We note that when a full periodic domain is used, the lower end of the domain is at y 1 = y * − 40. We can subsequently find x 2 such that η 0[1,2] (x 2 , y 2 ) = max(η 0 [1,2] ), that is where the argument of the sech 2 -function (45c) is zero, a [1,2] (x 2 , y 2 ) = 0. Following a similar argument, we can obtain x 1 by setting a [3,4] (x 1 , y 2 ) = 0 using (45d). Consequently we find These initial conditions for SP2 are shown in Fig. 6.   [1,2] = 1/4 = − tan θ [3,4] , so that both line solitons exactly pass through the respective corners of the domain

Simulation SP3 with Three-Soliton Interaction
The initial condition for the case with three interacting web-solitons in an x-periodic domain is set up as follows. Given the KPE solution u(X , Y , τ ) and K (X , Y , τ ) in (27) for the three interacting web-solitons, we define the initial condition for the BLE at t = t 0 = −200 (the corresponding value of τ is found from (8a)). Note that we start the simulation at a negative time t 0 < 0 to allow the solution to reach its maximum at t = 0, according to the results of Sect. 3.4. Furthermore, the six values k 1 , k 2 , k 3 , k 4 = −k 3 , k 5 = −k 2 , k 6 = −k 1 used in the simulations are defined in Table 2.

Discussion of Numerical Simulations
In discussing the simulations, we will address the following points: first, we demonstrate that the difference in simulations between the fully periodic and x-periodic domain is sufficiently minor to warrant usage of the smaller x-periodic domain; second, simulations with the higher-order polynomial resolution CG2 and CG3 are more suitable for the above numerical experiments than the simulation with CG1; third, BLE can explain short-crested high-amplitude crossing waves; and, finally, BLE numerical solutions are consistent with ones of KPE, that is, their difference is within order . An   η and the relative difference is shown in Fig. 8. The more computationally expensive simulation crashes before the end-time is reached, but the proximity between the two simulations shows that we are warranted to use the half-domain computational strategy for SP3.

SP2 Simulations
The evolution of the wave amplitude in time, the maximum amplification and the ratio of the two is presented in Fig. 9, for the three polynomial orders CG1/CG2/CG3. The profiles obtained with CG1 are overshot compared to those of CG2 and CG3, and compared with the exact solutions of the KPE. Hence, the simulation with CG1 seems unreliable with insufficient spatial and possibly temporal resolution, even though the energy oscillations appear to be bounded. Clearly, the results obtained using CG2 and CG3 are close to each other with their point-wise difference in L ∞ at most 10 −4 . We can therefore conclude that CG2 and CG3 have converged, and that these computations are therefore reliable. Because computations using CG3 are computationally much more expensive, we focus on explaining the CG2 results next. Starting from the initial conditions seen in Fig. 6, the solution evolves in time as seen in Fig. 11; top views of solutions at two different times t = 25, 50 are shown for (rescaled) KPE and BLE, as well as cross-sections of the BLE solution in y at the same times. The location and value of the solution maximum is also seen in the top and bottom panels. The evolution of the respective energy (2) is shown in Fig. 10, indicating bounded oscillations and hence a stable numerical integration. Furthermore, Fig. 9a shows that the maximum of η is maintained around 1.6, while Fig. 9b illustrates thatÃ is around 0.41. So their amplification turns out to about 3.8 as seen in Fig. 9c. We note that the oscillations are due to dispersion effects in BLE. Unlike the exact KPE solution η K P , the numerically-obtained η initially grows in amplitude. Subsequently, dispersion separates small ripples behind the main solitary waves, so both the maximum andÃ are seen to decrease for a while. By the time the maximum andÃ start to increase again, small ripples have circulated around the periodic domain and meet the main solitary waves. That is why the amplification is oscillating somewhat, but still remains around 3.9. Finally, we note that our result is consistent with the amplification of circa 3.6 found in Gidel et al. [11] and that of 3.9 found in Ablowitz and Curtis [1].

SP3 Simulations
In SP3 simulations with three soliton interactions, an amplification of 7.8 is attained despite the initial increasing, separation of ripples and oscillation of amplitudes, as seen in Fig. 12, similarly to the phenomena observed in SP2 where a 5% error was found between numerically-obtained and exact amplifications (3.8 vs. 4.0). The amplification of 7.8 for three-soliton interactions is within 7.25% error when compared to the amplification of 8.41 for the exact solution of the KPE, obtained for the same value of δ ≈ 0.0014 (see Fig. 5 and Table 2). The BLE numerical solutions are observed to be only different relative to the exact solution of the KPE within O( ), with = 0.05 here, see Figs. 14 and 15. This means that the BLE modelling is consistent with the KPE modelling, in that extreme-wave propagation and creation exists and is sustained in the BLE as well. The bounded energy fluctuations in Fig. 13 confirm that SP3-CG2 is stable. Finally, the time evolution of the free surface can be seen in the video available in [16], where all variables are transformed to dimensional units using  the time evolution and spatial location of the free-surface maximum, for both KPE and BLE. We note that running the simulation for much longer allows for repetition of the maximum amplification in time, due to the cnoidal-wave structure of the solution (see also video reference in [16]).

Summary and Conclusions
Extreme waves may arise randomly in crossing seas comprising waves aligned with two or more directions of travel. Linear superposition with weak nonlinearity has been proposed to explain so-called "everyday" extreme waves, those with an amplitude of at least twice that of those in the surrounding ambient sea. Extreme waves may alternatively arise through higher-order nonlinear effects in statistical distributions; these have been simulated and observed to lead to extreme waves. In the former case, it has been proposed to reserve the term "rogue waves" for exceptionally high and steep water waves. Thus motivated, we have investigated exact and numerical "rogue-wave" solu- tions of water-wave equations for crossing seas, but here in a deterministic manner. Two exact web-solitons have been analysed for the unidirectional Kadomtsev-Petviashvili equation (KPE), and numerical solutions have been simulated for the bi-directional and higher-order Benney-Luke equations (BLE) in two horizontal dimensions, the latter seeded at an initial time by either one of the two exact web-soliton solutions of the KPE. The first exact solution of the KPE is well known and consists of two main soliton branches of amplitudeÃ, interacting at an angle, and leading to a soliton branch with fourfold amplitude 4Ã. The second exact solution of the KPE is less well-known and consists of three main soliton branches, each of amplitudeÃ in the far field, involving waves coming from three directions, leading to a wave splash of extreme height at one point in space and time. We have analysed this exact three-line soliton solution in detail under a symmetric set-up and show in a novel analysis that its maximum amplification peaks at 9Ã for a certain angle between a main solitary wave travelling in a given direction (we use the positive Cartesian x-direction) and two other solitary waves whose directions of travel are symmetrically disposed thereto. However, due to a phase shift, it is currently possible to prove the theoretical ninefold amplification only analytically for the KPE, in terms of a suitably defined small parameter δ, and it is an open question how to set up the numerical BLE simulation in the limit δ → 0. This is because the spatial location of the maximum amplification is determined only for positive values of δ, as shown in Sect. 3.4. Finally, we note that the limit of δ → 0 also provides the optimal angle for maximum amplification which is given by tan θ = 8Ã for KPE (this is found from Eq. (26b) with λ = 1). The corresponding optimal angle for BLE (and the respective dimensional equations), θ BLE , can be calculated through tan θ BLE = √ ( 9 2 ) 1/6 tan θ = 12 A H 0 (by using (8) and (26)), with dimensional far-field soliton amplitude A and (rest) water depth H 0 .
Computational simulation of such solutions is challenging and potentially cumbersome given both the speed of wave propagation and infinite nature of the 2-D solution domain. Given the symmetries inherent in the two KPE web-soliton solutions, we have (necessarily artificially) imposed them as initial conditions on a 2-D computational domain that is sufficiently large (in an experimentally determined sense) and periodic in both x-and y-directions. Additionally, by using a symmetry relative to a fixed value of y, we have shown that it suffices to conduct computations on only a half-domain, bounded by two solid walls in the y-direction and periodic in the x-direction. By doing so, we have effectively created cnoidal-wave solutions of crossing seas (see, for example, Fig. 16). Hence, we have seeded simulations of the BLE with two exact solutions of the KPE at some initial time, and have used geometric or variational (finite-element) integrators to discretise the BLE in space and time. The numerical methodology is specifically designed to avert artificial numerical damping of wave amplitudes, as it allows conservation of phase-space volume and mass, and keeps energy oscillations small and bounded. Simulations at different resolutions (using polynomial or p-refinement) reveal that these two types of extreme waves or web-solitons reach maximum-amplification factors of ∼ 3.6-4.0 when the corresponding factor for the exact solutions of the KPE is 4.0, and of ∼ 7.8 when it is 8.4. Deviations from the exact KPE solutions have moreover emerged because of minor secondary waves, created after the initial time, which induce small-scale oscillations in the far-field soliton(s) of original amplitudeÃ; this has the effect of eroding the (computed) maximum amplification. Computed maximum amplitudes and amplifications are summarised in Table 4 and depicted in Figs. 9 and 12.
The final soliton splash examined in our simulations does indeed create a rogue wave in the sense of the definition thereof based on the abnormality-index calculation (see Appendix C). In the case of three interacting line solitons, we find an abnormality index of AI = max(η)/H s = 4.0, where H s is the "significant wave height", 1 and the maximum value of η is the yellow text annotating Fig. 15a. Details of this calculation are included in Appendix C.
Several research extensions based on the present study emerge naturally. An amalgamation of bespoke numerics and appropriate wave theory could seed simulations able to achieve (more closely) the theoretically predicted nine-fold maximum amplification factor. The three-line-soliton solution for the BLE could be simulated in a suitable channel geometry, as indicated in Fig. 1. The two web-soliton solutions could be simulated using potential-flow water-wave equations, possibly combined with a wave-breaking parametrisation scheme, in both periodic and channel-geometry settings. Such potential-flow simulations could offer a means of assessing how realistic these "KPE/BLE-rogue" waves are and whether they would endure sufficiently to attain their theoretical maximum amplitudes.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Conservation of Energy
The total energy E(t) as defined in Eq. (2) is conserved at all times, as shown below. Multiplying Eq. (1b) by (−∂ t Φ) and integrating in Ω h , gives where the other equations in (1) were used to get to the final result.
We define where the cosh arguments are assumed to satisfy b 1 ≤ b 2 ≤ b 3 . The derivatives of this function are so we find Since (55c) is the reciprocal of a summation of the cosh-functions, (55c) has the maximum at X = 0. Moreover, (55d) also has its maximum at X = 0 due to the following inequality Therefore, the maximum of ∂ X X ln G lies at X = 0. Finally, we claim that G(X ) = K (X , Y * , 0); this is clear by finding K (X , Y * , 0) from (50a) and comparing with the expression for G in (54), where we set such that b 1 < b 2 < b 3 is satisfied. This leads us to u(X , Y * , 0) = 2∂ X X ln K (X , Y * , 0) = 2∂ X X ln G(X ). As a consequence, the maximum of u(X , Y * , 0) lies at X = 0. Now that the two steps above are proved, we can conclude that max u(X , Y , τ ) = max

C Abnormality Index
Based on the results of simulation SP3, we can calculate the abnormality index AI defined as the ratio of the maximum wave amplification, max(η), to the significant wave height, H s . The maximum of η is attained at t = 0 as shown in Fig. 15a, and the significant wave height is defined as four times the standard deviation of the surface elevation [9]. Hence we need to determine H s at t = 0, i.e. at the time when the surface elevation is maximised.