Coupling Conditions for Water Waves at Forks

Abstract

We considered the propagation of nonlinear shallow water waves in a narrow channel presenting a fork. We aimed at computing the coupling conditions for a 1D effective model, using 2D simulations and an analysis based on the conservation laws. For small amplitudes, this analysis justifies the well-known Stoker interface conditions, so that the coupling does not depend on the angle of the fork. We also find this in the numerical solution. Large amplitude solutions in a symmetric fork also tend to follow Stoker’s relations, due to the symmetry constraint. For non symmetric forks, 2D effects dominate so that it is necessary to understand the flow inside the fork. However, even then, conservation laws give some insight in the dynamics.

1. Introduction

The propagation of nonlinear waves in a network is an important topic. As an example, consider a hydrological network which is prone to floods. Understanding the global dynamics of the network can help identify its most vulnerable sections and take the appropriate measures. Real networks are formed by long 2D or 3D channels of a small cross-section. To study the propagation of waves in such systems, a first step is to consider a simple fork as a model of elementary junctions. The final goal is to reduce the model to 1D channels connected by appropriate interface conditions. The study of such 1D systems is now well advanced, in particular for systems of conservation laws, see the review [1].

The type of PDE model describing the quantity propagating on the network is very important to derive the coupling conditions. Recently for the sine-Gordon nonlinear wave equation, we [2] introduced a homothetic reduction [3] where we averaged the operator over the fork region and consistently took the limit when the width tended to zero. Assuming continuity of the field, we obtained Kirchhoff’s law for the gradients. Comparing the 2D solution with the one for the reduced 1D equations gives excellent agreement. In this situation, the angle of the fork does not play a role. When considering networks of rivers, many authors, for example Stoker [4] and Jacovkis [5] assumed continuity of the water height and continuity of the flux so that again, the angle of the fork did not come in. In the close context of gas dynamics, Holden and Risebro [6] studied shocks in a pipe with an elbow. They showed that the Riemann problem had a unique solution when the angle was smaller than $\pi$. The angle is also important for classical hydrodynamics; in a fork, it sets the forces experienced by the pipes [7]. In fact, for large amplitude shallow water waves our numerical calculations show that the energy entering a branch can vary from 20% to 50% depending on the symmetry of the fork. These studies point out the importance of the angle.

A few authors addressed the problem of the angle of a junction. Schmidt [8] studied the 2D connection between 1D channels; he made no assumption on the size of the connecting domain. The flow in the junction was assumed linear so that the author used a variational method that gave the solution as a superposition of fields. The final result was a system of ordinary differential equations for the values at the ends of the branches coupled to the shallow water PDEs. Despite its formal beauty, it remains difficult to handle and does not give a simple picture. Shi et al. [9] studied experimentally and numerically the propagation of long waves in wide and narrow channels. They used the Boussinesq dispersive shallow water equations for narrow channels. They observed no angle dependence and a strong transmission. For the same equations, Nachbin and Simoes [10] obtained interface conditions containing implicitly the angles of the fork. These gave an excellent matching between the average of the 2D solution and the solution of the 1D effective model for angles smaller than $\pi /3$.

In this article, we consider the nonlinear shallow water equations. The system is very general because it only involves conservation laws. Also it is simple enough. We revisit the problem of shallow water propagation in 2D forks using our homothetic reduction procedure to obtain approximate conservation laws and compare them with the numerical solutions. We compute approximate conservation for the mass, momenta and energy laws for a general fork geometry. In the small amplitude limit we recover Stoker’s conditions, i.e., continuity of surface elevation and mass conservation (Kirchoff law). To our knowledge, this is a first formal justification of Stoker’s interface conditions. This angle independent reduction holds also for a general class of scalar nonlinear wave equations, for example the 2D sine-Gordon equation or the 2D reaction-diffusion equation; it confirms the results of [2]. We computed the 2D numerical solution for a simple T-fork geometry for small and large amplitudes. The wave was also launched in two different branches to see the effect of symmetry. We show that Stoker’s conditions hold for the symmetric case for small and large amplitudes. For the non-symmetric case, they hold for small amplitudes. When the amplitude is large, 2D effects dominate the fork region. Nevertheless the approximate conservation laws give an insight into the flow.

The article is organized as follows. Section 2 presents the fork geometry and shows the straightforward reduction for a general class of nonlinear wave equations. In Section 3 we recall the shallow water equations and their conserved quantities. Section 4 gives the integrals of these equations on the fork showing that the mass and energy laws do not involve the angles while the momenta laws do. Section 5 shows the 2D numerical solutions for symmetric and non symmetric configurations for small and large waves. There, we compare the numerical results with the conservation laws established in Section 4. We discuss these results and conclude in Section 6.

2. General Scalar Nonlinear Wave Equations

Before considering the nonlinear shallow water equations, we analyze the simpler case of a class of scalar 2D nonlinear wave equations. This large class includes hyperbolic wave equations like the sine-Gordon equation as well as reaction diffusion equations like the Fisher equation, to name a few. We consider equations of the form

