## Abstract

A long wave multi-dimensional approximation of shallow-water waves is the bi-directional Benney–Luke (BL) equation. It yields the well-known Kadomtsev–Petviashvili (KP) equation in a quasi one-directional limit. A direct perturbation method is developed; it uses underlying conservation laws to determine the slow evolution of parameters of two space-dimensional, non-decaying solutions to the BL equation. These non-decaying solutions are perturbations of recently studied web solutions of the KP equation. New numerical simulations, based on windowing methods which are effective for non-decaying data, are presented. These simulations support the analytical results and elucidate the relationship between the KP and the BL equations and are also used to obtain amplitude information regarding particular web solutions. Additional dissipative perturbations to the BL equation are also studied.

## 1. Introduction

The modelling of the propagation of waves on the surface of water is a classical problem and remains an active and important area of research. However, the mathematical formulation of the problem, also known as the Euler water wave equations, presents a number of difficulties such as strong nonlinearity and an unknown location of the boundary that makes finding solutions challenging. To deal with this problem, researchers have employed a number of approximate models to the full system of equations describing surface flow.

In weakly nonlinear water waves whose depth is shallow with respect to the wavelength of the wave, a maximally balanced multi-dimensional approximation to the Euler water wave equations is the Kadomtsev–Petviashvilli (KP) equation. There are two signs associated with the KP equation depending on surface tension (ST); these are termed the KP-I (large ST) and KP-II equations (small ST). Here, we assume that ST is sufficiently small such that the KP-II equation is the relevant approximation to the Euler water wave equations (cf. [1]). Furthermore, in order to derive the KP equation, one also assumes that the height of the wave is small compared with the depth, and that the flow is quasi one-dimensional; i.e. transverse variations are slower than those in the primary direction of propagation (cf. [2] for a derivation and details). While other two-dimensional, shallow-water approximations to the Euler equations can be found, the KP equation is the unique model that is maximally balanced. This, and the fact that it has a Lax Pair and many known closed form solutions makes the KP equation an important model.

In recent years, a novel class of solutions, called ‘web solutions’, of the KP equation have been discovered and analysed. Building on work in [3] and the resonant solutions found in [4] and [5], in [6] the web solutions were obtained in terms of Wronskians via Sato's formalism [7]. Important properties and pertinent theorems concerning web solutions can be found in [8]. Likewise, experimental and observational evidence exists (cf. [9–11]) that indicates that these web solutions are persistent in nature.

From a mathematical perspective, it is important to understand whether these solutions are robust to perturbations. In this regard, it is useful to analyse improved approximations to the Euler equations which themselves yield the KP equation. An approximation that is the object of study in this paper is the Benney–Luke (BL) equation (cf. [12]). The BL equation is also a shallow-water approximation to the Euler equation which arises naturally as an approximation to water waves. However, unlike the KP equation, the BL equation allows for two-directional waves. Also important in our study is an asymptotically equivalent modification of the BL equation, which unlike the KP equation, has a dispersion relationship in which the group velocity of arbitrarily large wavenumbers is bounded; we consider this a BL-type equation as well. The latter BL equation, aside from being a two-directional approximation to the Euler equations, is expected to be well posed over a broader class of initial conditions than the KP equation.

Previous studies of the BL equation have appeared in the literature. The analysis in [4] begins with the BL equation, and a variant of the BL equation was studied numerically in the context of the Mach-reflection problem in [13]. In each of these papers, particular web solutions, the so-called ‘Y-junction’ solutions, were studied. While this paper addresses the case of small ST in the BL equation, we note that the large ST case has also been studied. Rigorous results establishing the existence of lump soliton solutions to the BL equation appear in [14], and the dynamics of these lump solitons are investigated in [15]. Rigorous results on the connection between the KP and BL equations have appeared in [16].

Recently, we developed a perturbation theory of web-type solutions of the BL equation based on techniques from the theory of integrable systems. We established in [17] that the web solutions of the KP equation persist for asymptotically long periods of time in the BL equation. Thus, the web solutions are robust to the higher order effects caused by the BL equation.

The purpose of this paper is to: (i) develop a direct perturbation method based on conservation laws of the BL equation and to (ii) present new computations of the BL equation which are consistent with web-type initial data that is non-decaying. In one dimension, the conservation-law approach is widely used (cf. [18]). This is because of its inherent simplicity and does not require one to employ more complex integrable system methods. However, in two dimensions since conserved quantities are expressed in terms of integrals, and the web solutions do not decay in one of the spatial directions the integrals are infinite. To deal with this difficulty, we present a modified definition of the conserved quantities of KP that works for the web solutions.

