Preconditioned augmented Lagrangian method for mean curvature image deblurring

1.   Introduction

The dynamics of network flows is usually modelled by systems of Partial Differential Equations (briefly PDEs), most of time balance laws. The dynamics is defined on a topological graph with evolution on arcs given by system of PDEs, while additional conditions must be assigned at network nodes, e.g. conservation of mass and momentum. There is a large literature devoted to these problems and we refer to [5] for a extensive survey and for additional references.

In particular, here we focus on water flows on a oriented network of open canals and the model given by Saint-Venant or shallow water equations. The latter form a non linear system of balance laws composed by a mass and momentum balance laws. In water management problems, these equations are often used as a fundamental tool to describe the dynamics of canals and rivers, see [1] and papers in same volume, and various control techniques were proposed, see [2,3,17,15,19,23,26] and references therein. Moreover, the need of dynamic models in water management is well documented, see [25]. The shallow water system is hyperbolic (except when water mass vanishes) and has two genuinely nonlinear characteristic fields. Moreover, it exhibits two regimes: fluvial or sub-critical, when one eigenvalue is negative and one positive, and torrential or super-critical, when both eigenvalues are positive. This is captured by the so called Froude number, see (1). For a complete description of the physics of the problem one needs to supply the equations with conditions at nodes, which represent junctions. The junction conditions are originally derived by engineers in the modeling of the dynamic of canals and rivers. The first and most natural condition is the conservation of water mass which is expressed as the equality between the sum of fluxes from the incoming canals and that from outgoing ones. One single condition is not sufficient to isolate a unique solution, thus different additional condition were proposed in the literature. Physical reasons motivate different choices of conditions, among which the equality of water levels, of energy levels and conservation of energy. For the assessment of coupling conditions on canals networks and for more details on the existence of solutions in the case of subcritical flows, we refer the reader to [9,18,16,20,14,22,24]. For discussion on supercritical flow regimes, we refer the reader to [18] and references there in.

Then, to construct solutions one may resort to the concept of Riemann solver at a junction, see [13]. A Riemann solver at a junction is a map assigning solutions to initial data which are constant on each arc. Alternatively one may assign boundary conditions on each arc, but, due to the nonlinearity of equations, one has to make sure that boundary values are attained. This amounts to look for solutions with waves having negative speed on incoming channels and positive on outgoing ones: in other words waves do not enter the junction. A Riemann solver with such characteristics is called consistent, see also [12].

In this paper we are interested in transitions between different flow regimes, when the transition occurs at a junction of a canals network. We assume to have incoming canals which end at the junction and outgoing canals which start at the junction. Thus we formulate a left-half Riemann problem for incoming canals and a right-half Riemann problem for outgoing canals to define the region of admissible states such that waves do not propagate into the junction. This corresponds to identify the regions where Riemann solvers can take values in order to be consistent. Such regions are enclosed by the Lax curves (and inverted Lax curves) and the regime change curves. To help the geometric intuition, we developed pictures showing such curves and the regions they enclose.

The definitions described above and given in Section 4, are the necessary basis for an analysis on a complex network. Due to the complexity of the problem, we consider as case study the specific case of two identical canals interconnected at a junction (simple junction). We start focusing on conservation of water through the junction and equal height as coupling conditions. It is typically expected the downstream flow to be more regular, thus we consider three cases: fluvial to fluvial, torrential to fluvial and torrential to torrential. In the fluvial to fluvial case there exists a unique solution. However such solution may be different than the solution to the same Riemann problem inside a canal (without the junction) and may exhibit the appearance of a torrential regime. The torrential to fluvial case is more delicate to examine. Three different cases may happen: the solution propagates the fluvial regime upstream, the solution propagates the torrential regime downstream or no solution exits. Finally, in the torrential to torrential case, if the solution exists then it is torrential.

To illustrate the achieved results we perform simulations using a Runge-Kutta Discontinuous Galerkin scheme [6]. The RKDG method is an efficient, effective and compact numerical approach for simulations of water flow in open canals. Specifically, it is a high-order scheme and compact in the sense that the solution on one computational cell depends only on direct neighboring cells via numerical fluxes, thus allowing for easy handling the numerical boundary condition at junctions. In the first example we show a simulation where an upstream torrential regime is formed starting from special fluvial to fluvial conditions. The second example shows how a torrential regime may propagate downstream.

We conclude by discussing the possible solutions if the water height condition is replaced by the equal energy condition.

The paper is organized as follows: in Section 2, we present the model starting from the one-dimensional shallow water equations. In Section 3, we give useful notations and preliminary results that allow to determine the admissible states for the half-Riemann problems discussed in the following Section 4. In Section 5 we study possible solutions at a simple junction for different flow regimes and different junction conditions. Finally, in Section 6 we illustrate the results of the previous section with a couple of numerical tests.

2.   Flow classification and governing equations