$$\alpha u_{tt}+\beta u_t-\Delta u=N\left(u\right),$$

(1)

where $u(x,y,t)$ is a scalar, $\Delta$ is the usual 2D Laplacian and where $N\left(u\right)$ is a nonlinearity not containing derivatives. The boundary condition on the lateral domain is of Neumann type

$${\partial }_{n}u=\nabla u·\mathbf{n}=0.$$

(2)

Consider the fork domain shown in Figure 1. Far from the fork region, the solution can be assumed to be 1D so that we do not loose much information by approximating the 2D dynamics with a 1D equation. Inside the fork domain, a strong coupling occurs between the branches. To see this, we proceed as in [2] and integrate the operators on the fork region. Then we examine the behavior of the different terms as w, the width of the branches, goes to zero. We assume that domains that we consider behave in a regular way as we shrink w homothetically to zero, [3].

Consider the asymmetric Y-branch shown in the left panel of Figure 1. A first assumption is the continuity of u which is obvious for the 2D operator. The other condition comes from the integration of the operator (1) on the fork domain $\mathcal{F}=IABCDEFGHI$. We get

$$\int [\alpha u_{tt}+\beta u_t-N\left(u\right)]dxdy-{\int }_{\partial \mathcal{F}}(\nabla u)·\mathbf{n}ds=0.$$

(3)

The first integral is of order $O\left(w^2\right)$. On the exterior boundaries, $(\nabla u)·\mathbf{n}=0$ so the line integral reduces to

$${\int }_{IA}\cdots +{\int }_{CD}\cdots +{\int }_{FG}\cdots ,$$

which are $O\left(w\right)$. We then obtain for $w\to 0$

$$-{\partial }_s{u}_1+{\partial }_s{u}_2+{\partial }_s{u}_3=0,$$

(4)

where $u_i,i=1,2,3$ are respectively the values of the field at the end of branch 1 (IA) and at the beginning of branches 2 (FG) and 3 (CD). Relation (4) is Kirchhoff’s law [2]. When the widths of the branches are not equal, this Kirchoff relation becomes

$$-w_1{\partial }_s{u}_1+w_2{\partial }_s{u}_2+w_3{\partial }_s{u}_3=0.$$

(5)

Remark that in the result (4) the angle of the fork plays no role. The reduction leading from the flux equation to (5) is an asymptotic result that holds for $w\to 0$. It is then natural to approximate the 2D Equation (1) by a 1D equation in each branch together with the conditions of continuity and Kirchoff (4) at the junctions.

The result we obtain can be connected to a property of the Laplace operator with Neumann boundary conditions on a so-called “fat” graph [11]. Consider a graph where each edge has a transverse size w, assume Neumann boundary conditions on the transverse edge. Then the spectrum of the Laplacian converges to the one of the 1D Laplacian as $w\to 0$. This is true for compact and non compact graphs. See the article by Exner and Post [11] and the book by Post [12] for the details of the proof.

The validity of the reduction was confirmed numerically for the 2D sine-Gordon equation, (1) with $\alpha =1,\beta =0$ and $N\left(u\right)=-sin\left(u\right)$ in [2]. There we compared the 2D solutions to the ones of the 1D sine-Gordon equation in each branch, coupled by the interface conditions. For completeness, we recall the case of a sine-Gordon kink propagating in forks with angles 45 and 90 degrees. The kink is an exact solution in 1D, it is

$$u(x,t)=4arctan\left[exp(\frac{x-vt}{\sqrt{1-v^2}})\right],$$

(6)

where the velocity $0\le v<1$ is a free parameter. To compare the 2D and 1D solutions, we plot the energies in each branch

$$E_2^i={\int }_{{\Omega }_i}\left[\frac{1}{2}{u}_t^2+\frac{1}{2}{|\nabla u|}^2+(1-cosu)\right]dxdy,$$

(7)

and

$$E_1^i=\sum _{i=1,2,3}{\int }_{{\Omega }_i}\left[\frac{1}{2}{u}_t^2+\frac{1}{2}{\left|u_x\right|}^2+(1-cosu)\right]dx,$$

(8)

where ${\Omega }_i$ is branch i, abusively named the same in 1D and 2D. The kink is started in branch 1 with an initial velocity $v=0.75$, this gives a typical wavelength $\lambda \approx 4/\sqrt{1-v^2}=2.7$. The width of the branches is $w=0.7<<\lambda$. Figure 2 shows the time evolution of the energies $E_2^i$ for forks with angles 45 and 90 degrees and $E_1^i$, where $i=1,2$ corresponds to the branches. Initially the kink is in branch 1 so that $E_2^2=E_2^3=0$. As the kink crosses into branches 2 and 3, $E_2^1$ becomes very small. Note the excellent agreement between the two expressions $E_2^i$ and the expression $E_1^i$. This confirms that the angle of the fork plays no role for such a system.

The dynamics of kinks for the sine-Gordon equation is controlled by the energy: if the initial energy is enough, a kink in branch 1 gives rise to two kinks in branches 2 and 3. This gives a very simple picture. Other solutions like the breather have much more complicated dynamics, we refer the reader to [2] for more details.