With this modification, the conservation law approach now provides a direct method to perform perturbation analysis (via multiple scales methods) on the web solutions. As indicated above it does not require integrable systems machinery and readily achieves the same results. Using the convenience of this conservation-law approach, we are also able to analyse the impact of a typical type of dissipation (cf. [19]) on web solutions in the BL equation. We chose a linear local dissipative term which creates ‘shelves’ behind the web solution. In one dimension, due to mean term interactions, small amplitude shelves of long extent are known to arise from Korteweg-de Vries (KdV) models with such a dissipative term (cf. [18,19]).

While the type of dissipation we study is particular, the method we present can, in principle, be used for many other, more accurate, dissipation models (cf. [20]). In this paper, our aim is not to decide on the best dissipative model. Instead, we explain how the direct conservation law approach needs to be modified when dissipation is present. We further exhibit the growth of shelves from the linear local dissipative model. This method can also be extended to higher order approximations of water waves. However, carrying out this effort is outside the scope of this paper.

The asymptotic methods in this paper rely on separating the web solutions into what we define as near- and far-field components. Using our conservation law-based asymptotic method, the far-field components are studied analytically. The near-field component is not required in the leading order perturbation analysis. However, the near-field is where the various parts of a web-solution interact. Therefore, we employ new numerical simulations of the BL equation to investigate the near-field, or interaction region. In this paper, we in particular look at both Y-junctions and another class of web solutions, the X-waves (or O-type solutions in the terminology of Chakravarty & Kodama [8]).

We also examine the Y-junction solution further by studying the amplification of the so-called ‘Mach stem’ (cf. [5]). By varying the small parameter representing the ‘shallowness’ of the water and the angle of opening of the ‘Y’, we show how the amplification ratio (§3*b*) of the Mach stem in the BL equation varies. We find that the amplification ratio of the stem increases as the water becomes shallower and this result holds uniformly over several different angles. As the water becomes shallower or as the BL equation tends to the KP equation, the amplification ratios of the BL equation tend to increase to those of the KP equation. For the KP Y-junction solution, the largest amplification ratio is known to be four [5]. For the BL equation, we find that this ratio is reduced by *O*(*ϵ*). An interesting open question is: how much is this ratio reduced in the full water wave equations?

Thus the numerical study yields results consistent with asymptotic arguments and is also useful for predictive purposes. We also numerically solve the BL equation with our chosen dissipation term and show that small shelves of long extent form in the wake of the Y-junction; this is also consistent with previous one-dimensional results (cf. [18]). All of these results are new in the literature.

In summary, in this paper we

— develop a new, conservation law-based perturbation method applicable to any web solution. This significantly simplifies previous, integrable system-based approaches, and thus allows us to model perturbations to the KP equation more readily;

— examine, for the first time, the role dissipation has on web solutions. While we have chosen a particular dissipation model, our methods are, in principle, applicable to any dissipation model; and

— present new numerical results showing how higher order shallow-water effects influence the evolution of the interaction region of web solutions. In the case of Y-junctions, we do this over a large number of parameter values and investigate relative amplification ratios; these results are consistent with KP theory.

While in this paper we only study web solutions in the context of shallow-water flow, we believe that the combination of asymptotic and numerical methods presented here will be useful in other nonlinear problems involving non-decaying profiles.

### (a) The Benney–Luke equation and web solutions of the Kadomtsev–Petviashvili equation

The BL equation is given in a non-dimensional form by [2]
1.1where small dispersion, slow transverse variation and weak nonlinearity are all balanced. These small effects are denoted by *ϵ*,|*ϵ*|≪1; also *q* is the velocity potential, , and , where is related to the ST. As posed, the BL equation has the dispersion relation , where .

In this paper, we only consider small ST, hence , so that . However, when *α*<0, as seen from the dispersion relation, the BL equation is ill-posed. By replacing with , we use a regularized version of the BL equation which is asymptotically equivalent to the above BL equation (1.1) and is written in the form
1.2This regularized BL equation is linearly well-posed; it has the dispersion relationship .