The most common and interesting method of classifying open-channel flows is by dimensionless Froude number, which for a rectangular or very wide channel is given by the formula:

where $v$ is the average velocity and $h$ is the water depth. The three flow regimes are:

● $Fr < 1$ subcritical flow or fluvial regime;

● $Fr = 1$ critical flow;

● $Fr>1$ supercritical flow or torrential regime.

The Froude-number denominator $(gh)^{1/2}$ is the speed of an infinitesimal shallow-water surface wave. As in gas dynamics, a channel flow can accelerate from subcritical to critical to supercritical flow and then return to subcritical flow through a shock called a hydraulic jump, see [11] and references there in.

We are interested in the transition between different flow regimes when it occurs at a junction of a canals network. On each canal the dynamics of water flow is described by the following system of one-dimensional shallow water equations

The quantity $q = hv$ is often called discharge in shallow water theory, since it measures the rate of water past a point. We write the system as:

where

For smooth solutions, these equations can be rewritten in quasi-linear form

where the Jacobian matrix $f^\prime(u)$ is given by

The eigenvalues of $f^\prime(u)$ are

with corresponding eigenvectors

The shallow water equations are strictly hyperbolic away from $h = 0$ and both 1-th and 2-th fields are genuinely nonlinear ($\nabla\lambda_j(u)\cdot r_j(u)\ne 0$, $j = 1, 2$). Note that $\lambda_1$ and $\lambda_2$ can be of either sign, depending on the magnitude of $v$ relative to $\sqrt{gh}$, so depending on the Froude number.

Solutions to systems of conservation laws are usually constructed via Glimm scheme of wave-front tracking [10,21]. The latter is based on the solution to Rieman problems:

Here $u(x, 0) = (h(x, 0), q(x, 0))$ and $u_l = (h_l, q_l)$ and $u_r = (h_r, q_r)$. The solution always consists of two waves, each of which is a shock or rarefaction, thus we first describe these waves.

(R) Centered Rarefaction Waves. Assume $u^+$ lies on the positive $i$-rarefaction curve through $u^-$, then we get

where, for the 1-family

for $\xi\in[v^+-\sqrt{h^-}, v^-+2\sqrt{h^-})$, and for the second family

for $\xi\in[\lambda_2(u^-), \infty)$.

(S) Shocks. Assume that the state $u^+$ is connected to the right of $u^-$ by an $i$-shock, then calling $\lambda = \lambda_i(u^+, u^-)$ the Rankine-Hugoniot speed of the shock, the function

provides a piecewise constant solution to the Riemann problem. For strictly hyperbolic systems, where the eigenvalues are distinct, we have that

3.   The geometry of Lax and regime change curves

To determine a solution for problems on a network, we need to analyze in detail the shape of shocks and rarefaction curves and, more generally, of Lax curves (which are formed by joining shocks and rarefaction ones, see [4]). We start fixing notations and illustrating the shapes of curves.

For a given point $(h_0, v_0)$, we use the following notations for shocks and rarefaction curves:

Moreover, we define the inverse curves:

Similarly, we set:

We will also consider the regime transition curves: the 1-critical curve is given by

and the 2-critical curve by

In Figure 1 we illustrate the shape of these curves.

To construct a solution to a Riemann problem $(u_l, u_r)$, we define the Lax curves. Given a left state $u_l = (h_l, q_l)$ the Lax curve is given by:

For the right state $u_r = (h_r, q_r)$, we define the inverse Lax curve:

Remark 1. The Riemann problem for shallow water equations 2 with left state $u_l$ and right state $u_r$ has a unique solution if and only if the Lax curves $\phi_l(h)$ and $\phi_r(h)$ have a unique intersection. In that case, the intersection will be called the middle state $u_m$. As shown in Figure 2, the function $v = \phi_l(h)$ is strictly decreasing, unbounded and starting at the point $v_l+2\sqrt{gh_l}$ and, $v = \phi_r(h)$ is strictly increasing, unbounded, with minimum $v_r-2\sqrt{gh_r}$. Thus, the Riemann problem for shallow water has a unique solution in the region where

When working with $(h, q)$ variables we use the following notations for Lax curves and regime transition curves:

Moreover, for a given value $(h_i, v_i)$ (or $(h_i, q_i)$) we set

3.0.1.   The Lax curve $\tilde{\phi}_l(h)$

In this subsection we study in detail the properties of the function $q = \tilde{\phi}_l(h)$. For a given left state $u_l = (h_l, q_l)$,

with

By computing its first and second derivatives,

and