The dynamics of such waves can then be studied for general networks as we have done in [13].

3. The Nonlinear Shallow Water Equations

The shallow water equations in a 2D domain written in terms of the fluid velocity $\mathbf{u}(x,t)$

$$\mathbf{u}={(u,v)}^T$$

and the water height $h(x,t)$ read [4]

$$\begin{array}{c}\hfill h_t+\nabla ·\left(h\mathbf{u}\right)=0,\end{array}$$

(9)

$$\begin{array}{c}\hfill {\left(hu\right)}_t+\nabla ·\left(\begin{array}{c}h{u}^2+\frac{g{h}^2}{2}\\ huv\end{array}\right)=0,\end{array}$$

(10)

$$\begin{array}{c}\hfill {\left(hv\right)}_t+\nabla ·\left(\begin{array}{c}huv\\ h{v}^2+\frac{g{h}^2}{2}\end{array}\right)=0,\end{array}$$

(11)

where g is the gravitational acceleration. The wall boundary condition is

$$\mathbf{u}·\mathbf{n}=0.$$

(12)

We assume an even bottom of the channels $h=h_0$.

3.1. Conserved Quantities

We first recall the conserved quantities. Integrating Equations (9)–(11) over a 2D closed domain $\Omega$ and using the boundary condition (12) we get

$${\partial }_t{\int }_{\Omega }hdxdy=0,$$

(13)

$${\partial }_t{\int }_{\Omega }hudxdy+{\oint }_{\partial \Omega }\frac{g{h}^2}{2}{n}_{x}ds=0,$$

(14)

$${\partial }_t{\int }_{\Omega }hvdxdy+{\oint }_{\partial \Omega }\frac{g{h}^2}{2}{n}_{y}ds=0.$$

(15)

A localized wave will have as first conserved quantity the integral of the water elevation

$$M={\int }_{\Omega }hdxdy.$$

The total x and y momenta

$$P_x={\int }_{\Omega }hudxdy,P_y={\int }_{\Omega }hvdxdy$$

will not be conserved in the fork geometries.

A flux relation that can be deduced from the conservation laws (9)–(11) is the total energy flux

$$e_t+\nabla ·\left[\mathbf{u}(e+\frac{g{h}^2}{2})\right]=0.$$

(16)

where the total energy density is

$$e=\frac{1}{2}\left[g{h}^2+(u^2+v^2)h\right].$$

(17)

Integrating the energy flux relation over a volume $\Omega$ we obtain that a localized wave in $\Omega$ will have constant energy

$$\frac{dE}{dt}=\frac{d}{dt}{\int }_{\Omega }edxdy=0.$$

3.2. Small Amplitude Limit

It is well known that in the linear limit, Equations (9)–(11) reduce to the linear wave equation for the water height h. To see this, consider the steady state $h=h_0,u=v=0$, then the linearized system is

$$\begin{array}{c}\hfill h_t+h_0\nabla \mathbf{u}=0,\end{array}$$

(18)

$$\begin{array}{c}\hfill h_0{\mathbf{u}}_t+g\nabla h=0,\end{array}$$

(19)

Taking the time derivative of the first equation and plugging in the second equation, we get the wave equation

$$h_{tt}-g{h}_0(h_{xx}+h_{yy})=0.$$

(20)

The boundary conditions reduce to $\nabla h·\mathbf{n}=0$ as can be seen by projecting (19) on $\mathbf{n}$. This equation is in the class (1).

4. Reduction of the Shallow Water Equations

The shallow water equations cannot be reduced so simply as the nonlinear scalar wave equation. In fact, it is not clear what are the right interface conditions that should be implemented for a 1D effective model. Stoker, in his well-known book [4] introduces the following interface conditions for the water elevations $h_1,h_2,h_3$ and branch-oriented velocities $u_1^{\parallel },u_2^{\parallel },u_3^{\parallel }$

$$\begin{array}{c}\hfill h_1=h_2=h_3,\end{array}$$

(21)

$$\begin{array}{c}\hfill -h_1{u}_1^{\parallel }+h_2{u}_2^{\parallel }+h_3{u}_3^{\parallel }=0,\end{array}$$

(22)

and uses them to analyze the junction of the Mississippi and the Missouri rivers. These conditions were not justified by a formal argument. Note also that they do not depend on the angle of the junction.

Below, we will see that these conditions arise naturally in the limit of small amplitude for the shallow water equations. For general amplitudes, it is not clear that these apply. To analyze the problem, we proceed as in [2], integrate the governing equations on the bifurcation region and consider the limit of vanishing transverse width w.

4.1. Mass Flux

Integrating the Equation (9) over the closed region $\mathcal{F}\equiv ABCDEFGHIA$ yields

$${\int }_{\mathcal{F}}{h}_{t}dxdy+{\oint }_{\partial \mathcal{F}}h\mathbf{u}·\mathbf{n}ds=0.$$

Because of the boundary condition $\mathbf{u}·\mathbf{n}=0$ on ABC, DEF and GHI the expression above reduces to