We use the transformation *ξ*=*x*−*t*, *τ*=*ϵt*, so that ∂_{ξ}=∂_{x}, ∂_{t}=−∂_{ξ}+*ϵ*∂_{τ}. Introducing the scaling *A*=(−8*α*)^{1/5}, *β*=2/*A*^{2}, *γ*^{2}=6/*A*^{3}, *δ*=−2/*A*, and letting *q*=(*A*/*β*)*ϕ*(*βξ*,*γy*,*δτ*), (1.2) becomes
Defining the operator [21]
and letting *u*=*ϕ*_{ξ} and , the BL equation we consider is given by
1.3where *F*(*u*) is given by
Hereafter, we replace by *ϵ* (also note that for small ST *A*∼1.2). The definition of is chosen so that if , then is well defined, and
Alternatively, the perturbed BL equation is written as
1.4where .

As seen from the rescaled BL equation (1.3), the KP equation gives the leading-order behaviour of the BL equation. We now introduce the web solutions to KP which serve as the particular leading-order behaviour of interest. A web solution, say *w*(*ξ*,*y*,*τ*) to the KP equation in the Wronskian form is given by (cf. [8]; note for the general N-soliton solution in Hirota form see [3])
1.5where *Ω*(*ξ*,*y*,*τ*)=Wr(*f*_{1},…,*f*_{N}), where Wr denotes the Wronskian of the functions *f*_{i} (cf. [8]). The functions *f*_{i} have the particular form
1.6with *b*_{ij} being a constant, *N*<*M* and , where *k*_{1}<⋯<*k*_{n}. It was proved in [8] that as , there are *N* lines, and as , *M*−*N* lines, in the *ξ*,*y* plane with slopes *c*_{ij}=−(*k*_{i}+*k*_{j})^{−1}, *j*>*i*, such that
1.7where the phases *θ*_{ij} are constant in this pure KP solution.

We now introduce a perturbation scheme for studying the behaviour of solutions *u*(*ξ*,*y*,*τ*) to the BL equation of the form:
where *T*=*ϵτ*. The leading-order term *w* is a web solution to the KP equation in which we allow the terms *k*_{j} to vary in the slow time *T*. Motivated by the form of the web solution, it is convenient to separate the solution *w* into far-field (|*y*|≫1) and near-field () components (see also [17]). The far-field coordinates along a given ray are denoted as *X*_{ij} and *Y* _{ij}, where *j*>*i*, and *X*_{ij} and *Y* _{ij} are determined by the equations
1.8where in the large |*y*| limit derivatives with respect to *Y* _{ij} are considered asymptotically negligible. We also require that *k*_{i}(*T*)+*k*_{j}(*T*)=*C*_{ij}, where *C*_{ij} is independent of *τ* and *T*; hence ∂_{T}*X*_{ij}=0. This requirement means that there is no secular growth in ∂_{T}*w* for large *y*.

We assume along that *u*∼*w*(*X*_{ij},*T*)+*ϵs*(*X*_{ij},*τ*)+⋯ . In this limit, the equation for *w* is,
where *k*_{ij}=*k*_{j}−*k*_{i}. Integrating once in *X*_{ij}, we obtain
which has the soliton solution
1.9and which is consistent with the previous asymptotic arguments. The functions *k*_{ij}=*k*_{ij}(*T*) are functions of the slow time and is an arbitrary phase; these functions are determined later in the perturbation scheme.

At the next order, noting for convenience that we drop the subscripts on *X*_{ij}, the equation for *s*(*X*,*τ*) is
It is also noted that
1.10and
1.11Since is independent of the fast time *τ*, we separate the perturbation *s* into *s*=*h*(*X*)+*g*(*X*,*τ*), so that
1.12with *g*(*X*,0)=−*h*(*X*).

The slowly varying phases were determined in [17] via the additional orthogonality condition 1.13which removed unbounded growth as . We note that orthogonality conditions such as these usually correspond to physical constraints, i.e. conservation laws (cf. [18]).

In §2, from the conservation-law approach a direct and physically motivated derivation of the condition (1.13) is given; also in §2 the role of dissipation is examined. In §3, we present our numerical findings concerning the evolution of several different types of web solutions. These results indicate that web solutions are robust to these perturbative effects over asymptotically significant lengths of time. The numerical study of the amplification of the stem of a Y-Junction in the BL equation, and the effect of dissipation, is also presented in §3. The electronic supplementary material collects some of the technical details of the approach used in §§2 and 3.

## 2. Conservation-law scheme for arbitrary web solutions