we can conclude that $\tilde{\phi}_l\in C^2(]0, +\infty[)$. Specifically, in $]0, +\infty[$ $\tilde{\phi}_l^{\prime\prime}<0$, $\tilde{\phi}_l^\prime(h)$ is strictly decreasing and $\tilde{\phi}_l$ is a strictly concave function. As

we investigate two different cases:

Case 1. If $v_l\leq -2\sqrt{g h_l}$, $\tilde{\phi}_l$ is a strictly negative decreasing concave function and specifically $\tilde{\phi}_l(h) < \tilde{\mathcal{C}}^-(h)$ for $h>0$. Therefore, in this case $\tilde{\phi}_l$ never intersects the critical curves $\tilde{\mathcal{C}}^+$ and $\tilde{\mathcal{C}}^-$.

Case 2. If $v_l>-2\sqrt{g h_l}$, the function $\tilde{\phi}_l$ admits a maximum point $h_l^{max}$, so it is increasing in $(0, h_l^{max})$ an decreasing in $(h^{max}_l, +\infty)$. In this case the function $\tilde{\phi}_l$ intersects the two critical curves $\tilde{\mathcal{C}}^+$ and $\tilde{\mathcal{C}}^-$. The intersection points vary with the choice of the left state $u_l$, see Figure 3. To compute this points, we distinguish the following subcases: $u_l$ is such that $-2\leq\mathcal{F}_l\leq 1$ or such that $\mathcal{F}_l>1$.

Case 2.1. For $-2<\mathcal{F}_l \leq 1$, it is the rarefaction portion of $\tilde{\phi}_l$ to intersect the 1-critical curve $\tilde{\mathcal{C}}^+$ at $u^+_{l, \mathcal{R}} = (h^+_{l, \mathcal{R}}, \tilde{\mathcal{C}}^+(h^+_{l, \mathcal{R}}))$ with

and we have that the maximum point $h_l^{max}$ is such that $\tilde{\phi}_l(h_l^{max})\equiv \tilde{\mathcal{C}}^+(h_l^{max})$, i.e. $h_l^{max}\equiv h^+_{l, \mathcal{R}}$. In the special case $\mathcal{F}_l = 1$ we get $h_l^{max}\equiv h_l$. Moreover,

Notice that for $-2<\mathcal{F}_l\leq -1$ the rarefaction portion of $\tilde{\phi}_l$ intersects $\tilde{\mathcal{C}}^-$ at $u^-_{l, \mathcal{R}} = (h^-_{l, \mathcal{R}}, \tilde{\mathcal{C}}^-(h^-_{l, \mathcal{R}}))$ with

while for $-1<\mathcal{F}_l \leq 1$, the shock portion of $\tilde{\phi}_l$ intersects $\tilde{\mathcal{C}}^-$ at $u^-_{l, \mathcal{S}} = (h^-_{l, \mathcal{S}}, \tilde{\mathcal{C}}^-(h^-_{l, \mathcal{S}}))$ with $h^-_{l, \mathcal{S}}$ given by the following condition:

Case 2.2. For $\mathcal{F}_l >1$, the maximum value of $\tilde{\phi}_l$ is reached by the shock portion and so $h_l^{max}>h_l$. The shock portion intercepts the critical curve $\tilde{\mathcal{C}}^+$ at $u^+_{l, \mathcal{S}} = (h^+_{l, \mathcal{S}}, \tilde{\mathcal{C}}^+(h^+_{l, \mathcal{S}}))$ with

and the curve $\tilde{\mathcal{C}}^-$ at $u^-_{l, \mathcal{S}}$ defined in 21.

Notice that for $h\geq h_l$ the equation

has two solutions:

Moreover, for $\mathcal{F}_l>1$, we have $h_l^*>h^+_l$, where the height $h^+_l = (q_l^2/g)^\frac 13$ is given by the intersection between $\tilde{\mathcal{C}}^+$ and the horizontal line $q = q_l$. That is:

Indeed, using the relation $q_l^2 = gh_l^3\mathcal{F}_l^2$, if $\mathcal{F}_l>1$ then $2\mathcal{F}_l^{\frac{4}{3}} - \mathcal{F}_l^{\frac{2}{3}}-1 > 0.$ So, for $(h_l, q_l)$ supercritical, the point $(h^*_l, q_l)$ is subcritical, i.e. $|\mathcal{F}^*_l|<1$.

3.0.2.   The Lax curve $\tilde{\phi}_r(h)$

Here we study the properties of the function $\tilde{\phi}_r(h)$, given by:

By straightforward computations we get its derivatives:

Then, $\tilde{\phi}_r(h)\in C^2(]0, +\infty[)$ and in $]0, +\infty[$ $\tilde{\phi}_l^{\prime\prime}>0$, $\tilde{\phi}_l^\prime(h)$ is strictly increasing and $\tilde{\phi}_l$ is a strictly convex function. As

we investigate two different cases:

Case 1. If $v_r\geq 2\sqrt{gh_r}$, $\tilde{\phi}_r$ is a strictly positive increasing convex function and specifically $\tilde{\phi}_r(h) > \tilde{\mathcal{C}}^+(h)$, for $h>0$. Therefore, $\tilde{\phi}_r$ never intersects the critical curves $\tilde{\mathcal{C}}^+$ and $\tilde{\mathcal{C}}^-$.

Case 2. If $v_r< 2\sqrt{g h_l}$, the function $\tilde{\phi}_r$ admits a minimum point $h_r^{min}$. In this case the function $\tilde{\phi}_r$ intersects the two critical curves $\tilde{\mathcal{C}}^+$ and $\tilde{\mathcal{C}}^-$. As done before we distinguish two subcases:

Case 2.1. For $-1<\mathcal{F}_r \leq 2$, the rarefaction portion of $\tilde{\phi}_r$ intersects the 2-critical curve $\tilde{\mathcal{C}}^-$ at $u^-_{r, \mathcal{R}} = (h^-_{r, \mathcal{R}}, \tilde{\mathcal{C}}^-(h^-_{r, \mathcal{R}}))$, with

and we have that the minimum point $h_r^{min}$ is such that $\tilde{\phi}_r(h_r^{min})\equiv \tilde{\mathcal{C}}^-(h_r^{min})$, i.e. $h_r^{min}\equiv h^-_{r, \mathcal{R}}$. In the special case $\mathcal{F}_r = 1$ we get $h_r^{min}\equiv h_r$. Moreover,

For $1<\mathcal{F}_r\leq 2$ the rarefaction portion of $\tilde{\phi}_r$ intersects $\tilde{\mathcal{C}}^+$ at $u^+_{r, \mathcal{R}} = (h^+_{r, \mathcal{R}}, \tilde{\mathcal{C}}^+(h^+_{r, \mathcal{R}}))$, with

while for $-1<\mathcal{F}_r \leq 1$, the shock portion of $\tilde{\phi}_r$ intersects $\tilde{\mathcal{C}}^+$ at $u^+_{r, \mathcal{S}} = (h^+_{r, \mathcal{S}}, \tilde{\mathcal{C}}^+(h^+_{r, \mathcal{S}}))$, with $h^+_{r, \mathcal{S}}$ such that

Case 2.2. For $\mathcal{F}_r <-1$, the minimum value of $\tilde{\phi}_r$ is reached by the shock portion and so $h_r^{min}>h_r$. The shock portion intercepts the critical curves $\tilde{\mathcal{C}}^-$ at $u^-_{r, \mathcal{S}} = (h^-_{r, \mathcal{S}}, \tilde{\mathcal{C}}^-(h^-_{r, \mathcal{S}}))$, with

and curve $\tilde{\mathcal{C}}^+$ at $u^+_{r, \mathcal{S}}$ given by 26.

Notice that the equation

has two solutions:

Moreover, for $\mathcal{F}_r<-1$ we have $h_r^*>h^-_r = (q_r^2/g)^\frac 13$, where the height $h_r^-$ is given by the intersection between $\tilde{\mathcal{C}}^-$ and the horizontal line $q = q_r$. So, for $(h_r, q_r)$ supercritical, the point $(h^*_r, q_r)$ is subcritical.

4.   The half Riemann problems

4.1.   Left-half Riemann problem

[The case of an incoming canal] We fix a left state and we look for the right states attainable by waves of non-positive speed.

Fix $u_l = (h_l, q_l)$, we look for the set $\mathcal{N}(u_l)$ of points $\hat u = (\hat h, \hat q)$ such that the solution to the Riemann problem

contains only waves with non-positive speed. We distinguish three cases:

● Case A: the left state $u_l$ is such that $|\tilde{\mathcal{F}}_l| < 1$;

● Case B: the left state $u_l$ is such that $\tilde{\mathcal{F}}_l > 1$;

● Case C: the left state $u_l$ is such that $\tilde{\mathcal{F}}_l <- 1$.

4.1.1.   Case A

For this case we refer to Figure 4.We identify the set $\mathcal{N}^A(u_l)$ as the union of three regions $\mathcal{I}^A_1$, $\mathcal{I}^A_2$ and $\mathcal{I}^A_3$ defined in the plane $(h, q)$. The first region is identified by all points that belong to the curve $\tilde{\phi}_l(h)$ such that $\tilde{\mathcal{C}}^-(h)\leq\tilde{\phi}_l(h)\leq \tilde{\mathcal{C}}^+(h)$, i.e.

where the points $h^+_{l, \mathcal{R}}$ and $h^-_{l, \mathcal{S}}$ are given in 19 and 21 respectively. The second region is defined as follows

The last region $\mathcal{I}^A_3$ is defined by the set of all possible right states $\hat u$ that can be connected by a 2-shock with non-positive speed to an intermediate state $u_m$ lying on $\tilde{\phi}_l(h)$ curve such that $|\tilde{\mathcal{F}}_m|\leq 1$ and

To define this region, we have to look for values $q = hS_2(h_m, v_m;h)$ as $h<h_m$ such that $q\geq q_m$. That is,

This inequality is verified for $-1\leq \mathcal{F}_m <0$ and for all $h\leq h^*_m$ with $h^*_m$ given by

We obtain (see Figure 4),

4.1.2.   Case B

For this case we refer to Figure 5. It is always possible to connect the left value $u_l$ to a value $u^*_l$ by a 1-shock with zero speed (vertical shock). Specifically, we set $h^*_l>h_l$ such that $q^*_l = q_l$, i.e.

as previously computed in 23. Moreover, as previously observed at the end of subsection 3.0.1 the value $(h^*_l, q_l)$ falls in the subcritical region, then

where $\mathcal{N}^A(u_l)$ falls under the previous Case A.

4.1.3.   Case C

For this case we refer to Figure 6. We have that: if $\mathcal{F}_l \leq -2$ then

otherwise if $-2<\mathcal{F}_l<1$ the intersection point $u^-_{l, \mathcal{R}}$ between $h\mathcal{R}_1(h_l, v_l;h)$ and $\tilde{\mathcal{C}}_-(h)$ given in 20 is different from zero and defines the admissible region

4.2.   Right-half Riemann problem

[The case of an outgoing canal] We fix a right state and we look for the left states attainable by waves of non-negative speed. For sake of space the figures illustrating these cases will be postponed to the Appendix.

Fix $u_r = (h_r, q_r)$, we look for the set $\mathcal{P}(u_r)$ of points $\tilde u = (\tilde h, \tilde q)$ such that the solution to the Riemann problem

contains only waves with non-negative speed. As in the previous case we identify three cases:

● Case A: the right value $u_r$ is such that $|\tilde{\mathcal{F}}_r|<1$.

● Case B: the right value $u_r$ is such that $\tilde{\mathcal{F}}_r>1$

● Case C: the right value $u_r$ is such that $\tilde{\mathcal{F}}_r<-1$.

4.2.1.   Case A

For this case we refer to Figure 12 in the Appendix. We identify the set $\mathcal{P}^A(u_r)$ as the union of three regions $\mathcal{O}^A_1$, $\mathcal{O}^A_2$ and $\mathcal{O}^A_3$ defined in the plane $(h, q)$. The first region is defined by all points that belong to the curve $\tilde{\phi}_r(h)$ such that $\tilde{\mathcal{C}}^-(h)\leq\tilde{\phi}_r(h)\leq\tilde{\mathcal{C}}^+(h)$, i.e.

where the points $h^-_{r, \mathcal{R}}$ and $h^+_{r, \mathcal{S}}$ are given in 24 and 26 respectively.

The second region is such that

The third region is defined by the set of all possible left states $\tilde u$ that can be connected by a 1-shock with non-negative speed to a middle state $u_m$ lying on $\tilde{\phi}_r(h)$ curve such that $|\tilde{\mathcal{F}}_m|\leq 1$ and

To define this region we have to look for values $q = h\mathcal{S}_1^{-1}(h_m, v_m;h)$ for $h<h_m$ such that $q_m\geq q$. Then, following the same reasoning done for the left-half Riemann problem we get

4.2.2.   Case B

For this case we refer to Figure 13 in the Appendix. If $\mathcal{F}\geq 2$ then

otherwise, if $1\leq \mathcal{F} <2$ the intersection point $u^+_{r, \mathcal{R}}$ between $\mathcal{R}_2^{-1}(u_r;h)$ and $\mathcal{C}^+(h)$, given in 24 is different from zero and defines the admissible region

4.2.3.   Case C

For this case we refer to Figure 14 in the Appendix. It is always possible to connect the right value $u_r$ to a value $u^*_r$ by a 2-shock with zero speed (vertical shock). Specifically, we set $h^*_r>h_r$ such that $q^*_r = q_r$, i.e.

as done in 27. Moreover, as previously observed at the end of subsection 3.0.2, the point $(h^*_r, q_r)$ is subcritical and then

where $\mathcal{P}^A(u^*_r)$ falls under the previous Case A.

5.   A simple junction

We consider here as case study a fictitious network formed by two canals intersecting at one single point, which artificially represents the junction. The junction is straight and separates two equal canals, one is the continuation of the other. This simple scenario appears to be like considering a problem for one straight canal, but by adding a fictitious junction we mimic a network and we provide the first analysis necessary for addressing more complicated networks for which the solution strongly depends on the given assumptions at the junction.

We name the canals such that 1 is the incoming canal and 2 is the outgoing ones. We indicate by $u^b = (h^b, v^bh^b) = (h^b, q^b)$ the traces at the junction. The flow is given by the one-dimensional shallow-water equations 2 in each canal coupled with special conditions at the junction. Our aim is to define and solve Riemann problems at the junction.

A Riemann Problem at a junction is a Cauchy Problem with initial data which are constant on each canal incident at the junction. So, assuming constant initial conditions $u^0_1$, $u^0_2$ on canal 1 and 2 respectively, the Riemann solution consists of intermediate states $u^b_1$ and $u^b_2$ such that $u^b_1 \in\mathcal{N}(u^0_1)$ and $u^b_2\in\mathcal{P}(u^0_2)$ and verifying given coupling conditions. Here, we are interested in evaluating possible solutions for different types of initial data belonging to different flow regimes. We start assuming, as junction conditions, the conservation of mass

and equal heights

In the following we study the boundary solution $u^b$ at the junction in the following interesting cases: A$\rightarrow$A, the water flow is fluvial in both incoming and outgoing canals around the junction; B$\rightarrow$A, the water flow is torrential with positive velocity in the incoming canal while it is fluvial in the outgoing one; B$\rightarrow$B, the water flow is torrential with positive velocity in both incoming and outgoing canals.

Case A$\rightarrow$A (Fluvial $\rightarrow$ Fluvial). Here we assume to have a left state $u_l$ and a right state $u_r$ such that $|\mathcal{F}_l|< 1$ and $|\mathcal{F}_r|< 1$. The solution at the junction consists of two waves separated by intermediate states $u^b_1$ and $u^b_2$ such that

with $h^b_1, h^b_2>0$.

Proposition 1. Under the subcritical condition on $u_l$ and $u_r$, the system 42 admits a unique solution.

Proof. We distinguish two cases:

Case 1. The two curves $\tilde{\phi}_l$ and $\tilde{\phi}_r$ intersect inside the subcritical region (Figure 7). In this case the solution is trivially the intersection point.

Case 2. The two curves $\tilde{\phi}_l$ and $\tilde{\phi}_r$ do not intersect inside the subcritical region. If $h_r<h_l$ (specifically $h^+_{r, \mathcal{S}}<h^+_{l, \mathcal{R}}$, see Figure 8) then the only point that verifies the junction conditions 42 is the critical point $u^+_{l, \mathcal{R}}$ given in 19, thus

If $h_r>h_l$ (specifically $h^-_{r, \mathcal{R}}>h^-_{l, \mathcal{S}}$, see Figure 9) the only point that verifies the junction conditions 42 is the critical point $u^-_{r, \mathcal{R}}$ given in 24, then

Remark 2. Notice that the proposed procedure may give a solution which is different from the classical solution of the Riemann problem on a single channel, given by the intersection point of $\phi_l$ and $\phi_r$ curves, see Remark 1.

Case B$\rightarrow$A (Torrential $\rightarrow$ Fluvial). Here we assume to have a left state $u_l$ and a right state $u_r$ with $\mathcal{F}_l>1$ and $|\mathcal{F}_r|<1$. The solution at the junction consists of two waves separated by intermediate states $u^b_1$ and $u^b_2$ such that

with $h^b_1, h^b_2>0$. Next Proposition provides the results about solutions, while the illustrating figures are postponed to the Appendix.

Proposition 2. System 45 admits a solution if the two regions $\mathcal{N}^B$ and $\mathcal{P}^A$ intersect in the subcritical set $\{(h, q): h>0, \ \tilde{\mathcal{C}}^-(h) \leq q \leq \tilde{\mathcal{C}}^+(h)\}$ or if $u_l\in\mathcal{P}^A{(u_r)}$.

Proof. We distinguish two cases:

Case 1. The two curves $h\mathcal{S}_1(h_l^*, v_l^*;h)$ and $\tilde{\phi}_r(h)$ intersect inside the subcritical region (Figure 15 of the Appendix). If $u_l\notin\mathcal{P}^A(u_r)$ the solution is trivially the intersection point. On the contrary, if $u_l\in\mathcal{P}^A(u_r)$ we have two possible solutions: the intersection point inside the subcritical region or the starting supercritical value $u_l$. So, we may obtain both fluvial or torrential regime. Coherently with the solution that would be obtained in a single channel we choose the fluvial regime assigning as solution the intersection point inside the subcritical region.

Case 2. The two curves $h\mathcal{S}_1(h_l^*, v_l^*;h)$ and $\tilde{\phi}_r(h)$ do not intersect inside the subcritical region. We distinguish the two subcases shown in Figure 16 and in Figure 17 of the Appendix.

Case 2.1. Referring to Figure 16, the point $u_l$ may fall or not in the $\mathcal{P}^B(u_r)$ region. If $u_l\in\mathcal{P}^A(u_r)$ the only possible solution is $u_l$ self and $u^b_1 = u^b_2 = u_l$; if on the contrary $u_l\notin\mathcal{P}^A(u_r)$, system 45 does not admit a solution.

Case 2.2. Referring Figure 17, if $u_l\notin\mathcal{P}^A(u_r)$ the only possible solution is the critical point $u^-_{r, \mathcal{R}}$ defined in 24, then $u^b_1 = u^b_2 = u^-_{r, \mathcal{R}}$; If $u_l\in\mathcal{P}^A(u_r)$, both $u_l$ self and the critical point $u^-_{r, \mathcal{R}}$ are admissible solutions. Coherently with the solution obtained in the previous subcase, we select again $u_l$ as solution and we assign $u^b_1 = u^b_2 = u_l$. So, the torrential regime propagates on the outgoing canal.

Case B$\rightarrow$B (Torrential $\rightarrow$ Torrential). As shown in Figure 18 of the Appendix, the two admissible regions $\mathcal{N}^B$ and $\mathcal{P}^B$ may have empty intersection. So, in this case there exists a solution at the junction which verifies the conditions of conservation of the mass and of equal heights, if and only if $u_l\in\mathcal{P}^B(u_r)$, i.e. $u^b_1 = u^b_2 = u_l$.

5.1.   Other possible conditions at the junction

Assuming different conditions at the junction give rise to new possible solutions. In canals network problems, it is usual to couple the conservation of the mass with the conservation of energy at the junctions. The specific energy $E$ is a useful parameter in channel flow and it is defined as

For a given flow rate, there are usually two states possible for the same specific energy. Studying $E$ as a function of $h$ for constant $q$, there is a minimum value of $E$ at a certain value of $h$ called the critical depth,

Critical depth $h_c$ corresponds to some critical channel velocity $v_c$ defined by $\mathcal{F}_c = 1$. For $E < E_{min}$ solution does not exists, and thus such a flow is impossible physically. For $E > E_{min}$ solutions are possible: large depth with $|\mathcal{F}|< 1$ subcritical flow, and small depth with $|\mathcal{F}| > 1$ supercritical flow.

In our case, assuming equal energy at the junction gives

Moreover, assuming $q^b_1 = q^b_2 = q^b$ to be constant (and known) we get from 48

where $\mathcal{F}_1 = v_1/\sqrt{g h_1}$ and where we used the following relations

Then, we have two possible solution for the heights values at the junction:

or

So, for $h^b_1\ne h^b_2$ we get new possible solutions at the junction with $(h^b_1, q^b)$ subcritical and $(h^b_2, q^b)$ supercritical or vice-versa. Specifically, case Torrential$\rightarrow$Fluvial and case Torrential$\rightarrow$Torrential may admit solution even if their admissible regions have empty intersection in the subcritical set.

Remark 3. In the case of a simple junction, the natural assumption (consistent with the dynamic of shallow-water equations) should be to assume the conservation of the momentum. With our notation, the relation 49 or 48 sholud be replaced by the following:

By the same reasoning used before in the case of the conservation of energy, from 51 we get

Then, we have again two possible relations for the heights values at the junction:

or

So again, for $h^b_1\ne h^b_2$ we get new possible solutions at the junction. Case Torrential$\rightarrow$ Fluvial and case Torrential$\rightarrow$Torrential may have solution and specifically we have that in Figure 16, points $(h^*_l, q_l)$ and $(h_l, q_l)$ verify 53 (see 23) and if $u_l\in\mathcal{P}^A(u_r)$ they are candidate to be the new solution. Even in the case described in Figure 18, if $u_l\in\mathcal{P}^B(u_r)$ the points $(h^*_l, q_l)$ and $(h_l, q_l)$ are the solution at the junction.

Let us conclude observing that for appropriate values of $(h_l, q_l)$, for Torrential$\rightarrow$ Fluvial we would get the same solution considering our simple network as a simple canal, i.e. a stationary shock called hydraulic jump characterized indeed by the conservation of the momentum in the transition from a supercritical to subcritical flow [11].

6.   Numerical tests

In this Section we illustrate the results of Section 5 by means of numerical simulations. We first give a sketch of the adopted numerical procedure and then we focus on two numerical tests which illustrate the regime transitions from fluvial to torrential and viceversa. The latter depend on well chosen initial conditions for Riemann problems at the junction.

We consider again a network formed by two canals intersecting at one single point, which represents the junction. Following [6], we use a high order Runge-Kutta Discontinuous Galerkin scheme to numerically solve system 3 on both canals 1 and 2:

The 1D domain of each canal is discretized into cells $C_m = [x_{m-\frac12}, x_{m+\frac12}]$, with $x_m = (x_{m-\frac12}+x_{m+\frac12})/2$, $m = 1, \ldots, M$, being $M$ the total number of computational cells. The DG method for 3 is formulated by multiplying the equation system by some test functions $w$, integrating over each computational cell, and performing integration by parts. Specifically, we seek the approximation $U = (u_1, u_2)$ with $u_i\in{W}_{\Delta x} = \{w: w|_{C_m}\in P^k(C_m), m = 1, \ldots, M\}$, $i = 1, 2$, where $P^k(C_m)$ is the space of polynomials of degree at most $k$ on cell $C_m$, such that $\forall\ w\in W_{\Delta x}$

Terms $w^{\pm}_{m+\frac12}$ denote left and right limits of the function values and the numerical fluxes $\hat f_{m\pm\frac12} = f(U^-_{m\pm\frac12}, U^+_{m\pm\frac12})$ are approximate Riemann solvers. In our simulations, we use the Lax-Friedrich flux. For implementation, each component of the approximate solution $U$, e.g. $u_1$, on mesh $C_m$ is expressed as

where $\{\psi_m^l(x)\}_{l = 0}^k$ is the set of basis functions of $P^k(C_m)$. Specifically, we choose the Legendre polynomials as local orthogonal basis of $P^k(C_m)$ and take the test function $w(x)$ in eq. 55 exactly as the set of basis functions $\psi_m^l(x)$, $l = 0 \cdots k$, assuming the polynomial degree $k = 2$. The equation system 55 can then be evolved in time via the method of lines approach by a TVD RK method. More implementation details for RKDG methods can be found in the original paper [7] and the review article [8].

Once the numerical procedure on both canals has been settled, the two systems in 55 have to be coupled with boundary conditions. At the junction the boundary values is settled as follows: at each time step and at each RK stage via the method-of-line approach, we set as left state in 42 (or 45) the approximate solution from canal 1 at the left limit of the junction, i.e.

with $[x_{M-\frac12}, x_{M+\frac12}]$ being the right-most cell in the 1D discretization of the incoming canal 1, and as right state in 42 (or 45) the approximate solution from canal 2 at the right limit of the junction, i.e.

with $[x_{\frac12}, x_{\frac32}]$ being the left-most cell in the 1D discretization of the outgoing canal 2. With this two values, we compute the intermediate states $u^b_1$ and $u^b_2$ by solving 42 (or 45) and, preserving the mass we directly assign the numerical fluxes at the junction as

Finally, in our simulations we assume Neumann boundary conditions at the free extremity of the channels.

Applying this numerical procedure, in Figure 10 and 11 we give two examples which illustrate the solution that is obtained in the regime transitions from fluvial to torrential and viceversa. In Figure 10, we assume to have a starting configuration given by the following subcritical constant states: $u_l = (0.25, 0.025)$ on canal 1 (left) and $u_r = (2.5, 0.25)$ on canal 2 (right). We are in the situation showed in Figure 9 and the boundary value at the junction is defined by the critical value $u^-_{r\mathcal{R}}\approx(1.088, -3.55)$. The solution is a backward water movement along canal 1 with a torrential regime. In Figure 11, we assume as initial state on canal 1 the supercritical constant value $u_l = (0.2, 3)$ and on canal 2 the subcritical constant value $u_r = (1.8, 4)$. We are in the configuration given in Figure 16 and the boundary value at the junction is defined by the starting point $u_l$. The solution is a forward movement with positive velocity, so the torrential regime propagates on canal 2.

7.   Conclusions

This paper deals with open canal networks. The interest stems out of applications such as irrigation channels water management. We base our investigations on the well-known Saint-Venant or shallow water equations. Two regimes exist for this hyperbolic system of balance laws: the fluvial, corresponding to eigenvalues with different sign, and the torrential, corresponding to both positive eigenvalues. Most authors focused the attention on designing and analysing network dynamics for the fluvial regime, while here we extend the theory to include regime transitions. After analyzing the Lax curves for incoming and outgoing canals, we provide admissibility conditions for Riemann solvers, describing solutions for constant initial data on each canal. Such analysis allows to define uniquely dynamics according to a set of conditions at junctions, such as conservation of mass, equal water height or equal energy. More precisely, the simple case of one incoming and outgoing canal is treated showing that, already in this simple example, regimes transitions appear naturally at junctions. Our analysis is then visualized by numerical simulations based on Runge-Kutta Discontinuous Galerkin methods.

Acknowledgments

M. Briani is a member of the INdAM Research group GNCS.

Appendix

Here we collect additional figures illustrating attainable regions for half-riemann problems and solutions for a simple channel. Figures 12-14 refer to the right-half Riemann problem described in Section 4.2. They show the regions of admissible states such that waves on the outgoing canals do not propagate into the junction, given a right state $u_r$ such that $|\tilde{\mathcal{F}}_r|<1$, $\tilde{\mathcal{F}}_r>1$ and $\tilde{\mathcal{F}}_r<-1$ respectively.

Figures 15-18 refer to Section 5 in which we study the possible solutions at a simple junction for different flow regimes, assuming the conservation of mass and equal heights at the junction. Specifically, Figures 15-17 illustrate the possible configurations and their associated solution that may occur during the transition from torrential to fluvial regime. The last Figure 18 shows instead the only possible configuration that admits a solution for the torrential flow regime.