$${\int }_{\mathcal{F}}{h}_{t}dxdy+{\int }_{AI}h\mathbf{u}·\mathbf{n}ds+{\int }_{CD}h\mathbf{u}·\mathbf{n}ds+{\int }_{FG}h\mathbf{u}·\mathbf{n}ds=0.$$

The first integral is $O\left(w^2\right)$ while the three other integrals are $O\left(w\right)$. Dividing the equation by w and taking the limit $w\to 0$ we get from these three terms

$$-h_1{u}_1^{\parallel }+h_2{u}_2^{\parallel }+h_3{u}_3^{\parallel }=0,$$

(23)

where we have introduced the local branch-oriented velocities $u^{\parallel },u^{\perp }$ such that

$$\left(\begin{array}{c}{u}^{\parallel }\\ u^{\perp }\end{array}\right)=\left(\begin{array}{cc}cos\theta & sin\theta \\ -sin\theta & cos\theta \end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)$$

(24)

and where the indices 1,2 and 3 refer to the branches. Of course, when the transverse widths $w_1,w_2,w_3$ are different, with the condition that the ratios $w_2/w_1,w_3/w_1$ remain finite, the relation (23) becomes

$$-w_1{h}_1{u}_1^{\parallel }+w_2{h}_2{u}_2^{\parallel }+w_3{h}_3{u}_3^{\parallel }=0.$$

4.2. Energy Flux

The energy flux (16) can be consistently reduced to a 1D relation. As for the mass relation, we integrate Equation (17) over the domain $\mathcal{F}=$ ABCDEFGHIA to obtain

$${\int }_{\mathcal{F}}{e}_{t}dxdy+{\oint }_{\partial \mathcal{F}}(e+\frac{g{h}^2}{2})\mathbf{u}·\mathbf{n}ds=0.$$

Because of the boundary condition $\mathbf{u}·\mathbf{n}=0$ on ABE, the expression above reduces to

$${\int }_{\mathcal{F}}{e}_{t}dxdy+{\int }_{AI}(e+\frac{g{h}^2}{2})\mathbf{u}·\mathbf{n}ds+{\int }_{CD}(e+\frac{g{h}^2}{2})\mathbf{u}·\mathbf{n}ds+{\int }_{FG}(e+\frac{g{h}^2}{2})\mathbf{u}·\mathbf{n}ds=0.$$

The first integral is $O\left(w^2\right)$ while the three other integrals are $O\left(w\right)$. Dividing the equation by w and taking the limit $w\to 0$ we get from these three terms

$$-(e_1+\frac{g{h}_1^2}{2})u_1^{\parallel }+(e_2+\frac{g{h}_2^2}{2})u_2^{\parallel }+(e_3+\frac{g{h}_3^2}{2})u_3^{\parallel }=0.$$

(25)

To conclude, Equation (9) gives in the 1D limit, the balance of mass (23). The same happens for the energy flux (16) which yields (25). The natural matching conditions for 1D shallow water equations on a network are then

$$-h_1{u}_1^{\parallel }+h_2{u}_2^{\parallel }+h_3{u}_3^{\parallel }=0,$$

(26)

$$-u_1^{\parallel }(g{h}_1^2+h_1\frac{{u_1^{\parallel }}^2}{2})+u_2^{\parallel }(g{h}_2^2+h_2\frac{{u_2^{\parallel }}^2}{2})+u_3^{\parallel }(g{h}_3^2+h_3\frac{{u_3^{\parallel }}^2}{2})=0.$$

(27)

For the mass and the energy balance laws, we have a similar situation to the one of the nonlinear scalar wave equation, the angles of the fork do not play any role. In the small amplitude limit, the speeds $u_1,u_2,u_3$ are small and the squares can be neglected in the energy relation. Then, we recover the Stoker interface conditions (21).

4.3. Momentum Flux for a General Fork

Contrary to the mass and the energy, the momentum Equations (10) and (11) cannot be consistently reduced to a 1D condition involving $h,u^{\parallel }$ at each end of $\mathcal{F}$.

To see this, integrate the horizontal momentum Equation (10) over the domain $\mathcal{F}$ and get

$${\int }_{\mathcal{F}}{\left(hu\right)}_{t}dxdy+{\oint }_{\partial \mathcal{F}}\left(\begin{array}{c}h{u}^2+\frac{g{h}^2}{2}\\ huv\end{array}\right)·\mathbf{n}ds=0,$$

where the first integral is a surface integral and the second one a line integral. In the integrand of the latter, we have

$$\left(\begin{array}{c}h{u}^2\\ huv\end{array}\right)·\mathbf{n}=hu\left(\begin{array}{c}u\\ v\end{array}\right)·\mathbf{n}=0$$

on the exterior boundaries of $\partial \mathcal{F}$ because of the boundary condition (2). Then, only the potential term $\frac{g{h}^2}{2}$ will contribute to these terms.

The $O\left(w\right)$ terms (line integrals) reduce to