In order to determine the evolution of the slow coefficients, we use a method [18], which employs the conserved quantities of the leading-order integrable problem to determine the evolution of slow variables. We use conservation of energy which is defined by the following iterated integral:
The average along the *y*-axis is taken because the line solutions to the KP equation do not decay along the *y*-axis. It is also important to point out that for web solutions and their perturbations that
and so the placement of the limit is important.

Using the above definition of energy, it follows that
The analytical issues concerned with differentiating the energy functional are discussed in the electronic supplementary material. From the definition of *K*(*u*), we expect that
2.1which would then give
2.2Note, owing to the average, in order to show (2.1), we must establish that
2.3To show this, some care must be taken to accommodate for the presence of the average which makes integration by parts more complicated. This is worked out in the electronic supplementary material for web solutions to the KP equation. We further assume that it holds for the perturbation *u*=*w*+*ϵs*+⋯ , so that taking (2.1) to be true is valid to the order of the asymptotics.

Since we are assuming *u*(*ξ*,*y*,*τ*) is, in particular regions of the plane, a solution of the form:
2.4and is otherwise at most outside the regions of interest, the ansatz for u implies *w*_{τ}=*ϵw*_{T}. Furthermore, based on the one-dimensional problem, it is natural to assume that
2.5(this is discussed in more detail below), then we obtain
2.6

We now demonstrate, in some detail, how to use the solvability conditions associated with (1.12) in order to reduce (2.2) to (2.5) and (2.6). As , using the coordinate change given by (1.8), *w* depends only on *X*_{ij}, hence becoming stationary with respect to *τ*. The perturbation *s*(*X*_{ij},*τ*) is written as *s*(*X*_{ij},*τ*)=*h*(*X*_{ij})+*g*(*X*_{ij},*τ*) in region *R*_{ij} (figure 1), where *h* solves the following equation:
2.7

Note, (2.7) is the same as (1.12). Then in *R*_{ij}, where we ignore the *Y* _{ij} coordinate and thus reduce the problem to one dimension, i.e. *X*_{ij}, as shown in the introduction, we have (*K*′(*w*))^{†}*w*=0. In order to ensure that a solution *h* exists to (2.7), the Fredholm alternative enforces to leading order that
2.8This result gives us immediately that the terms *k*_{ij} are independent of the slow time *T*. To show this, we first note that
2.9since *w* is even in *X*, and from (1.10), so is . Using (2.8) then shows that
On the other hand, using (1.11), we have
Since , this implies that ∂_{T}*k*_{ij}=0. Having now shown that the coefficients *k*_{ij} do not vary slowly, and given our assumption that *k*_{j}+*k*_{i} does not depend on *T*, evidently the coefficients *k*_{j} are also independent of the slow time *T*.

We now present an argument that establishes (2.6) using (2.8). First, we have the identity
Taking *L*_{a}≫1, and letting , we can write the integrals
as sums over the regions *R*_{ij} since in between the regions, *w* and its derivatives are exponentially small. Then, over any region *R*_{ij}, we have
where *J*_{ij}=1+(*k*_{i}+*k*_{j})^{−2} is the Jacobian of the coordinate transformation and *Y* _{a} is some constant value of *Y* _{ij} that only depends on the choice of *L*_{a}. Using (2.8) then establishes that
Given that
since *L*_{a} is constant, we finally have
or that (2.6) holds.

We now turn to determining the phases . To do this, using (2.6), (2.2) becomes, noting also that *s*(*ξ*,*y*,0)=0,
2.10This global relation can be separated into near- and far-field components; note the integrals in the interaction region [−*L*_{a},*L*_{a}] are negligible in comparison with the far-field. Then in the far-field, we can write this integral over regions *R*_{ij} (again see figure 1 for clarification). If we assume that after a short time the dominant part of *s*(*X*_{ij},*τ*) is given by *h*(*X*_{ij}), then since *k*_{j} is constant , so that (1.12) becomes
One obtains from the global condition (2.10) that in *R*_{ij}, again ignoring *Y* _{ij},
2.11

This is exactly the orthogonality condition derived in [17] (i.e. (1.13)) that was used to determine the evolution of the difference between phases. To compute this phase evolution, first introduce the transformation , and then let , which turns (1.12) after one integration in into
where the coefficients , *p*=1,2,3, which implicitly depend on the region of interest *R*_{i,j}, are given by