$$\begin{array}{c}-\frac{g}{2}(|AB|h_{AB}^2-\left|HI\right|h_{HI}^2)-sin{\theta }_2\frac{g}{2}(|BC|h_{BC}^2-\left|DE\right|h_{DE}^2)-sin{\theta }_3\frac{g}{2}(|EF|h_{EF}^2-\left|HG\right|h_{HG}^2)\\ -w{h}_1{u}_1{v}_1+w\left[(h_2{u}_2^2+g\frac{h_2^2}{2})cos{\theta }_2+h_2{u}_2{v}_{2}sin{\theta }_2\right]+w\left[(h_3{u}_3^2+g\frac{h_3^2}{2})cos{\theta }_3+h_3{u}_3{v}_{3}sin{\theta }_3\right]=0.\end{array}$$

(28)

Using the branch oriented velocities (24) we get the approximate law

$$\begin{array}{c}-\frac{g}{2}\left(\right|AB|h_{AB}^2-\left|HI\right|h_{HI}^2)-sin{\theta }_2\frac{g}{2}\left(\right|BC|h_{BC}^2-\left|DE\right|h_{DE}^2)-sin{\theta }_3\frac{g}{2}\left(\right|EF|h_{EF}^2-\left|GH\right|h_{GH}^2)\\ -w{h}_1{u}_1{v}_1+wcos{\theta }_2\left[h_2{u_2^{\parallel }}^2+g\frac{{h_2}^2}{2}\right]+wcos{\theta }_3\left[h_3{u_3^{\parallel }}^2+g\frac{{h_3}^2}{2}\right]=0,\end{array}$$

(29)

where we neglected the velocity components $u^{\perp }$.

Similarly for the vertical momentum equation we obtain

$$\frac{g}{2}cos{\theta }_2\left(\right|BC|h_{BC}^2-\left|DE\right|h_{DE}^2)+\frac{g}{2}cos{\theta }_3\left(\right|EF|h_{EF}^2-\left|GH\right|h_{GH}^2)-w\left[h_1{v}_1^2+g\frac{{h_1}^2}{2}\right]$$

$$+w\left[(h_2{u}_2^2+g\frac{h_2^2}{2})sin{\theta }_2+h_2{u}_2{v}_{2}cos{\theta }_2\right]+w\left[(h_3{u}_3^2+g\frac{h_3^2}{2})sin{\theta }_3+h_3{u}_3{v}_{3}cos{\theta }_3\right]=0.$$

(30)

Using the branch velocities and neglecting the transverse components we get

$$\frac{g}{2}cos{\theta }_2(|BC|h_{BC}^2-\left|DE\right|h_{DE}^2)+\frac{g}{2}cos{\theta }_3(|EF|h_{EF}^2-\left|GH\right|h_{GH}^2)$$

$$-w\left[h_1{v}_1^2+g\frac{{h_1}^2}{2}\right]+wsin{\theta }_2\left[h_2{u_2^{\parallel }}^2+g\frac{h_2^2}{2}\right]+wsin{\theta }_3\left[h_3{u_3^{\parallel }}^2+g\frac{h_3^2}{2}\right]=0.$$

(31)

4.4. Momentum Flux for the T-Fork

Consider now the T-geometry shown in the right panel of Figure 1. The calculations are simpler so that we used this geometry to validate the approach numerically. The general fork domain $\mathcal{F}$ can be reduced to the square $ADFIA$ by taking ${\theta }_2=\pi ,{\theta }_3=0$ and $B\to C\to A,G\to H\to I$. Then the Equations (28) and (30) reduce to

$$\begin{array}{c}\hfill -h_1{u}_1{v}_1-(h_2{u}_2^2+g\frac{h_2^2}{2})+h_3{u}_3^2+g\frac{h_3^2}{2}=0,\end{array}$$

(32)

$$\begin{array}{c}\hfill -(h_1{v}_1^2+g\frac{h_1^2}{2})-h_2{u}_2{v}_2+g\frac{h_{23}^2}{2}+h_3{u}_3{v}_3=0,\end{array}$$

(33)

where the term $h_{23}$ is

$$h_{23}^2\equiv \frac{1}{w}{\int }_{DF}{h}^{2}ds.$$

(34)

We will see that it can be obtained by interpolation of $h_2$ and $h_3$.

4.5. Effective 1D Model for the T-Fork

The pseudo-conservation laws (26), (27), (29) and (31) established in the previous section in the limit $w\to 0$ provide a formal connection between $h,u^{\parallel }$ in branches 1,2 and 3. In principle, they enable to approximate the 2D problem (9)–(11) by three 1D shallow water equations

$$\begin{array}{c}\hfill H_t^i+{\left(H^i{U}^i\right)}_x=0,\end{array}$$

(35)

$$\begin{array}{c}\hfill {(H^i{U}_i)}_t+{(H^i{U^i}^2+\frac{g{H}_i^2}{2})}_x=0,\end{array}$$

(36)

where $i=1,2,3$ correspond to the different branches. These 1D shallow water equations can be solved using a standard finite difference scheme, see for example [14]. The discretization is shown in Figure 3 where the first nodes in each branch have values $H=h_i,U=u_i$. The coupling equations between these three nodes given by (26), (27) and (29) would be solved using a Newton iteration.

5. Numerical Solutions of the 2D Shallow Water Equations

The approximation described in the previous section holds if the error remains small. We now evaluate this error by solving numerically the 2D problem (9)–(11), compute $h,u^{\parallel }$ and see how these values agree with the pseudo-conservation laws (26), (27), (29) and (31). We chose the T geometry shown in the right panel of Figure 1 for simplicity and considered symmetric and non symmetric initial conditions. We also increased the wave amplitude to estimate the effect of the non linearity.

The Equations (9)–(11) were discretized using as space unit the depth d. The time unit was $\sqrt{\frac{d}{g}}$. The variables and fields was rescaled as

$$x^{\prime }=\frac{x}{d},t^{\prime }=t\sqrt{\frac{g}{d}},h^{\prime }=\frac{x}{d},u^{\prime }=\frac{u}{\sqrt{gd}}.$$

(37)

This amounts to taking $d=1,g=1$ in (9)–(11).

We solved the nonlinear shallow water equations using a first order finite volume scheme on an unstructured triangular mesh produced with the Gmsh meshing software (see details in [15]). We used the width $w=0.125$ and the typical size of the triangles is $0.02$. The time advance used a variable order Adams–Bashforth–Moulton multistep solver (implemented in Matlab under ode113 subroutine [16]). The relative and absolute tolerances were set to ${10}^{-5}$.

The initial condition is taken as a travelling solitary wave of velocity c. This is an exact solution for the mass conservation law. We used a solitary wave inspired by the Serre theory [17], (see [18] for the modern variational derivation)

$$\begin{array}{c}\hfill h(x,y,t=0)=d+\eta \left(y\right),\end{array}$$

(38)

$$\begin{array}{c}\hfill v(x,y,t=0)=c\frac{\eta \left(y\right)}{d+\eta \left(y\right)},\end{array}$$

(39)

$$\begin{array}{c}\hfill \eta \left(y\right)=a{\mathrm{sech}}^2(\frac{1}{2}k(y-y_0)),\end{array}$$

(40)

where the speed is

$$c=\sqrt{g(d+a)}.$$

The other parameters were

$$g=1,k=1,d=1,a=1,x_0=y_0=2.5.$$

The wave was chosen so that its extension $2/k=2$ is much larger than the width $w=0.125$. Below we discuss the effect of the width.

The four pseudo-conservation laws for the mass, momenta and energy (26), (27) and (29) on the fork domain $ADFIA$ are

$$\begin{array}{c}\hfill \delta m\equiv -h_1{v}_1-h_2{u}_2+h_3{u}_3=0,\end{array}$$

(41)

$$\begin{array}{c}\hfill \delta p_x\equiv -h_1{u}_1{v}_1-(h_2{u}_2^2+g\frac{h_2^2}{2})+h_3{u}_3^2+g\frac{h_3^2}{2}=0.\end{array}$$

(42)

$$\begin{array}{c}\hfill \delta p_y\equiv -(h_1{v}_1^2+g\frac{h_1^2}{2})-h_2{u}_2{v}_2+g\frac{h_{23}^2}{2}+h_3{u}_3{v}_3=0,\end{array}$$

(43)

$$\begin{array}{c}\hfill \delta e\equiv -v_1(g{h}_1^2+h_1\frac{v_1^2}{2})-u_2(g{h}_2^2+h_2\frac{u_2^2}{2})+u_3(g{h}_3^2+h_3\frac{u_3^2}{2})=0,\end{array}$$

(44)

where we introduced the residuals $\delta m,\delta p_x,\delta p_y$ and $\delta e$.

We considered a symmetric situation where the wave is incident from branch 1 and a non symmetric situation where the wave was send into the fork from branch 3. In both cases, the number of unknowns was the same; see

Table 1. The two different dynamic problems for the T-branch.

Type	Known	Unknown
wave in branch 1	$h_1,v_1$	$h_2,u_2,h_3,u_3$
wave in branch 3	$h_3,u_3$	$h_1,v_1,h_2,u_2$

.

The wave mass and wave energy in each branch have been calculated. They are defined as

$$M_w={\int }_{\Omega }(h-d)dxdy,$$

$$E_w={\int }_{\Omega }\frac{1}{2}\left[g{(h-d)}^2+(u^2+v^2)h\right]dxdy.$$

Energy will propagate very differently in problems 1 and 2. In the next sections we examine in detail the two types of problems and use the conservation laws to establish jump conditions for the 1D effective model.

To verify the approximation given by the relations (41)–(44), we also computed the time evolution of the quantities $h_1,h_2,h_3,v_1,u_2,u_3$ from the 2D direct numerical simulations. We used a scattered linear interpolation to estimate these physical variables along the four different segments of the fork region from the unstructured triangular mesh data.

5.1. Wave Incident into Branch 1

5.1.1. Small Amplitude Waves $a/d=0.1$