This leads to the solution
Using (2.11) gives or
2.12With this, we can now compute the speed, *v*_{i,j}, of each part of the far-field via the formula
2.13where ∂_{τ}*X*_{ij} is given in (1.8).

Therefore, using a conservation-law approach, we have seen how the results of [17] can be derived in a way that makes no use of integrable systems. This makes the perturbation method presented in this paper more broadly applicable than that presented in [17].

### (a) Dissipation

The BL equation conserves energy, and is a Hamiltonian system [14]. Of course, real physical systems such as water waves are not perfectly conservative. So it is useful and instructive to study the effect of adding typical damping terms to the conservative model. While the dissipative model we use is special, it is physically relevant [19], and it shows how to deal with methodological difficulties created by dissipation models. It is also instructive to study this dissipation model since it causes the formation of small-amplitude shelves with long extent.

For the BL equation, we add a linear local term which leads us to study an equation of the form
where *F*(*u*) represents dispersive terms and the constant *γ*>0 represents the magnitude of damping. We can then define the function and repeat the analysis from above. Matching terms of after expanding gives
and repeating the Fredholm alternative argument that led to (2.6), yields
2.14Now we separate (2.14) over the regions *R*_{ij}. We showed in the previous section, i.e. equation (2.9), that

Then, if we naively use the asymptotic form of the solution given by (1.9) and the corresponding representation for *w*_{T} in *R*_{ij}, we obtain
so that we have
Coupling this condition with the requirement that *k*_{i}(*T*)+*k*_{j}(*T*)=*C*_{ij} creates a difficulty. For example, if we have *k*_{1}, *k*_{2} and *k*_{3}, then a web solution the slopes of which are determined by say *k*_{1} and *k*_{2}, and *k*_{1} and *k*_{3}, would have two different values for *k*_{3}(*T*). Thus, we see that the ansatz given by equation (1.9) for the leading order behaviour is not appropriate globally.

To deal with these inconsistencies, we use a modified ansatz for *w* (the index *i*,*j* on *w* and *ψ*_{l}, *l*=1,2, is understood) of the form:
where , and
with
Note, the functions *η*(*T*), *ψ*_{1}(*T*) and *ψ*_{2}(*T*) vary in each region *R*_{ij}.

To determine , again noting that *s*(*ξ*,*y*,0)=0, we employ the global relationship
To make use of this identity, we first assume that on *R*_{ij} that the dominant contribution to *s* is given by the solution to the stationary equation
2.15which is the same equation as (1.12) except now taking the dissipation into account. Using the transformations from the previous section, we solve for *h* and derive the phase equation

## 3. Computation of the Benney–Luke equation

We first note that the BL equation is second order in time. It is convenient for the numerics to transform to a system by introducing the variable *v*=*u*_{τ}. In order to simulate the BL equation numerically, we use a windowing method developed in [22]. This method introduces a smooth function of compact support, say *V* (*y*), such that *V* (*y*) is nearly one over some interval in *y*, say [−*L*_{y}+*δ*,*L*_{y}−*δ*], and *V* has support in [−*L*_{y},*L*_{y}]. We use the function
where eps is of the order of 10^{−16} (approximately machine precision) and write the solution to the BL equation as
We use the leading-order asymptotic solution computed in the previous sections to evaluate *u* and *v* in the far-field with distances greater than for times . We denote these far-field solutions as *u*_{asy} and *v*_{asy}.

We now wish to determine *V* (*y*)*u*. Define the functions *u*_{nr}=*V* *u*, *v*_{nr}=*V* *v*, so that
Using the fact that *u*_{asy} satisfies the KP equation in the far-field, and keeping only terms through *O*(*ϵ*), by substituting *u* and *v* into the BL equation, it then follows that *u*_{nr} satisfies the equation
where denotes the BL equation. The term *F*_{f} is given by
which has support only in the intervals [−*L*_{y},−*L*_{y}+*δ*] and [*L*_{y}−*δ*,*L*_{y}]. Note, we ignore the window effects from the nonlinearities of since these terms are smaller and isolated to have support in the region [−*L*_{y},−*L*_{y}+*δ*] and [*L*_{y}−*δ*,*L*_{y}]. If we work on a large enough domain, the effect of the error should be nominal far from the boundary since the error propagates with finite speed.