The time evolution of the wave mass and energy is presented in Figure 4. Consider the wave mass, at $t=0$: $M_1^0=57{10}^{-3},M_2^0=M_3^0=0$. After the wave has passed, at $t=6.5,M_1=0,M_2=M_3=26$. We have $2\times 26=52$ which shows the conservation of mass. Notice the depression in the mass in branch 1 after the wave passes. Almost all energy is transferred to branches 2 and 3.

Here our balance laws hold well for both the mass and the energy, see Figure 5.

We can use them to obtain $u_2,h_2$. Assume symmetry $h_2=h_3,u_2=-u_3$. The balance laws reduce to

$$\begin{array}{c}\hfill -h_1{v}_1-2h_2{u}_2=0,\end{array}$$

(45)

$$\begin{array}{c}\hfill -v_1(g{h}_1^2+h_1{v}_1^2/2)-2u_2(g{h}_2^2+h_2{u}_2^2/2)=0.\end{array}$$

(46)

Since $v_1^2,u_2^2\ll g{h}^2$ we can neglect the terms $v_1^2,u_2^2$ of the second equation. The resulting relations are satisfied by

$$h_2=h_1,u_2=-v_1/2,$$

(47)

which are the Stoker conditions. These are in good agreement with the simulations as shown by Figure 6.

5.1.2. Very Large Amplitude Waves $a/d=2$

In this case, 2D effects start to appear. Figure 7 shows a snapshot of the surface elevation h for a wave such that $a/d=2$. Notice the lump $h\approx 2$ on the edge of the domain.

Figure 8 shows the time evolution of the wave mass and energy. Despite the evidence of 2D effects, the overall transfer of wave mass and wave energy from branch 1 to branches 2 and 3 does not vary significantly as $a/d$ changes from $0.1$ to 2.

Figure 9 shows the time evolution of $\delta m$ and $\delta e$. Notice that the mass relation is better satisfied than the energy relation.

Again the Stoker relations (47) give a good approximation as shown by Figure 10 which show that $h_2\approx h_1$ and $u_2\approx v_1/2$. The price to pay to approximate the 2D situation by a 1D effective model is an energy loss at the junction.

Also remark that for the approximation to hold it is crucial that the wave be wider than w and not too fast. If these conditions are not met, $h_2$ and $u_2$ will be delayed from $h_1,v_1$ and will need to describe what happens in the fork. We observed this for a larger channel $w=1$ and the same parameters.

5.2. Wave Incident into Branch 3

For this configuration, we observe a significant difference in behavior as the wave amplitude increases. Figure 11 shows the time evolution of the wave mass and wave energy for $a/d=0.1$ (top panels) and $a/d=2$ (bottom panels). Small amplitude waves get transmitted to branch 1 as much as to branch 2. On the other hand, large amplitude waves are predominantly transmitted to branch 2. The mass entering branch 2 is three times larger than the one entering branch 1; for energies, the factor is six.

5.2.1. Small Amplitude Waves $a/d=0.1$

First observe that $u_1$ is non zero and close to $v_1$. Nevertheless, the mass and energy residuals $\delta m$ and $\delta e$ are small as seen in Figure 12. The wave elevation h does not vary much from one branch to the other as seen in the top panel of Figure 13. The velocities $u_2$ and $v_1$ verify $u_2\approx u_3/2,v_1\approx u_3/2$. These two results show that the Stoker conditions hold for this small amplitude.

5.2.2. Large Amplitude Waves $a/d=1$

Figure 14 shows $h(t=0.8)$ for a wave incident in branch 3 for $a/d=2$. Notice the complex structure of the flow at the junction. There is some recirculation so that the flow is essentially 2D and not amenable to a 1D reduction.

Nevertheless, for a smaller amplitude $a/d=1$, the balance laws (41)–(44) give some insight into the flow. Figure 15 shows the mass $\delta m$ and energy $\delta e$. The mass is much better conserved than the energy.

The momenta (42) and (43) are plotted in Figure 16.

When the wave is coming from branch 3, an obvious solution is

$$v_1=0,u_2=u_3,h_2=h_3,h_1=h_2.$$

(48)

This is simplistic, in reality $v_1\ne 0$ but remains small. The horizontal component $u_1$ is non zero and close to $u_2$ as shown in Figure 17.

The mass equation and y momentum equations allow to extract relations between $v_1,h_1,h_2,u_2,v_1,h_3,u_3$. Assuming $v_1,u_2,u_3$ smaller than $h_1^2,h_2^2,h_3^2$, we have

$$\begin{array}{c}\hfill v_1=\frac{h_3{u}_3-h_2{u}_2}{h_1},\end{array}$$

(49)

$$\begin{array}{c}\hfill h_1=h_{23}.\end{array}$$

(50)

The quantity (34) in the y component of the momentum is computed from the numerical solution. It is plotted as a function of time together with the estimate

$$h_{23}^i=\sqrt{\frac{1}{2}(h_2^2+h_3^2)},$$

(51)

in the left panel of Figure 18. As can be seen, the agreement is very good.

The velocity $v_{1m}$ given by the mass conservation relation agrees semi-quantitatively with the value $v_1$ estimated from the 2D numerical solution. Both quantities are plotted as a function of time in Figure 19. Note the delay due to the time the wave needs to propagate from one interface to the other. The y momentum conservation law is not satisfied so that there is no additional equation to estimate $u_2$.

6. Discussion and Conclusions

The results of the previous section show that for large amplitudes and an asymmetric fork Stoker’s interface conditions do not hold and the angle of the fork plays a role. This seems to contradict the findings of Shi et al. [9]. Two reasons show that there is no contradiction. First, the amplitude of our waves ($a/d\approx 1$) are much larger than the ones presented in [9] ($a/d\approx 0.3$) so that nonlinear effects are much stronger in our study. The other point is that the ${\mathrm{sech}}^2$ initial condition is an exact solution of the Boussinesq equations, but not of the nonlinear shallow water equations. For the Boussinesq equations, we also expect an angle dependence, even for narrow channels, when the amplitude becomes large. To see this, we examine the reduction of the equations for a fork.

The Boussinesq equations read

$$\begin{array}{c}\hfill h_t+\nabla ·[(1+h)\nabla \phi ]=0,\end{array}$$

(52)

$$\begin{array}{c}\hfill {\phi }_t+\frac{1}{2}{(\nabla \phi )}^2+h-\frac{1}{3}{(\Delta \phi )}_t=0,\end{array}$$

(53)

where $h(x,y,t)$ is the water elevation. The velocity potential $\phi (x,y,t)$ is such that ${(u,v)}^T=\nabla \phi$. The boundary conditions are non slip $\nabla \phi ·\mathbf{n}=0$. Integrating the equations on the fork domain $\mathcal{F}$ (left panel of Figure 1) we get

$$\begin{array}{c}\hfill {\partial }_t{\int }_{\mathcal{F}}hdxdy-{\int }_{IA\cup CD\cup FG}(1+h)\nabla \phi ·\mathbf{n}ds,\end{array}$$

(54)

$$\begin{array}{c}\hfill {\partial }_t{\int }_{\mathcal{F}}(\phi -\frac{1}{3}\Delta \phi )dxdy+{\int }_{\mathcal{F}}(\frac{1}{2}{(\nabla \phi )}^2+h)dxdy=0.\end{array}$$

(55)

Neglecting the time evolution in the fork region, we get the following interface conditions

$$\begin{array}{c}\hfill (1+h_1)u_1^{\parallel }+(1+h_2)u_2^{\parallel }+(1+h_3)u_3^{\parallel }=0,\end{array}$$

(56)

$$\begin{array}{c}\hfill {\int }_{\mathcal{F}}\frac{1}{2}\left[{(\nabla \phi )}^2+h\right]dxdy=0.\end{array}$$

(57)

Note how the first equation reduces to Kirchhoff’s law for small h. The second equation contains is an integral over the whole domain and depends on the angle of the fork. For small angles, we can assume that $\nabla \phi =u^{\parallel }$ so that the conditions reduce to

$$\begin{array}{c}\hfill (1+h_1)u_1^{\parallel }+(1+h_2)u_2^{\parallel }+(1+h_3)u_3^{\parallel }=0,\end{array}$$

(58)

$$\begin{array}{c}\hfill \frac{1}{2}{(u_1^{\parallel })}^2+h_1+\frac{1}{2}{(u_2^{\parallel })}^2+h_2+\frac{1}{2}{(u_3^{\parallel })}^2+h_3=0.\end{array}$$

(59)

Not surprisingly, these conditions are very close to the ones obtained by Nachbin and Simoes [10], except for the Jacobian of the conformal transformation.

To conclude, we studied the propagation of shallow water waves in a fork between three narrow channels. We considered both the 2D numerical solution and a homothetic reduction procedure that gives coupling conditions at the interface. For such narrow widths, the delay experienced by the wave is negligible so that one can envision describing the junction by an effective 1D PDE model.

Our reduction enabled us to derive balance laws for the mass, momenta and energy of the flow across a general junction. For small amplitude waves, these laws reduce to the commonly used Stoker jump conditions, giving these a formal justification. We verified these Stoker conditions on the 2D numerical solutions of the shallow water equations for symmetric and non symmetric conditions. Then, the angle of the junction does not play any role. This happens also for a general nonlinear wave equation; we had seen this a previous study for the particular case of the sine-Gordon equation [2].

For large amplitude shallow water waves, the situation depends on the symmetry of the fork. For a symmetric fork, the Stoker conditions are approximately verified. This is explained by the strong constraint imposed by the symmetry. Then, the only solution of the balance laws corresponds to the Stoker conditions. When the fork is non symmetric as in our case 2, more information is needed about what happens inside the fork. The quantities $u_i^{\parallel },i=1,2,3$ are velocities projected in the direction of the branches and this projection leads to a loss of information. Far from the junction, the flow is quasi-1D so that not much is lost. On the contrary, inside the junction, the flow is full 2D. A possible solution, to be studied in the future would be to use the full conservation law including the time dependent term. Then we would introduce a fictitious node inside the junction and couple it to the boundaries using average differential equations obtained by integrating (9)–(11) on the fork domain.