The functions *u*_{nr} and *v*_{nr} satisfy at *τ*=0 the boundary conditions
At later times, we enforce periodic boundary conditions in the *y*-coordinate and choose the window large enough to ensure that
Given the exponential decay in *ξ*, we also restrict the domain in *ξ* to the interval [−*L*_{ξ},*L*_{ξ}] and enforce periodic boundary conditions in *ξ* and *y* in order to numerically solve for *u*_{nr}. This gives us the advantage of being able to use pseudo-spectral methods.

However, we have to find the operator in terms of Fourier series. Writing a typical periodic function, say , as
we calculate the inverse derivative of via
3.1where we have used . In the electronic supplementary material, we show that the pseudo-spectral representation of the (*j*,*l*) mode of the inverse derivative is given by
3.2where and *N*_{T} is the number of modes used in the pseudo-spectral approximation.

### (a) Evolution of the Y-junction and X-wave in the Benney–Luke equation

The following figures show top-down surface plots of the numerical evolution of various web-solution profiles. The time-stepping algorithm used was the ETDRK4 method [23], which is a variant of the fourth-order Runge–Kutta method.

The figures represent numerical simulations run on time scales long enough to allow the next order terms that distinguish the BL equation from the KP equation to have a significant impact. This effect manifests itself in weak dispersive tails; see figure 3 for clarification. In each figure, the initial conditions were chosen as *w*(*ξ*,*y*,0) and *w*_{τ}(*ξ*,*y*,0) as defined earlier, see (1.5) and (1.6). The *X*-wave solutions in figure 2*a*–*d* correspond to the classical two-soliton solutions discussed in [1,5]. In the language used earlier, we choose the matrix of coefficients *B* with entries *b*_{ij} as (see (1.5) and (1.6)) follows:
This gives the function *Ω*(*ξ*,*y*,*τ*) as
We then choose *b*_{12}=(*k*_{4}−*k*_{1})/(*k*_{4}−*k*_{2}),*b*_{24}=(*k*_{3}−*k*_{2})/(*k*_{4}−*k*_{2}), and we set *k*_{1}=−*k*_{4}, *k*_{2}=−*k*_{3}, *k*_{3}>0 and *k*_{4}=1+*k*_{3}. Given this choice of coefficients *b*_{ij} and *k*_{l}, by letting , one can let the length of the stem grow as shown in figure 2*c*,*d*.

A different case can be investigated by setting *k*_{3}=*k*_{2}. Then *Ω* becomes
If we then choose *b*_{24}=(*k*_{2}−*k*_{1})/(*k*_{4}−*k*_{1}), *b*_{12}=(*k*_{4}−*k*_{1})/(*k*_{4}−*k*_{2}), and relabel *θ*_{4} as *θ*_{3}, we then get , where . We note that , and so from this point on we use as a solution to the KP equation. The function corresponds to the *B* matrix
which shows that the *X*-wave has passed to a *M*=3, *N*=2 web solution. This is an example of a Y-junction. We choose the parameters *k*_{j} such that
3.3where 0≤*k*≤1. Alternatively, we can use the formalism in [3] to obtain these initial conditions.

These Y-junctions have a ray along with maximum amplitude , another ray along the angle arctan(−*k*) of amplitude , and finally a ray of amplitude *k*^{2}/2 moving at an angle of *π*/4 where the angles are measured relative to the *y*-axis. Thus, as , the KP equation predicts a fourfold amplification of the ray along relative to the amplitude of the ray along the angle arctan(−*k*). Likewise, as we vary *k*, we move between a one-dimensional line soliton solution when *k*=0 up to a Y-junction with maximal amplification ratio of a factor four when *k*=1. As to the effect of the BL equation, we refer to figure 3 in which one sees that the radiation is more significant in figure 3*a*, where *ϵ*=0.2, than in figure 3*b*, where *ϵ*=0.1. We note that the decrease in the magnitude of the dispersive tails as *ϵ* decreases is a feature found in the X junctions as well.

In figure 2, it can be seen that the primary wave remains essentially unchanged over time scales that are the reciprocal of the order of magnitude of the perturbation. We also point out that the figures show that as *ϵ* decreases, or as the water becomes shallower, the amount of radiation decreases. These figures indicate that the web solutions studied here are stable with respect to evolution in the BL equation. Hence, this result provides quantitative evidence that the class of web solutions examined in this paper are robust with respect to the perturbations introduced by the BL equation.

### (b) Amplification ratio for the Y-junction in the Benney–Luke equation

The problem of the amplification ratio of the stem in a Y-junction is an interesting issue that was first studied in [5] in the context of the KP equation. We define the amplification ratio to be the magnitude of the ray along and the ray along the angle arctan(−*k*), which for a KP solution is given by (1+*k*)^{2}. We now discuss how the BL equation affects the amplification ratio, and how this class of KP solutions behaves under the influence of the BL equation.

As indicated above, taking a Y-junction with parameter values given by (3.3) and its corresponding derivative with respect to time as Cauchy data, figure 4 shows the amplification ratio (vertical axis) of a Y-junction in the BL equation for different values of *ϵ* and *k*. For each curve, the simulations were run for times up to *τ*=1/*ϵ*, thus allowing time for the BL equation to have an asymptotically significant impact. The domain size was chosen to be *L*_{ξ}=50 and *L*_{y}=48 with 512 modes used for the pseudo-spectral approximation. The size of the rays of the Y-junction were measured at *L*_{y}= ± 30 so as to avoid any effects from the windowing. In order to find the wave amplitude, we averaged the height of the wave over seven mesh points around the grid point corresponding to *L*_{y}= ± 30; this minimizes the impact of any spurious oscillations because of the numerical method. Finally, we chose *k*=0.1,…,1 in intervals of a tenth. Figure 2*e*,*f* are plots of the case *k*=1 for *ϵ*=0.2 and *ϵ*=0.1, respectively.

We also note that the phase, or speed correction, along the ray is given by so that the speed, see (2.13), of the ray in the far field should be −(1+*ϵ*((1+*k*)^{4}/96)), where the negative sign indicates that the ray moves to the left. While this correction goes into the windowing approximation, as a consistency check, we also numerically computed the speed of the ray at *y*=−30 for *ϵ*=0.1 and *k*=1. After averaging to take into account errors from discretization, we got a computed speed *of*−1.0090, while the predicted speed from the asymptotics *is*−1.0167, so the agreement is very good.

Figure 4 indicates that the numerically generated solutions tend to the exact results for the KP equation as *ϵ* decreases and shows that increasing *ϵ* decreases the amplification ratio for each *k*. As indicated above, we note that larger values of *ϵ* introduce more dispersive radiation. This also implies lower pulse heights due to energy conservation which in turn suggests that there will be a somewhat smaller amplification ratio than four at *k*=1 of *O*(*ϵ*) (i.e. about 3.9 for *ϵ*=0.1), but still much larger than linear theory. This result also agrees with the fact that the BL equation contains the KP equation as an asymptotic limit.

The variation of the amplification as a function of *ϵ* holds uniformly for all *k*. This indicates that the numerics is well behaved as we vary the shape of the initial conditions. Therefore, we find that the combination of asymptotic evaluation of the far-field, pseudo-spectral method, ETDRK4 and windowing is an effective way of approaching problems with non-decaying solutions.

### (c) Dissipation in the Y-junction

We next provide numerical results corresponding to the linear dissipation model discussed earlier. As a leading-order solution, we use the Y-junction from the previous section with *k*=1 so that *k*_{13}=2, and we set the coefficient of dissipation *γ*=2. However, we also take into account the dissipation in the far-field using the asymptotic method shown earlier; this is how we find *u*_{asy} and *v*_{asy}. The size of the domain is *L*_{ξ}=50 and *L*_{y}=48, with 512 modes used in the pseudo-spectral method. We measure the maximum height of the numerical solution in the *region*−30≤*y*≤30 in order to avoid any effects from the windowing. The plot in figure 5*a* gives a log plot of the amplitude decay.

Thus, at *τ*=1/*ϵ*, the asymptotic theory predicts that the maximum amplitude of the ray should be 2×10^{−0.66}=1.033, while the numerics gives the maximum amplitude as 2×10^{−0.55}=1.15. Thus, the numerics confirm the asymptotic prediction since the discrepancy between the two results is . We look at a cross section of the ray at *y*=−30 in figure 5*b*. As shown in the figure, to the right of the peak a shelf forms. This is consistent with the dissipative theory of KdV for linear dissipative models of wave amplitude [18]; this is another novel feature of web solutions not yet reported in the literature.

## Acknowledgements

This research was partially supported by the National Science Foundation under grant no. DMS-0905779. We thank Douglas Baldwin for many stimulating and enlightening conversations.

- Received November 20, 2012.
- Accepted January 18, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.