<figure><img src="image_1.png" /></figure>

HAL
open science

# Second-order well-balanced Lagrange-projection schemes for blood flow equations

Alessia Del Grosso, Christophe Chalons

► To cite this version:

Alessia Del Grosso, Christophe Chalons. Second-order well-balanced Lagrange-projection schemes for blood flow equations. Calcolo, 2021, 58 (4), 10.1007/s10092-021-00434-5. http-03993025

HAL Id: hal-03993025

https://hal.science/hal-03993025

Submitted on 28 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Copyright
# Second-order well-balanced Lagrange-Projection schemes for Blood Flow Equations

A. Del Grosso* and C. Chalons†

December 7, 2020

**Abstract.** We focus on the development of well-balanced Lagrange-projection schemes applied to the one-dimensional blood flow system of balance laws. Here we neglect the friction forces and the source term is due to the presence of varying parameters as the cross-sectional area at the equilibrium and the arterial stiffness. By well-balanced we mean that the method preserves the "man at eternal rest" solution. For this purpose we present two different strategies: the former requires a consistent definition of the source term based on an approximate Riemann solver, while the second one exploits the well-established hydrostatic reconstruction. Subsequently we explain how to reach the second-order of accuracy for both procedures. Numerical simulations are carried out in order to show the right order of accuracy and the good behaviour of the schemes.

## 1 Introduction

This paper focuses on the construction of second-order well-balanced Lagrange-projection schemes applied to the 1D Blood Flow Equations (BFE). This model results to be extremely useful when dealing with the study of the cardiovascular system and related diseases. Indeed, it proved to be effective in the computation of averaged quantities as the cross-sectional area $A$ of the vessel, the blood flow $q$ and internal pressure $p$. Hence, there is a huge amount of works about this system, and we refer the reader to [17, 33] and the references therein for details about it. Here we study the model as applied to arteries, in the particular case in which the cross-sectional area at equilibrium and the wall stiffness could be not constant. See for instance [15, 20, 36, 31, 32]. Indeed, there exist physiological and pathological situations in which geometrical and mechanical parameters can vary locally, as in presence of stenoses or aneurysm and tapering of blood vessels. However, to consider non-constant parameters leads to the presence of a non-null source term and consequently we aim to develop numerical schemes for hyperbolic system of balance laws.

We are also interested in preserving the so-called "man at eternal rest" stationary solution, namely in the well-balancedness of the numerical method. As a matter of fact, it reveals itself to be an important property as a non well-balanced scheme could produce non-physical spurious oscillations in certain cases, especially when the solution is near to a steady state. In particular, the "man at eternal rest" condition is characterized by zero-velocity; if the numerical scheme preserves also all the stationary solutions, it is called fully well-balanced. Many studies have been done about well-balanced methods as in [4, 6, 7, 8], while for application to hyperbolic systems as the shallow water equations see [1, 3, 21, 22]. As far as the blood flow equations are concerned, we refer for instance to the works of Delestre and collaborators [15], in which variations in the values of parameters as the cross-sectional area at equilibrium is considered. In particular they developed a first-order well-balanced scheme basing themselves on the well-known hydrostatic reconstruction procedure, introduced for the first time by Audusse et al. in [1] in the context of the shallow water equations. In [20] Delestre et al. expanded their work considering varying values for the arterial wall rigidity as well. In [26] Müller et al. followed the generalized hydrostatic reconstruction to build a high-order well-balanced path-conservative numerical method for blood flow equations with varying mechanical properties. Then, once again in [27] Müller and Toro presented a high-order well-balanced path-conservative

<footnote>*Laboratoire de Mathématiques de Versailles, UMR 8100, Université de Versailles Saint-Quentin-en-Yvelines, UFR des Sciences, bâtiment Fermat, 45 avenue des Etats-Unis, 78035 Versailles cedex, France. E-mail: alessia.del-grosso@ens.uvsq.fr</footnote>

<footnote>†Laboratoire de Mathématiques de Versailles, UMR 8100, Université de Versailles Saint-Quentin-en-Yvelines, UFR des Sciences, bâtiment Fermat, 45 avenue des Etats-Unis, 78035 Versailles cedex, France. E-mail: christophe.chalons@uvsq.fr</footnote>
numerical scheme for BFE but with also discontinuous values for the cross-sectional area at rest and external pressure and not only for the wall stiffness.

Last but not least, we design numerical schemes based on the Lagrange-projection formalism, which allows us to split up the system into two different ones, and in particular to take separately into account the acoustic (Lagrangian step) and the transport (remap step) part of the system. This also implies the decoupling of fast and slow waves, which plays an important role in determining the CFL time step. Indeed, a fast wave leads to a very restrictive time step, contrarily to a slow wave. However, as far as the human arteries are concerned, the ratio between the velocity and the wave speed, i.e. the Shapiro number (or equivalently the Froude number for the shallow water equations), is in general of order $10^{-2}$ and thus, it could be interesting to develop a numerical method in which the Lagrangian step is solved implicitly. In this paper we start describing explicit schemes but we aim to develop implicit-explicit schemes in future works. For the Lagrangian-projection scheme we refer the reader to [12, 9, 10] and the references therein, while for well-balanced Lagrange-projection methods see [5, 11, 25].

Concluding, we present two second-order well-balanced Lagrange-projection schemes for the 1D BFE. We start introducing the numerical method in the case of constant parameters, and thus for a system of conservation laws. Then, when varying properties are considered, two different ways of preserving the "man at eternal rest" stationary solution are described. On the one hand, referring to the work of Suliciu [29], we relax the Lagrangian system introducing a new variable, which stands for a linearization of the pressure term, and then, following the theory of Gallice [18, 19], we easily solve the associated Riemann problem. Alternatively, we exploit the hydrostatic reconstruction as well.

*Outline of the paper.* In the next section we present the 1D mathematical model for the blood flow equations. We also introduce the Lagrange-projection decomposition which leads to two different systems, the acoustic and transport one. An approximate Riemann solver is described for the acoustic system as well. In sections 3 and 4 we respectively present the first and second-order well-balanced schemes. In both sections, we explain two different strategies in order to preserve the "man at eternal rest" solution. In section 5 numerical simulations are carried out. Finally conclusions and perspectives are drawn in section 6.

## 2 The mathematical model

Given the axial coordinate $x$ along the longitudinal axis of the vessel and the time $t > 0$, the general one-dimensional blood flow model consists of two equations, the mass conservation and momentum balance equation, namely

$$ \left\{ \begin{aligned} \partial_t A + \partial_x q &= 0 \\ \partial_t q + \partial_x \left( \hat{\alpha} \frac{q^2}{A} \right) + \frac{A}{\rho} \partial_x p &= f, \end{aligned} \right. \qquad (2.1) $$

where $A(x, t) > 0$ is the cross-sectional area of the vessel, $q = Au$ the blood flow, with $u(x, t)$ the averaged velocity of blood at cross section, and finally $p(x, t)$ is the averaged internal pressure at cross section. Furthermore, $\rho$ represents the blood density and it is assumed to be constant, while $\hat{\alpha}$ is determined by the velocity profile and it is considered to be flat in this work, thus we take $\hat{\alpha} = 1$. At last, $f$ accounts for the friction forces but we will neglect it in the rest of the paper. We assume that the initial area $A(x, t = 0)$ and initial velocity $u(x, t = 0)$ are given at time $t = 0$. For more details about the derivation of system (2.1), we refer to [17] and [33].

Since in system (2.1) there are three unknowns but only two equations, we need a closure condition, namely a tube law or more specifically a relation between the internal pressure and the cross-sectional area. In this paper we refer to [20] and we consider the blood vessels to be purely elastic arteries and, as such, the tube law reads

$$ p(x, t) = p_{\text{ext}} + K(x)(\sqrt{A(x, t)} - \sqrt{A_0(x)}), \qquad (2.2) $$

where $p_{\text{ext}}$ is the constant external pressure, $A_0(x)$ the cross-sectional area at equilibrium and $K(x)$ a parameter related to the arterial stiffness. In particular $K$ is a positive function depending on the vessel thickness $h_0(x)$ and the Young modulus $E(x)$, refer again to [20]. Note also that (2.2) is valid only for blood flow in arteries and not in veins, for a more general tube law refer for instance to [17, 33]. Equipped with closure condition (2.2), we can now show the system of balance laws we
will investigate in the rest of the paper, namely

$$ \left\{ \begin{aligned} \partial_t A + \partial_x q &= 0 \\ \partial_t q + \partial_x \left( \frac{q^2}{A} + \gamma A^{\frac{3}{2}} \right) &= s, \end{aligned} \right. \qquad (2.3) $$

with $\gamma = \frac{K}{3\rho}$ and

$$ s = s(A; A_0, K) = \frac{A}{\rho} \partial_x (K \sqrt{A_0}) - \frac{2A}{3\rho} \sqrt{A} \partial_x K. \qquad (2.4) $$

In compact form this system reads

$$ \partial_t \mathbf{Q} + \partial_x \mathbf{F}(\mathbf{Q}) = \mathbf{S}(\mathbf{Q}; A_0, K) \qquad (2.5) $$

where

$$ \mathbf{Q} = \begin{pmatrix} A \\ q \end{pmatrix}, \quad \mathbf{F}(\mathbf{Q}) = \begin{pmatrix} q \\ \frac{q^2}{A} + \gamma A^{\frac{3}{2}} \end{pmatrix} \quad \text{and} \quad \mathbf{S}(\mathbf{Q}; A_0, K) = \begin{pmatrix} 0 \\ s \end{pmatrix}. $$

Let us note that if both $K(x)$ and $A_0(x)$ are constant, then $s = 0$ and (2.5) is reduced to a system of conservation laws.

It is not difficult to see that the two eigenvalues of system (2.5) are $\lambda^{\pm} = u \pm c$, where $c$ is the wave speed defined by

$$ c = \sqrt{\frac{3}{2}\gamma\sqrt{A}}. $$

Consequently the convective part of (2.5) is strictly hyperbolic as $\lambda^{\pm}$ are real and distinct, namely as long as the vector of unknowns $\mathbf{Q}$ belongs to the phase space $\Omega = \{(A, Au)^t \in \mathbf{R}^2 | A > 0\}$. Finally, both the two characteristic fields are genuinely non-linear and the Riemann invariants associated with $\lambda^{\pm}$ are respectively given by $I^{-} = u+4c$ and $I^{+} = u-4c$. For more details refer to [35].

In this paper we are specially interested in developing second-order well-balanced Lagrange-projection methods, thus hereafter we introduce both the well-balanced property and the Lagrange-projection decomposition.

**The well-balanced property.** A numerical scheme is well-balanced if it is able to preserve the smooth stationary solutions of the system, that is to say the steady states which satisfy the ordinary differential equations

$$ \partial_x \mathbf{F}(\mathbf{Q}) = \mathbf{S}(\mathbf{Q}; A_0, K), $$

and hence

$$ q = q_0 = \text{constant}, \quad \frac{q_0^2}{2A^2} + \frac{K}{\rho}(\sqrt{A} - \sqrt{A_0}) = \text{constant}, \qquad (2.6) $$

where the quantity $E = \frac{q^2}{2A^2} + \frac{K}{\rho}(\sqrt{A} - \sqrt{A_0})$ can be referred as the energy discharge. In particular, a scheme able to preserve the steady states (2.6) is called fully well-balanced, while a method which conserves only the stationary solutions with zero velocity ($u = 0$) is defined well-balanced. We are interested in a scheme endowed with the latter property, and thus in the "man at eternal rest" solution,

$$ q = 0, \quad K(\sqrt{A} - \sqrt{A_0}) = \text{constant}. \qquad (2.7) $$

For more details about well-balanced schemes for blood flow equations refer for instance to [15, 20, 27, 36].

**The Lagrangian coordinates.** Observing that system (2.3) is given in Eulerian coordinates, now we want to express it using the Lagrangian coordinates, which describe the flow following the fluid motion. While with the Eulerian coordinates the viewer has a fixed position and watches the flow from the exterior, with the Lagrangian coordinates he focus on a single "fluid particle" $\xi$, for which we introduce the characteristic curves

$$ \left\{ \begin{aligned} \frac{\partial x}{\partial t}(\xi, t) &= u(x(\xi, t), t) \\ x(\xi, 0) &= \xi. \end{aligned} \right. \qquad (2.8) $$
Then, given the trajectory: $t \to x(\xi, t)$, any function: $(x, t) \to Q(x, t)$ in Eulerian coordinates can be written in Lagrangian coordinates,

$$ \bar{Q}(\xi, t) = Q(x(\xi, t), t). $$

Moreover, defining the volume ratio

$$ L(\xi, t) = \frac{\partial x}{\partial \xi}(\xi, t) \qquad (2.9) $$

which satisfies

$$ \left\{ \begin{array}{l} \frac{\partial L}{\partial t}(\xi, t) = \partial_{\xi} u(x(\xi, t), t) \\ L(\xi, 0) = 1 \end{array} \right. \qquad (2.10) $$

and

$$ \partial_t L(\xi, t) = \partial_{\xi} u(x(\xi, t), t) = \partial_{\xi} \bar{u}(\xi, t), $$

we easily find that

$$ \partial_{\xi} \bar{Q}(\xi, t) = L(\xi, t) \partial_x Q(x, t) \quad \text{and} \quad \partial_t \bar{Q}(\xi, t) = \partial_t Q(x, t) + u(x, t) \partial_x Q(x, t). $$

Hence, using the chain rule and defining $\tilde{p} = \tilde{p}(A; K) = \gamma A^{\frac{3}{2}}$, from system (2.3) we write

$$ \left\{ \begin{array}{l} \partial_t A + A \partial_x u + u \partial_x A = 0 \\ \partial_t (Au) + u \partial_x (Au) + Au \partial_x u + \partial_x \tilde{p} = s, \end{array} \right. \qquad (2.11) $$

and multiplying by $L(\xi, t)$, we obtain

$$ \left\{ \begin{array}{l} L \partial_t \bar{A} + \bar{A} \partial_t L = 0 \\ L \partial_t (\bar{A}u) + \bar{A}u \partial_t L + \partial_{\xi} \tilde{p} = \bar{s}, \end{array} \right. $$

where $\bar{s} = \frac{A}{\rho} \partial_{\xi} (\overline{K \sqrt{A_0}}) - \frac{2A}{3\rho} \sqrt{A} \partial_{\xi} \bar{K}$. We finally find that in Lagrangian coordinates system (2.3) reads,

$$ \left\{ \begin{array}{l} \partial_t (L \bar{A}) = 0 \\ \partial_t (L \bar{A}u) + \partial_{\xi} \tilde{p} = \bar{s}. \end{array} \right. \qquad (2.12) $$

Hence, the Lagrange-projection algorithm consists of two steps:

1. Solve system (2.12) written in Lagrangian coordinates;

2. Project the solution of system (2.12) in Eulerian coordinates.

For more details about the Lagrangian-projection decomposition, once again we refer the reader to [5, 9, 10, 11, 16, 25].

## 2.1 The Lagrangian-projection splitting and the relaxation formulation

At this stage, we present the Lagrangian-projection decomposition in a different way, which will prove to be extremely useful for one of the two numerical methods we are going to describe in sections 3.2.1 and 4.2.1. In particular we split system (2.3) into two different ones, the Lagrangian/acoustic system and the projection/transport system. The former takes into account the acoustic effects and parameters variations, while the latter the transport phenomena, see [5, 11].

In particular, starting from formulation (2.11), we find that the acoustic system reads

$$ \left\{ \begin{array}{l} \partial_t A + A \partial_x u = 0 \\ \partial_t (Au) + Au \partial_x u + \partial_x \tilde{p} = s, \end{array} \right. \qquad (2.13) $$

while the transport system is given by

$$ \left\{ \begin{array}{l} \partial_t A + u \partial_x A = 0 \\ \partial_t q + u \partial_x q = 0. \end{array} \right. \qquad (2.14) $$
System (2.14) can also simply be seen as $\partial_t X + u \partial_x X = 0$, with either $X = A$ or $X = q$.

Let us rewrite system (2.13) as

$$ \begin{cases} -\frac{1}{A^2}(\partial_t A + A\partial_x u) = 0 \\ A\partial_t u + u\partial_t A + Au\partial_x u + \partial_x \tilde{p} = s \end{cases} $$

and defining $\tau = \frac{1}{A}$ and the mass variable $m$ such that $\frac{1}{A}\partial_x = \partial_m$, the Lagrangian system also reads

$$ \begin{cases} \partial_t \tau - \partial_m u = 0 \\ \partial_t u + \partial_m \tilde{p} = \tilde{s}, \end{cases} \qquad (2.15) $$

with $\tilde{s} = \frac{A}{\rho}\partial_m(K\sqrt{A_0}) - \frac{2A}{3\rho}\sqrt{A}\partial_m K$. Observe that (2.15) is equivalent to (2.12). System (2.15) has two eigenvalues $\lambda_{\pm} = \pm Ac$, and it is strictly hyperbolic in the same phase space of system (2.3), namely when $A > 0$, with the two characteristic fields genuinely non-linear.

**Relaxation formulation.** At this stage we are interested in finding an approximate solution of a Riemann problem for system (2.15). For this purpose we exploit the Suliciu relaxation approach, which allows us to enlarge (2.15) to a strictly hyperbolic system with only linearly degenerate characteristic fields, which is well-known to be easier to solve. For the Suliciu relaxation approach and related applications refer to [29, 4, 2, 12, 13, 14] and the references therein.

Thus, we introduce the relaxation parameter $\lambda$ and the new variable $\Pi$ such that at least formally

$$ \lim_{\lambda \to \infty} \Pi = \tilde{p}, $$

where $\Pi$ can be interpreted as a linearization of the pressure $\tilde{p}$. We observe that $\partial_{\tau}\tilde{p}(\tau) = -\frac{1}{A^2}\partial_A\tilde{p}(A)$ and $\partial_t\tilde{p} = \partial_{\tau}\tilde{p}(\tau)\partial_t\tau$ so that, multiplying the first equation of system (2.15) by $\partial_{\tau}\tilde{p}(\tau)$, we find that $\partial_t\tilde{p} + A^2c^2\partial_m u = 0$. The latter motivates the relaxation system

$$ \begin{cases} \partial_t \tau - \partial_m u = 0 \\ \partial_t u + \partial_m \Pi = \tilde{s} \\ \partial_t \Pi + a^2 \partial_m u = \lambda(\tilde{p}(\tau) - \Pi) \end{cases} \qquad (2.16) $$

where $a^2$ is a constant which linearizes $A^2c^2$ and that should be taken such that $a^2 \ge A^2c^2$ according to the sub-characteristic condition. Indeed, this condition entails that the information in the relaxation model (2.16) propagates faster than in the original one (2.15). Refer for instance to the works [29, 11, 12].

Considering that the initial data for $\Pi$ is well-prepared in the sense that $\Pi = \tilde{p}$, it is natural to introduce a more compact notation for system (2.16), which reads

$$ \partial_t \mathbf{U} + \partial_m \mathbf{G}(\mathbf{U}) = \tilde{\mathbf{S}} $$

with

$$ \mathbf{U} = \begin{pmatrix} \tau \\ u \\ \Pi \end{pmatrix}, \quad \mathbf{G}(\mathbf{U}) = \begin{pmatrix} -u \\ \Pi \\ a^2 u \end{pmatrix} \quad \text{and} \quad \tilde{\mathbf{S}} = \begin{pmatrix} 0 \\ \tilde{s} \\ 0 \end{pmatrix}, $$

and where we note that the relaxation source term in the evolution equation for $\Pi$ is not present anymore. Rewriting this system in quasi-linear form, it gives

$$ \partial_t \mathbf{U} + \mathbf{A}(\mathbf{U})\partial_m \mathbf{U} = \tilde{\mathbf{S}} $$

where $\mathbf{A}(\mathbf{U})$ is the Jacobian matrix of the flux vector $\mathbf{G}(\mathbf{U})$, that is

$$ \mathbf{A}(\mathbf{U}) = \frac{\partial \mathbf{G}}{\partial \mathbf{U}} = \begin{pmatrix} 0 & -1 & 0 \\ 0 & 0 & 1 \\ 0 & a^2 & 0 \end{pmatrix}. $$

The eigenvalues of the matrix $\mathbf{A}(\mathbf{U})$ are

$$ \lambda_{-} = -a, \quad \lambda_{0} = 0 \quad \text{and} \quad \lambda_{+} = a; $$
note that they can be seen as constant approximations of the eigenvalues of system (2.15). Therefore system (2.16) is strictly hyperbolic as long as a is real and $a \neq 0$. The corresponding right eigenvectors are

$$ \mathbf{R}_{-} = \begin{pmatrix} 1 \\ a \\ -a^2 \end{pmatrix}, \quad \mathbf{R}_{0} = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} \quad \text{and} \quad \mathbf{R}_{+} = \begin{pmatrix} 1 \\ -a \\ -a^2 \end{pmatrix}. $$

All the three associated characteristic fields are linearly degenerate, therefore the three waves will be contact discontinuities and, as such, the Riemann problem can be solved using the Riemann invariants. The equations in phase space associated to $\lambda_{\pm}$ are

$$ \frac{d\tau}{1} = \frac{du}{\mp a} = \frac{d\Pi}{-a^2} $$

and they lead to the Riemann invariants $RI_{\pm,1} = a\tau \pm u$ and $RI_{\pm,2} = \Pi \mp au$. Finally, the corresponding equations to $\lambda_0$ in phase space are

$$ du = 0 \quad \text{and} \quad d\Pi = 0 $$

and the associated Riemann invariants are $RI_{0,1} = u$ and $RI_{0,2} = \Pi$.

At this stage we can look for an approximate solution of system (2.16) for a Riemann problem. We will use the theory of Gallic [18, 19] which consists in an extension of the well-known Harten, Lax and van Leer formalism [24] for systems of conservation laws.

## 2.2 Approximate solution of a RP for the acoustic system

Let us start giving some general notions about the approximate Riemann solver and consistency with the integral form as described by Gallic [18, 19]. We will then focus on our specific case, namely system (2.16).

We briefly consider a general system of form

$$ \partial_t \mathbf{U} + \partial_m \mathbf{G}(\mathbf{U}) = \tilde{\mathbf{S}} \quad (2.17) $$

where $\mathbf{U}$ is the vector of $h$ unknowns, $\mathbf{G}(\mathbf{U})$ the physical flux and $\tilde{\mathbf{S}}$ the source term. We want to solve the Riemann Problem (RP) with the following initial condition,

$$ \mathbf{U}(m, t=0) = \begin{cases} \mathbf{U}_L & \text{if } m < 0 \\ \mathbf{U}_R & \text{if } m > 0, \end{cases} \quad (2.18) $$

for any given $\mathbf{U}_L, \mathbf{U}_R$ in the phase space. Supposing to have $h$ discontinuities with velocities $\lambda_k, k = 1, \dots, h$, the solution of the Riemann problem consists of $h+1$ states separated by the discontinuities, hence

$$ \mathbf{U}\left(\frac{m}{t}, \mathbf{U}_L, \mathbf{U}_R\right) = \begin{cases} \mathbf{U}_1 = \mathbf{U}_L & \text{if } \frac{m}{t} < \lambda_1 \\ \vdots & \vdots \\ \mathbf{U}_k & \text{if } \lambda_k < \frac{m}{t} < \lambda_{k+1} \\ \vdots & \vdots \\ \mathbf{U}_{h+1} = \mathbf{U}_R & \text{if } \frac{m}{t} > \lambda_h \end{cases} \quad (2.19) $$

Given the space and time steps $\Delta m$ and $\Delta t$, the approximate solution (2.19) of the RP (2.17)-(2.18) is said to be consistent with the integral form of (2.17) if in the interval $[0, \Delta m]$ we have

$$ \mathbf{G}(\mathbf{U}_R) - \mathbf{G}(\mathbf{U}_L) - \Delta m \tilde{\mathbf{S}}(\Delta m, \Delta t; \mathbf{U}_L, \mathbf{U}_R) = \sum_{k=1}^{h} \lambda_k (\mathbf{U}_{k+1} - \mathbf{U}_k), \quad (2.20) $$

where also the source term $\tilde{\mathbf{S}}(\Delta m, \Delta t; \mathbf{U}_L, \mathbf{U}_R)$ has to satisfy a consistency property, namely

$$ \lim_{\mathbf{U}_L, \mathbf{U}_R \to \mathbf{U}; \Delta m, \Delta t \to 0} \tilde{\mathbf{S}}(\Delta m, \Delta t; \mathbf{U}_L, \mathbf{U}_R) = \tilde{\mathbf{S}}(\mathbf{U}). \quad (2.21) $$
### 2.2.1 Riemann solver for the system of conservation laws

Next, we focus on system (2.16). We start neglecting the source term, hence we want to solve the following Riemann problem,

$$ \left\{ \begin{array}{l} \partial_t \mathbf{U} + \partial_m \mathbf{G}(\mathbf{U}) = 0 \\ \mathbf{U}(m, t=0) = \begin{cases} \mathbf{U}_L & \text{if } m < 0 \\ \mathbf{U}_R & \text{if } m > 0 \end{cases} \end{array} \right. \qquad (2.22) $$

where in particular

$$ \mathbf{U}_L = \begin{pmatrix} \tau_L \\ u_L \\ \Pi_L \end{pmatrix} \quad \text{and} \quad \mathbf{U}_R = \begin{pmatrix} \tau_R \\ u_R \\ \Pi_R \end{pmatrix}. $$

In this specific case, the solution can be computed exactly and takes the following form

$$ \mathbf{U}\left(\frac{m}{t}, \mathbf{U}_L, \mathbf{U}_R\right) = \begin{cases} \mathbf{U}_L & \text{if } \frac{m}{t} < \lambda_- = -a \\ \mathbf{U}_{*,L} & \text{if } \lambda_- < \frac{m}{t} < 0 \\ \mathbf{U}_{*,R} & \text{if } 0 < \frac{m}{t} < \lambda_+ = a \\ \mathbf{U}_R & \text{if } \frac{m}{t} > \lambda_+ \end{cases} \qquad (2.23) $$

with

$$ \mathbf{U}_{*,L} = \begin{pmatrix} \tau_{*,L} \\ u_{*,L} \\ \Pi_{*,L} \end{pmatrix} \quad \text{and} \quad \mathbf{U}_{*,R} = \begin{pmatrix} \tau_{*,R} \\ u_{*,R} \\ \Pi_{*,R} \end{pmatrix}. \qquad (2.24) $$

In order to find the values $\mathbf{U}_{*,L}$ and $\mathbf{U}_{*,R}$, we exploit the consistency conditions (2.20) together with the Rankine-Hugoniot conditions associated to the velocities $\pm a$, namely

$$ a(\mathbf{U}_{*,L} - \mathbf{U}_L) + \mathbf{G}(\mathbf{U}_{*,L}) - \mathbf{G}(\mathbf{U}_L) = 0 \quad \text{and} \quad -a(\mathbf{U}_R - \mathbf{U}_{*,R}) + \mathbf{G}(\mathbf{U}_{*,R}) - \mathbf{G}(\mathbf{U}_R) = 0, \qquad (2.25) $$

and to the discontinuity of null velocity, i.e. $u_{*,L} = u_{*,R} = u_*$, $\Pi_{*,L} = \Pi_{*,R} = \Pi_*$. Thus, it is straightforward to find the following algebraic system of 6 relations,

$$ \left\{ \begin{array}{l} a\tau_{*,L} - u_{*,L} = a\tau_L - u_L \\ au_{*,L} + \Pi_{*,L} = au_L + \Pi_L \\ u_{*,L} = u_{*,R} = u_* \\ \Pi_{*,L} = \Pi_{*,R} = \Pi_* \\ a\tau_{*,R} + u_{*,R} = a\tau_R + u_R \\ au_{*,R} - \Pi_{*,R} = au_R - \Pi_R, \end{array} \right. $$

whose resolution leads us to the values in the star regions,

$$ \left\{ \begin{array}{l} \tau_{*,L} = \tau_L + \frac{1}{a}(u_* - u_L) = \tau_L + \frac{1}{2a}(u_R - u_L) - \frac{1}{2a^2}(\Pi_R - \Pi_L) \\ \tau_{*,R} = \tau_R - \frac{1}{a}(u_* - u_R) = \tau_R + \frac{1}{2a}(u_R - u_L) + \frac{1}{2a^2}(\Pi_R - \Pi_L) \\ u_* = \frac{1}{2}(u_L + u_R) - \frac{1}{2a}(\Pi_R - \Pi_L) \\ \Pi_* = \frac{1}{2}(\Pi_L + \Pi_R) - \frac{a}{2}(u_R - u_L). \end{array} \right. \qquad (2.26) $$

### 2.2.2 Riemann solver for the system of balance laws

Let us now consider the source term due to the variations of the cross-sectional area at equilibrium $A_0$ and the arterial stiffness $K$ as well. Therefore, we look for an approximate solution of the following Riemann Problem,

$$ \left\{ \begin{array}{l} \partial_t \mathbf{U} + \mathbf{A}(\mathbf{U})\partial_m \mathbf{U} = \tilde{\mathbf{S}}(\mathbf{Q}; A_0, K) \\ \mathbf{U}(m, t=0) = \begin{cases} \mathbf{U}_L & \text{if } m < 0 \\ \mathbf{U}_R & \text{if } m > 0. \end{cases} \end{array} \right. \qquad (2.27) $$
Once again we assume that the $m-t$ plane divided in four different zones by the three waves and the solution has the form given by (2.23) and (2.24). In order to find such a solution made of 6 unknowns, we need to impose different conditions. Starting from the consistency relations (2.20), we have

$$ egin{cases} -u_R + u_L = -a(\tau_{*,L} - \tau_L) + a(\tau_R - \tau_{*,R}) \\ \Pi_R - \Pi_L - \Delta m \tilde{s} = -a(u_{*,L} - u_L) + a(u_R - u_{*,R}) \\ a^2(u_R - u_L) = -a(\Pi_{*,L} - \Pi_L) + a(\Pi_R - \Pi_{*,R}) \end{cases} \quad (2.28) $$

then, the Rankine-Hugoniot conditions associated to the mass equation read

$$ egin{cases} -(u_R - u_{*,R}) = a(\tau_R - \tau_{*,R}) \\ u_L - u_{*,L} = -a(\tau_{*,L} - \tau_L) \\ u_{*,L} = u_{*,R} = u_* \end{cases} \quad (2.29) $$

which give us only two additional conditions, as the first relation of (2.28) is a linear combination of (2.29). Consequently we have found only five conditions for six unknowns and thus, they are not sufficient to define the approximate RP solution. Note also that $\tilde{s} = \tilde{s}(\Delta m, \Delta t; U_L, U_R)$ in (2.28) has to be specified such that (2.21) holds true. In particular, $\tilde{s}$ should be determined such that it is equal to zero if both $A_0$ and $K$ are constant. With this request, the solution of (2.27) would degenerate to the solution of the homogeneous Riemann problem (2.22). Finally, two relations are still missing to define our solution.

Let us suppose to have two other equations, namely $\partial_t A_0 = 0$ and $\partial_t K = 0$. For these two equations the solutions are $A_{0,L}$ if $m < 0$, $A_{0,R}$ if $m > 0$ and $K_L$ if $m < 0$, $K_R$ if $m > 0$. Hence, we ask for the jump condition across the middle discontinuity associated to the momentum equation, that is to say

$$ \Pi_{*R} - \Pi_{*L} + \mathcal{M} = 0, \quad (2.30) $$

where the function

$$ \mathcal{M} = - \frac{\Delta m_L + \Delta m_R}{2} \tilde{s}(\Delta m, \Delta t; U_L, U_R) $$

has to be defined such that it satisfies $\mathcal{M} = 0$ if $A_{0,L} = A_{0,R}$ and $K_L = K_R$. Thus, after few algebraic computations, we obtain the following solution for the Riemann problem (2.27),

$$ egin{cases} \tau_{*,L} = \tau_L + \frac{1}{a}(u_* - u_L) \\ \tau_{*,R} = \tau_R - \frac{1}{a}(u_* - u_R) \\ u_* = \frac{1}{2}(u_L + u_R) - \frac{1}{2a}(\Pi_R - \Pi_L) - \frac{M}{2a} \\ \Pi_{*L} = \frac{1}{2}(\Pi_L + \Pi_R) - \frac{a}{2}(u_R - u_L) + \frac{M}{2} \\ \Pi_{*R} = \frac{1}{2}(\Pi_L + \Pi_R) - \frac{a}{2}(u_R - u_L) - \frac{M}{2} \\ \Pi_* = \frac{1}{2}(\Pi_{*L} + \Pi_{*R}) = \frac{1}{2}(\Pi_L + \Pi_R) - \frac{a}{2}(u_R - u_L). \end{cases} \quad (2.31) $$

It remains to define $\mathcal{M}$ in a consistent way and it is clear that

$$ \mathcal{M} = - \left\{ \frac{A}{\rho} \right\} \left( K_R \sqrt{A_{0,R}} - K_L \sqrt{A_{0,L}} \right) + \left\{ \frac{2A}{3\rho} \sqrt{A} \right\} (K_R - K_L) $$

is relevant provided that $\left\{ \frac{A}{\rho} \right\}$ and $\left\{ \frac{2A}{3\rho} \sqrt{A} \right\}$ are consistent approximations of $\frac{A}{\rho}$ and $\frac{2A}{3\rho} \sqrt{A}$ respectively. Let us discuss this issue. We first recall that since $U_L$ and $U_R$ are taken to be well prepared, we have $\Pi_L = \frac{K_L}{3\rho} A_L^{\frac{3}{2}}$ and $\Pi_R = \frac{K_R}{3\rho} A_R^{\frac{3}{2}}$. Since we are interested in preserving the "man at eternal rest" solution, we also ask for the well-balanced condition, $U_{*,L} = U_L$ and $U_{*,R} = U_R$ if $U_L$ and $U_R$ verify the "man at eternal rest" solution, namely

$$ u_L = u_R = 0, \quad K_L(\sqrt{A_L} - \sqrt{A_{0,L}}) = K_R(\sqrt{A_R} - \sqrt{A_{0,R}}). \quad (2.32) $$
If we have (2.32), then

$$
\begin{align*}
\mathcal{M} &= -\left\{\frac{A}{\rho}\right\}(K_R \sqrt{A_{0,R}} - K_L \sqrt{A_{0,L}}) + \left\{\frac{2A}{3\rho}\sqrt{A}\right\}(K_R - K_L) \\
&= -\left\{\frac{A}{\rho}\right\}(K_R \sqrt{A_R} - K_L \sqrt{A_L}) + \left\{\frac{2A}{3\rho}\sqrt{A}\right\}(K_R - K_L) \\
&= -\left\{\frac{A}{\rho}\right\}\frac{K_R + K_L}{2}(\sqrt{A_R} - \sqrt{A_L}) - \left\{\frac{A}{\rho}\right\}\frac{\sqrt{A_R} + \sqrt{A_L}}{2}(K_R - K_L) + \left\{\frac{2A}{3\rho}\sqrt{A}\right\}(K_R - K_L)
\end{align*}
$$

and thanks to (2.30), using formulae $[XY]_L^R = \vec{X}[Y]_L^R + \overleftarrow{Y}[X]_L^R$, where $[X]_L^R = X_R - X_L$, $\vec{X} = \alpha X_L + (1-\alpha)X_R$, $\overleftarrow{X} = \alpha X_R + (1-\alpha)X_L$ and $\alpha = \frac{1}{2}$,

$$
\begin{align*}
\mathcal{M} &= -(\Pi_R - \Pi_L) \\
&= -\frac{1}{3\rho}(K_R A_R \sqrt{A_R} - K_L A_L \sqrt{A_L}) \\
&= -\frac{1}{6\rho}(K_R + K_L)(A_R \sqrt{A_R} - A_L \sqrt{A_L}) - \frac{1}{6\rho}(A_L \sqrt{A_L} + A_R \sqrt{A_R})(K_R - K_L) \\
&= -\frac{1}{3\rho}\frac{K_R + K_L}{2}(A_L + \sqrt{A_R}\sqrt{A_L} + A_R)(\sqrt{A_R} - \sqrt{A_L}) - \frac{1}{6\rho}(A_L \sqrt{A_L} + A_R \sqrt{A_R})(K_R - K_L).
\end{align*}
$$

Therefore, it is clearly sufficient to set

$$
\left\{ \frac{A}{\rho} \right\} = \frac{A_L + A_R + \sqrt{A_L} \sqrt{A_R}}{3\rho}
$$

and

$$
\begin{align*}
\left\{ \frac{2A}{3\rho} \sqrt{A} \right\} &= \frac{\sqrt{A_L} + \sqrt{A_R}}{2} \left\{ \frac{A}{\rho} \right\} - \frac{1}{6\rho} (A_L \sqrt{A_L} + A_R \sqrt{A_R}) \\
&= \frac{\sqrt{A_L} \sqrt{A_R} (\sqrt{A_L} + \sqrt{A_R})}{3\rho}.
\end{align*}
$$

Finally, we get

$$
\begin{align}
\mathcal{M} &= \mathcal{M}((A_L, A_{0,L}, K_L); (A_R, A_{0,R}, K_R)) \nonumber \\
&= -\frac{A_L + A_R + \sqrt{A_L \sqrt{A_R}}}{3\rho} \left( K_R \sqrt{A_{0,R}} - K_L \sqrt{A_{0,L}} - \frac{\sqrt{A_L \sqrt{A_R}} (\sqrt{A_L} + \sqrt{A_R})}{A_L + A_R + \sqrt{A_L \sqrt{A_R}}} (K_R - K_L) \right). \tag{2.33}
\end{align}
$$

Lastly, let us note that the definition of $\mathcal{M}$ is consistent, indeed

$$
\lim_{A_L, A_R \to A} \mathcal{M} = -\frac{A}{\rho} \left( K_R \sqrt{A_{0,R}} - K_L \sqrt{A_{0,L}} \right) + \frac{2A\sqrt{A}}{3\rho} (K_R - K_L).
$$

3 First-order well-balanced scheme

The next step consists in presenting the first-order well-balanced scheme. We start giving a first-order approximation of
the homogeneous version of system (2.3) and then we modify it in order to satisfy the well-balanced property and include
the source term at the same time. In particular, for the latter step, we show two different ways of preserving the stationary
state (2.7), one of them exploiting the approximate solution of the Riemann problem (2.27) and the other the well-known
hydrostatic reconstruction procedure, for which we respectively refer to [11] and [25].

Let us now introduce some notations. First of all, we define the constant space step $\Delta x$ and constant time step $\Delta t$. The mesh interfaces are given by $x_{j+1/2} = j\Delta x$ for $j \in \mathbb{Z}$ and the intermediate times by $t^n = n\Delta t$ for $n \in \mathbb{N}$. As usual in the
finite volume framework, we seek at each time $t^n$ for an approximation $\mathbf{Q}_j^n$ of the solution in the interval $[x_{j-1/2}, x_{j+1/2})$, $j \in \mathbb{Z}$. Therefore, a piecewise constant approximate solution $x \to \mathbf{Q}_{\Delta t, \Delta x}(x, t^n)$ of the solution $\mathbf{Q}$ is given by

$$ \mathbf{Q}_{\Delta t, \Delta x}(x, t^n) = \mathbf{Q}_j^n \quad \text{for all } x \in C_j = [x_{j-1/2}; x_{j+1/2}), \quad j \in \mathbb{Z}, \quad n \in \mathbb{N}. $$

As far as the variable $\xi$ is concerned, we use the same space discretization of $x$, hence $\Delta x = \Delta\xi$, $x_{j+1/2} = \xi_{j+1/2}$ and $x_j = \xi_j \forall j$.

## 3.1 First-order approximation with constant parameters

We start describing the first-order scheme by assuming $A_0$ and $K$ constant parameters, therefore the source term $s$ disappears. As we already pointed out, the Lagrange-projection scheme is composed of two steps:

1. Solve system (2.12) written in Lagrangian coordinates, or equivalently update $\mathbf{Q}^n$ to $\mathbf{Q}^{n+1^-}$ approximating the solution of system (2.13);

2. Project the solution of system (2.12) in Eulerian coordinates, namely update $\mathbf{Q}^{n+1^-}$ to $\mathbf{Q}^{n+1}$ by solving system (2.14).

For more details about this procedure in the shallow water context refer to [5, 11, 25].

*Lagrangian step.* Given a system of form $\partial_t U + \partial_m G(U) = 0$ as in (2.22), the first-order Godunov-type scheme associated with the Riemann solver of section 2.2.1 reads

$$ \mathbf{U}_j^{n+1^-} = \mathbf{U}_j^n - \frac{\Delta t}{\Delta m_j} (\mathbf{G}_{j+\frac{1}{2}}^n - \mathbf{G}_{j-\frac{1}{2}}^n), \quad (3.1) $$

with

$$ \mathbf{G}_{j+\frac{1}{2}}^n = \mathbf{G}(\mathbf{U}_j^n, \mathbf{U}_{j+1}^n), \quad (3.2) $$

and

$$ \mathbf{G}(\mathbf{U}_L, \mathbf{U}_R) = \frac{1}{2} \left( \mathbf{G}(\mathbf{U}_L) + \mathbf{G}(\mathbf{U}_R) - \sum_{k=1}^{h} |\lambda_k|(\mathbf{U}_{k+1} - \mathbf{U}_k) \right), \quad (3.3) $$

where $\lambda_k$ the speeds of the discontinuities and $\mathbf{U}_k$ the intermediate states. In more details, we get that the natural discretization of the homogeneous version of the first two equations of system (2.16) is given by

$$ \left\{ \begin{aligned} \tau_j^{n+1^-} &= \tau_j^n + \frac{\Delta t}{\Delta m_j} (u_{j+\frac{1}{2}}^* - u_{j-\frac{1}{2}}^*), \\ u_j^{n+1^-} &= u_j^n - \frac{\Delta t}{\Delta m_j} (\Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^*), \end{aligned} \right. \quad (3.4) $$

with $\Delta m_j = \frac{\Delta x}{\tau_j^n}$ and $a_{j+\frac{1}{2}}^n = \max((Ac)_j^n, (Ac)_{j+1}^n)$. The numerical fluxes are given by $u_{j+\frac{1}{2}}^*$ and $\Pi_{j+\frac{1}{2}}^*$ at time $t^n$ and in particular, exploiting formulae (2.26),

$$ \begin{aligned} u_{j+\frac{1}{2}}^* &= u_{j+\frac{1}{2}}^* (\mathbf{U}_j^n, \mathbf{U}_{j+1}^n) = \frac{1}{2}(u_{j+1}^n + u_j^n) - \frac{1}{2a_{j+\frac{1}{2}}^n}(\Pi_{j+1}^n - \Pi_j^n), \\ \Pi_{j+\frac{1}{2}}^* &= \Pi_{j+\frac{1}{2}}^* (\mathbf{U}_j^n, \mathbf{U}_{j+1}^n) = \frac{1}{2}(\Pi_{j+1}^n + \Pi_j^n) - \frac{a_{j+\frac{1}{2}}^n}{2}(u_{j+1}^n - u_j^n). \end{aligned} \quad (3.5) $$

Note that the first equation $\tau_j^{n+1^-} = \tau_j^n + \frac{\Delta t}{\Delta m_j}(u_{j+\frac{1}{2}}^* - u_{j-\frac{1}{2}}^*)$ of system (3.4) is equivalent to

$$ A_j^n = A_j^{n+1-} \left( 1 + \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* - u_{j-\frac{1}{2}}^*) \right). \quad (3.6) $$

For more details we refer again to [5, 9, 10, 11, 18, 19, 25].
<figure><img src="image_2.png" /><figcaption>Figure 1: Connection between Lagrangian and Eulerian coordinates, see [25].</figcaption></figure>

It is useful to observe that a numerical discretization for the Lagrangian step starting from (2.12) with no source term, namely

$$ \left\{ \begin{aligned} L_j^{n+1-} A_j^{n+1-} &= L_j^n A_j^n \\ L_j^{n+1-} (Au)_j^{n+1-} &= L_j^n (Au)_j^n - \frac{\Delta t}{\Delta x} (\Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^*) \end{aligned} \right. \qquad (3.7) $$

with $\Pi_{j \pm \frac{1}{2}}^*$ given by (3.5), is equivalent to (3.4).

At last, observe from system (2.10) that a natural discretization for $L(\xi, t)$ is the following,

$$ L_j^{n+1-} = L_j^n + \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* - u_{j-\frac{1}{2}}^*) \quad \text{with} \quad L_j^n = 1, \qquad (3.8) $$

$u_{j \pm \frac{1}{2}}^*$ given by (3.5) and, from its definition (2.9),

$$ L_j(t) = \frac{1}{\Delta x} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} \frac{\partial x}{\partial \xi}(\xi, t) d\xi = \frac{x_{j+\frac{1}{2}}^* - x_{j-\frac{1}{2}}^*}{\Delta x} $$

with $x_{j \pm \frac{1}{2}}^* = x(\xi_{j \pm \frac{1}{2}}, t) \simeq \xi_{j \pm \frac{1}{2}} + t u_{j \pm \frac{1}{2}}^*$ according to (2.8), see also figure 1 for a better understanding.

*Projection step.* As already explained, the subsequent step consists in projecting the solution of system (2.12) (here with no source term), or equivalently the solution given by (3.7), in Eulerian coordinates. With this in mind, let us first note that by definition of $L$,

$$ \int_{x(\xi_l, t)}^{x(\xi_r, t)} \mathbf{Q}(x, t) dx = \int_{\xi_l}^{\xi_r} L(\xi, t) \bar{\mathbf{Q}}(\xi, t) d\xi. $$

This relation allows to make a simple link between $\mathbf{Q}$ in Eulerian and Lagrangian coordinates. In order to use it at the discrete level, it just remains to define $\hat{\xi}_{j+\frac{1}{2}}(t)$ such that for all $j$ (see again figure 1)

$$ x(\hat{\xi}_{j+\frac{1}{2}}(T), T) = x_{j+\frac{1}{2}}, \quad \text{with} \quad T \ge 0, $$

and the corresponding trajectories

$$ \left\{ \begin{aligned} \frac{\partial x}{\partial t}(\hat{\xi}_{j+\frac{1}{2}}(T), t) &= u(x(\hat{\xi}_{j+\frac{1}{2}}(T), t), t) \\ x(\hat{\xi}_{j+\frac{1}{2}}(T), 0) &= \hat{\xi}_{j+\frac{1}{2}}(T), \end{aligned} \right. $$

leading to

$$ \mathbf{Q}_j(t) = \frac{1}{\Delta x} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} \mathbf{Q}(x, t) dx = \frac{1}{\Delta x} \int_{x(\hat{\xi}_{j-\frac{1}{2}}, t)}^{x(\hat{\xi}_{j+\frac{1}{2}}, t)} \mathbf{Q}(x, t) dx = \frac{1}{\Delta x} \int_{\hat{\xi}_{j-\frac{1}{2}}}^{\hat{\xi}_{j+\frac{1}{2}}} L(\xi, t) \bar{\mathbf{Q}}(\xi, t) d\xi. \qquad (3.9) $$

Note that we can approximate $x_{j+\frac{1}{2}}$ at first-order:

$$ x_{j+\frac{1}{2}} = x(\hat{\xi}_{j+\frac{1}{2}}(T), T) \simeq x(\hat{\xi}_{j+\frac{1}{2}}(T), 0) + T \partial_t x(\hat{\xi}_{j+\frac{1}{2}}(T), 0) \simeq \hat{\xi}_{j+\frac{1}{2}} + T u_{j+\frac{1}{2}}^*, \qquad (3.10) $$
for a fixed time $T \ge 0$ and observe that $x_{j+\frac{1}{2}} - x_{j+\frac{1}{2}}^* = \hat{\xi}_{j+\frac{1}{2}} - \xi_{j+\frac{1}{2}}$. In order to define $\mathbb{Q}_j^{n+1}$ using the knowledge of $(L\mathbb{Q})_j^{n+1}$, we suggest to split the integral (3.9) to obtain

$$ \mathbb{Q}_j^{n+1} = \frac{1}{\Delta x} \int_{\hat{\xi}_{j-\frac{1}{2}}}^{\xi_{j-\frac{1}{2}}} L(\xi, t^{n+1-}) \bar{\mathbb{Q}}(\xi, t^{n+1-}) dx + \frac{1}{\Delta x} \int_{\xi_{j-\frac{1}{2}}}^{\xi_{j+\frac{1}{2}}} L(\xi, t^{n+1-}) \bar{\mathbb{Q}}(\xi, t^{n+1-}) dx + \frac{1}{\Delta x} \int_{\xi_{j+\frac{1}{2}}}^{\hat{\xi}_{j+\frac{1}{2}}} L(\xi, t^{n+1-}) \bar{\mathbb{Q}}(\xi, t^{n+1-}) dx. \quad (3.11) $$

Then, the middle integral can be clearly replaced by $\frac{1}{\Delta x} \int_{\xi_{j-\frac{1}{2}}}^{\xi_{j+\frac{1}{2}}} L(\xi, t) \bar{\mathbb{Q}}(\xi, t) dx = (L\bar{\mathbb{Q}})_j^{n+1}$, while the others two integrals can be approximated as in the following. Let us state $\frac{1}{\Delta x} \int_{\hat{\xi}_{j-\frac{1}{2}}}^{\xi_{j-\frac{1}{2}}} L(\xi, t) \bar{\mathbb{Q}}(\xi, t) dx = \frac{\xi_{j-\frac{1}{2}} - \hat{\xi}_{j-\frac{1}{2}}}{\Delta x} (L\bar{\mathbb{Q}})_{j-\frac{1}{2}}^{n+1}$, where we set

$$ (L\mathbb{Q})_{j-\frac{1}{2}}^{n+1-} = \begin{cases} (L\mathbb{Q})_{j-1}^{n+1-} & \text{if } \xi_{j-\frac{1}{2}} > \hat{\xi}_{j-\frac{1}{2}} \\ (L\mathbb{Q})_{j}^{n+1-} & \text{if } \xi_{j-\frac{1}{2}} \le \hat{\xi}_{j-\frac{1}{2}} \end{cases} $$

or equivalently

$$ (L\mathbb{Q})_{j-\frac{1}{2}}^{n+1-} = \begin{cases} (L\mathbb{Q})_{j-1}^{n+1-} & \text{if } x_{j-\frac{1}{2}}^{*,n+1-} > x_{j-\frac{1}{2}} \\ (L\mathbb{Q})_{j}^{n+1-} & \text{if } x_{j-\frac{1}{2}}^{*,n+1-} \le x_{j-\frac{1}{2}} \end{cases} \quad (3.12) $$

with

$$ x_{j+\frac{1}{2}}^{*,n+1-} = x_{j+\frac{1}{2}} + \Delta t u_{j+\frac{1}{2}}^* \quad \text{and} \quad x_{j+\frac{1}{2}} = \hat{\xi}_{j+\frac{1}{2}} + \Delta t u_{j+\frac{1}{2}}^* $$

and analogously for the third integral appearing in (3.11). Consequently, the projection step reads

$$ \begin{aligned} \mathbb{Q}_j^{n+1} &= \frac{\xi_{j-\frac{1}{2}} - \hat{\xi}_{j-\frac{1}{2}}}{\Delta x} (L\mathbb{Q})_{j-\frac{1}{2}}^{n+1-} + (L\mathbb{Q})_j^{n+1-} + \frac{\hat{\xi}_{j+\frac{1}{2}} - \xi_{j+\frac{1}{2}}}{\Delta x} (L\mathbb{Q})_{j+\frac{1}{2}}^{n+1-} \\ &= (L\mathbb{Q})_j^{n+1-} - \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* (L\mathbb{Q})_{j+\frac{1}{2}}^{n+1-} - u_{j-\frac{1}{2}}^* (L\mathbb{Q})_{j-\frac{1}{2}}^{n+1-}), \end{aligned} \quad (3.13) $$

with $u_{j+\frac{1}{2}}^*$ defined as in the Lagrangian step. Let us also remark that (3.13) can be seen as a discretization of system (2.14).

*Overall discretization.* It can be easily proved that the whole scheme is conservative using formulae (3.7) and (3.13). Indeed, one can recover the following final form:

$$ \begin{cases} A_j^{n+1} = A_j^n - \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* A_{j+\frac{1}{2}}^{n+1-} - u_{j-\frac{1}{2}}^* A_{j-\frac{1}{2}}^{n+1-}) \\ q_j^{n+1} = q_j^n - \frac{\Delta t}{\Delta x} (u_{j+\frac{1}{2}}^* q_{j+\frac{1}{2}}^{n+1-} + \Pi_{j+\frac{1}{2}}^* - (u_{j-\frac{1}{2}}^* q_{j-\frac{1}{2}}^{n+1-} + \Pi_{j-\frac{1}{2}}^*)), \end{cases} \quad (3.14) $$

with

$$ X_{j+\frac{1}{2}}^{n+1-} = \begin{cases} (LX)_j^{n+1-} & \text{if } x_{j+\frac{1}{2}}^{*,n+1-} > x_{j+\frac{1}{2}} \\ (LX)_{j+1}^{n+1-} & \text{if } x_{j+\frac{1}{2}}^{*,n+1-} \le x_{j+\frac{1}{2}}, \end{cases} $$

and $X$ either $A$ or $q$.

## 3.2 First-order approximation with varying parameters

As far as the source term is concerned, it is different from zero only when considering non-constant parameters, as the derivatives of $A_0$ and $K$ appears in it. Since the source term is taken into consideration at the level of the Lagrangian step, in order to include it into the numerical scheme we have to modify only the first step and not the remap one. Note however that, in general, the projection step has to be modified in order to obtain a fully well-balanced numerical method, see [5, 25]. We show that this is not required in our particular case as we want to preserve only the "man at eternal rest" solution. In the rest of this section, we present two different ways to obtain a well-balanced Lagrangian step in which the source term is included.
### 3.2.1 Based on the approximate Riemann solver

The first approach we present requires the use of the approximate Riemann solver described in section 2.2.2. In particular the well-balanced property is achieved approximating the source term in a special way, namely exploiting formula (2.33).

The first-order Godunov-type method associated with the Riemann solver proposed in section 2.2.2 now reads

$$ U_j^{n+1-} = U_j^n - \frac{\Delta t}{\Delta m_j} (G_{j+\frac{1}{2}}^n - G_{j-\frac{1}{2}}^n) + \Delta t \tilde{S}_j^n \quad (3.15) $$

with the numerical flux $G_{j+\frac{1}{2}}^n$ as in (3.2)-(3.3) and the source term defined as the average of the source at the interfaces,

$$ \tilde{S}_j^n = \frac{1}{2} \left( \frac{\Delta m_{j+\frac{1}{2}}}{\Delta m_j} \tilde{S}_{j+\frac{1}{2}}^n + \frac{\Delta m_{j-\frac{1}{2}}}{\Delta m_j} \tilde{S}_{j-\frac{1}{2}}^n \right), \quad \tilde{S}_{j+\frac{1}{2}}^n = \tilde{S}(\Delta m_{j+\frac{1}{2}}, \Delta t; U_j^n, U_{j+1}^n) \quad (3.16) $$

with $\Delta m_{j+1/2} = \frac{\Delta m_j + \Delta m_{j+1}}{2}$. For more details we refer again to Gallice [18, 19]. Hence, we approximate system (2.15) by

$$ \left\{ \begin{aligned} \tau_j^{n+1-} &= \tau_j^n + \frac{\Delta t}{\Delta m_j} (u_{j+\frac{1}{2}}^* - u_{j-\frac{1}{2}}^*) \\ u_j^{n+1-} &= u_j^n - \frac{\Delta t}{\Delta m_j} (\Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^*) + \Delta t \tilde{s}_j^n \end{aligned} \right. \quad (3.17) $$

or equivalently the Lagrangian system (2.12) by

$$ \left\{ \begin{aligned} L_j^{n+1-} A_j^{n+1-} &= L_j^n A_j^n \\ L_j^{n+1-} (Au)_j^{n+1-} &= L_j^n (Au)_j^n - \frac{\Delta t}{\Delta x} (\Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^*) + \Delta t s_j^n, \end{aligned} \right. \quad (3.18) $$

with $L_j^{n+1-}$ defined as in (3.8). Similarly, we observe that the numerical fluxes are now given by

$$ \left\{ \begin{aligned} u_{j+\frac{1}{2}}^* &= \frac{1}{2}(u_{j+1}^n + u_j^n) - \frac{1}{2a_{j+\frac{1}{2}}^n}(\Pi_{j+1}^n - \Pi_j^n) + \frac{\Delta m_{j+1/2}}{2a_{j+1/2}} \tilde{s}_{j+1/2}^n \\ \Pi_{j+\frac{1}{2}}^* &= \frac{1}{2}(\Pi_{j+1}^n + \Pi_j^n) - \frac{a_{j+\frac{1}{2}}^n}{2}(u_{j+1}^n - u_j^n). \end{aligned} \right. \quad (3.19) $$

As far as the source term (2.4) is concerned, we state

$$ \tilde{s}_j^n = \frac{1}{2} \left( \frac{\Delta m_{j+\frac{1}{2}}}{\Delta m_j} \tilde{s}_{j+1/2}^n + \frac{\Delta m_{j-\frac{1}{2}}}{\Delta m_j} \tilde{s}_{j-1/2}^n \right) \quad \text{with} \quad \tilde{s}_{j+1/2}^n = - \frac{M_{j+1/2}^n}{\Delta m_{j+1/2}}, \quad (3.20) $$

and $M_{j+1/2}^n = M((A_j^n; A_{0,j}, K_j); (A_{j+1}^n; A_{0,j+1}, K_{j+1}))$ given by (2.33). Note that, in (3.18), $s_j^n = L_j^n A_j^n \tilde{s}_j^n = A_j^n \tilde{s}_j^n = \frac{\Delta m_j^n}{\Delta x} \tilde{s}_j^n$ or equivalently

$$ s_j^n = \frac{1}{2} (s_{j+1/2}^n + s_{j-1/2}^n) \quad \text{with} \quad s_{j+1/2}^n = - \frac{M_{j+1/2}^n}{\Delta x} \quad \forall j. \quad (3.21) $$

**Theorem 1.** The Lagrangian step (3.18)-(3.21) (or equivalently (3.17)-(3.20)) here described is well-balanced under the "man at eternal rest" condition (2.7).

*Proof.* Let us suppose that $Q_j^n$ satisfies the "man at eternal rest" condition (2.7). Thus, $u_j^n = 0 \forall j$ and, thanks to condition (2.30), we also have that $\Pi_{j+1}^n - \Pi_j^n + M_{j+1/2} = 0$ and, consequently, $-(\Pi_{j+1}^n - \Pi_j^n) + \Delta m_{j+1/2} \tilde{s}_{j+1/2}^n = 0$. Hence, $u_{j+\frac{1}{2}}^* = 0 \forall j$ and $\Pi_{j+\frac{1}{2}}^* = \frac{1}{2}(\Pi_{j+1}^n + \Pi_j^n)$. Similarly $s_j^n$ compensates $-\frac{1}{\Delta x}(\Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^*)$ and we find that $A_j^{n+1-} = A_j^n$ and $q_j^{n+1-} = q_j^n$.

□
### 3.2.2 Based on the hydrostatic reconstruction

Let us consider once again a discretization of the form (3.18) for the Lagrangian system (2.12). Now we want to exploit the well-known hydrostatic reconstruction approach [1, 5], thus we are going to modify the numerical flux and source to be employed in either system (3.17) or (3.18). Hence, we perform a reconstruction of the variables values $Q$ at the cells interfaces. Considering locally relations (2.7), we first write for all $j$,

$$ \begin{cases} q(x) = 0, \\ (K(\sqrt{A} - \sqrt{A_0}))(x) = K_j(\sqrt{A_j^n} - \sqrt{A_{0,j}}) \end{cases} $$

which can be understood as a reconstruction procedure of a space dependent stationary solution inside the $j$-th cell. Then, following [20], we define cross-sectional areas at the cell interfaces ($x = x_{j \pm \frac{1}{2}}$) by

$$ \begin{cases} (K\sqrt{A})_{j+\frac{1}{2},L}^n = \max(K_j(\sqrt{A_j^n} - \sqrt{A_{0,j}}) + (K\sqrt{A_0})_{j+\frac{1}{2}}, 0) \\ (K\sqrt{A})_{j+\frac{1}{2},R}^n = \max(K_{j+1}(\sqrt{A_{j+1}^n} - \sqrt{A_{0,j+1}}) + (K\sqrt{A_0})_{j+\frac{1}{2}}, 0), \end{cases} \quad (3.22) $$

where the maximum has been added in order to preserve the non-negativity, and

$$ \left\{ \begin{aligned} \sqrt{A_{j+\frac{1}{2},L}^n} &= \max\left(\frac{1}{K_{j+\frac{1}{2}}^*}(K_j(\sqrt{A_j^n} - \sqrt{A_{0,j}}) + (K\sqrt{A_0})_{j+\frac{1}{2}}), 0\right) \\ \sqrt{A_{j+\frac{1}{2},R}^n} &= \max\left(\frac{1}{K_{j+\frac{1}{2}}^*}(K_{j+1}(\sqrt{A_{j+1}^n} - \sqrt{A_{0,j+1}}) + (K\sqrt{A_0})_{j+\frac{1}{2}}), 0\right) \end{aligned} \right. $$

where we have imposed

$$ (K\sqrt{A_0})_{j+\frac{1}{2}} = \min(K_j\sqrt{A_{0,j}}, K_{j+1}\sqrt{A_{0,j+1}}) $$

and

$$ K_{j+\frac{1}{2}}^* = \max(K_j, K_{j+1}). \quad (3.23) $$

Hence, the reconstructed values at cell interfaces are defined by

$$ \mathbf{Q}_{j+\frac{1}{2},L}^n = \begin{pmatrix} A_{j+\frac{1}{2},L}^n \\ A_{j+\frac{1}{2},L}^n u_i^n \end{pmatrix} \quad \text{and} \quad \mathbf{Q}_{j+\frac{1}{2},R}^n = \begin{pmatrix} A_{j+\frac{1}{2},R}^n \\ A_{j+\frac{1}{2},R}^n u_{i+1}^n \end{pmatrix}. \quad (3.24) $$

At this stage we can define the star values for the velocity $u$ and linearized pressure $\Pi$ as in (3.5) exploiting the values at interfaces (3.24), therefore

$$ u_{j+\frac{1}{2}}^* = u^*(\mathbf{Q}_{j+\frac{1}{2},L}^n, \mathbf{Q}_{j+\frac{1}{2},R}^n) \quad \text{and} \quad \Pi_{j+\frac{1}{2}}^* = \Pi^*(\mathbf{Q}_{j+\frac{1}{2},L}^n, \mathbf{Q}_{j+\frac{1}{2},R}^n) $$

with

$$ \left\{ \begin{aligned} u_{j+\frac{1}{2}}^* &= \frac{1}{2}(u_{j+1}^n + u_j^n) - \frac{1}{2a_{j+\frac{1}{2}}^n}(\Pi_{j+\frac{1}{2},R}^n - \Pi_{j+\frac{1}{2},L}^n) \\ \Pi_{j+\frac{1}{2}}^* &= \frac{1}{2}(\Pi_{j+\frac{1}{2},R}^n + \Pi_{j+\frac{1}{2},L}^n) - \frac{a_{j+\frac{1}{2}}^n}{2}(u_{j+1}^n - u_j^n) \end{aligned} \right. \quad (3.25) $$

and

$$ a_{j+\frac{1}{2}} = \max(A_{j+\frac{1}{2},L}^n c_{j+\frac{1}{2},L}^n, A_{j+\frac{1}{2},R}^n c_{j+\frac{1}{2},R}^n). $$

In particular, (3.25) has to be used in (3.17)-(3.18) instead of (3.19). Finally, since we want to preserve the steady states with zero velocity, namely the ones satisfying $u = 0$ and $\partial_x \tilde{p} = s$, for the source term $s_j^n$ we suggest the following,

$$ s_j^n = \frac{1}{\Delta x} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} s(\mathbf{Q}; A_0, K)(x, t) dx = \frac{1}{\Delta x} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} \partial_x \tilde{p}(x, t) dx = \frac{1}{\Delta x} (\tilde{p}(x_{j+\frac{1}{2}}) - \tilde{p}(x_{j-\frac{1}{2}})) = \frac{\tilde{p}_{j+\frac{1}{2},L} - \tilde{p}_{j-\frac{1}{2},R}}{\Delta x}. \quad (3.26) $$

Note that the first equality comes from space dependent reconstruction of a solution inside the cell, while the second equality holds as we are exploiting the reconstructed values (3.24) to define the pressures $\tilde{p}$ at the interfaces. More precisely, let us
remark that to define $\tilde{p}_{j+\frac{1}{2},L}^n = \tilde{p}(A_{j+\frac{1}{2},L}^n, K)$, $\tilde{p}_{j+\frac{1}{2},R}^n = \tilde{p}(A_{j+\frac{1}{2},R}^n, K)$ in (3.26), and $\Pi_{j+\frac{1}{2},L}^n, \Pi_{j+\frac{1}{2},R}^n$ at equilibrium in (3.25), namely $\Pi_{j+\frac{1}{2},L}^n = \tilde{p}(A_{j+\frac{1}{2},L}^n, K)$, $\Pi_{j+\frac{1}{2},R}^n = \tilde{p}(A_{j+\frac{1}{2},R}^n, K)$, we can use either the values $K = K(x_{j+\frac{1}{2}})$ or $K = K_{j+\frac{1}{2}}^*$ as in (3.23) with the only request to define the source term (3.26) accordingly.

**Theorem 2.** The Lagrangian step (3.18), (3.25)-(3.26) is well-balanced under the "man at eternal rest" condition (2.7).

*Proof.* Let us suppose that $Q_j^n$ satisfies the "man at eternal rest" condition (2.7). Thus, $u_j^n = 0 \forall j$, $A_{j+\frac{1}{2},L}^n = A_{j+\frac{1}{2},R}^n$ and consequently, $Q_{j+\frac{1}{2},L}^n = Q_{j+\frac{1}{2},R}^n$. Hence, $u_{j+\frac{1}{2}}^* = 0 \forall j$ and $-\frac{1}{\Delta x}(\Pi_{j+\frac{1}{2}}^* - \Pi_{j-\frac{1}{2}}^*) + \{s\}_j^n = 0$. Therefore, we obtain $A_j^{n+1-} = A_j^n$ and $q_j^{n+1-} = q_j^n$. $\square$

*Remark.* As it has been shown in the proofs of theorems 1 and 2, $u_{j+\frac{1}{2}}^* = 0$ under the "man at eternal rest" condition. Hence, the projection step (3.13) preserves the stationary solution (2.7) and the whole numerical scheme (Lagrangian plus projection step) has not to be further modified.

# 4 Second-order well-balanced scheme

At this stage we are interested in second (or higher) order extension of the Lagrangian-projection schemes. We proceed as for the first-order scheme: we explain how to reach the second-order of accuracy in the case of constant parameters, which is already non-trivial due to the presence of two steps in the Lagrange-Projection procedure, and then we extend it to the case of varying parameters $A_0$ and $K$. Once again, in the latter case we pay attention to the well-balanced property.

Here we focus on a second-order simplified version of the scheme applied to (2.3) but it can be easily extended to higher order of accuracy following [25, 6]. In particular, we make use of polynomial reconstruction and Runge-Kutta TVD scheme [23] in order to reach second-order of accuracy respectively in space and time.

## 4.1 Second-order approximation with constant parameters

First of all we explain how to reach the second-order of accuracy in space in both the Lagrangian and projection step, and then we comment on the Runge-Kutta TVD scheme for the second-order approximation in time. Hereafter the time is assumed to be left continuous for the sake of clarity.

*Space discretization of the Lagrangian step.* Given a time $t$, a $j$-th cell and the cell value $Q_j(t)$, this step aims at defining evolved values $Q_{j+\frac{1}{2},L}(t)$ at the cell interface $x_{j+\frac{1}{2}}$ by means of polynomial data reconstructions. More precisely, using for each cell $I_j$ a reconstructed polynomial vector $\mathbf{P}_j^t(x)$ such as

$$ \mathbf{P}_j^t(x) = \mathbf{Q}_j(t) + \Delta_j^t(x - x_j), \quad (4.1) $$

where $\Delta_j^t = \Delta_j^t(\mathbf{Q}_{j-1}(t), \mathbf{Q}_j(t), \mathbf{Q}_{j+1}(t))$ is the slope, either the ENO [34] or the MINMOD [30] one, we define

$$ \mathbf{Q}_{j+\frac{1}{2}L}(t) = \mathbf{P}_j^t(x_{j+\frac{1}{2}}) \quad \text{and} \quad \mathbf{Q}_{j+\frac{1}{2}R}(t) = \mathbf{P}_{j+1}^t(x_{j+\frac{1}{2}}). \quad (4.2) $$

Therefore, once again we use formulae (3.5) computed in $(\mathbf{Q}_{j+\frac{1}{2}L}(t), \mathbf{Q}_{j+\frac{1}{2}R}(t))$, namely

$$ u_{j+\frac{1}{2}}^*(t) = u_{j+\frac{1}{2}}^*(\mathbf{Q}_{j+\frac{1}{2}L}(t), \mathbf{Q}_{j+\frac{1}{2}R}(t)) \quad \text{and} \quad \Pi_{j+\frac{1}{2}}^*(t) = \Pi_{j+\frac{1}{2}}^*(\mathbf{Q}_{j+\frac{1}{2}L}(t), \mathbf{Q}_{j+\frac{1}{2}R}(t)), $$

with

$$ a_{j+\frac{1}{2}} = \max(A_{j+\frac{1}{2},L}^n c_{j+\frac{1}{2},L}^n, A_{j+\frac{1}{2},R}^n c_{j+\frac{1}{2},R}^n). $$

Note that the polynomial $\mathbf{P}_j^t(x)$ should satisfy the conservation property, which reads

$$ \frac{1}{\Delta x} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} \mathbf{P}_j^t(x) = \mathbf{Q}_j(t). $$
Space discretization of the remap step. In order to obtain the second-order of accuracy in space we exploit relations (3.11)-(3.13) seen in section 3.1. Once again the middle integral in (3.11) can be substituted by $(L\mathbf{Q})_j(t)$ thanks to the conservation property. Then, for the other two integrals, instead of considering the values $(L\mathbf{Q})_j(t)$, we reconstruct them using the polynomial $\mathbf{P}_j^t(x)$. Thus, we introduce

$$ (L\mathbf{P})_j^t(\xi) = (L\mathbf{Q})_j(t) + \Delta_j^t(\xi - \xi_j), \quad (4.3) $$

with $\Delta_j^t = \Delta_j^t((L\mathbf{Q})_{j-1}(t), (L\mathbf{Q})_j(t), (L\mathbf{Q})_{j+1}(t))$ and

$$ (L\mathbf{P})_{j-\frac{1}{2}}^t(\xi) = \begin{cases} (L\mathbf{P})_{j-1}^t(\xi) & \text{if } \xi_{j-\frac{1}{2}} > \hat{\xi}_{j-\frac{1}{2}} \\ (L\mathbf{P})_j^t(\xi) & \text{if } \xi_{j-\frac{1}{2}} \le \hat{\xi}_{j-\frac{1}{2}} \end{cases} $$

in place of $(L\mathbf{Q})_{j-\frac{1}{2}}(t)$, where we recall that $x_{j+\frac{1}{2}} = \hat{\xi}_{j+\frac{1}{2}} + \Delta t u_{j+\frac{1}{2}}^*$. Finally, since $(L\mathbf{P})^t(\xi)$ is not constant, in order to be able to evaluate its integral $\frac{1}{\Delta x} \int_{\hat{\xi}_{j-\frac{1}{2}}}^{\xi_{j-\frac{1}{2}}} (L\mathbf{P})_{j-\frac{1}{2}}^t(\xi) dx$, one can exploit either the mid-point rule (only for second-order of accuracy) or a Gauss quadrature formula with nodes $\xi_{j-\frac{1}{2},k}$ and weights $\omega_k$ for $k = 1, \dots, m$. In the former case we find

$$ \begin{aligned} \mathbf{Q}_j(t) &= (L\mathbf{Q})_j(t) + \frac{\xi_{j-\frac{1}{2}} - \hat{\xi}_{j-\frac{1}{2}}}{\Delta x} (L\mathbf{P})_{j-\frac{1}{2}}^t \left( \frac{\xi_{j-\frac{1}{2}} + \hat{\xi}_{j-\frac{1}{2}}}{2} \right) + \frac{\hat{\xi}_{j+\frac{1}{2}} - \xi_{j+\frac{1}{2}}}{\Delta x} (L\mathbf{P})_{j+\frac{1}{2}}^t \left( \frac{\xi_{j+\frac{1}{2}} + \hat{\xi}_{j+\frac{1}{2}}}{2} \right) \\ &= (L\mathbf{Q})_j(t) - \frac{\Delta t}{\Delta x} \left( u_{j+\frac{1}{2}}^* (L\mathbf{P})_{j+\frac{1}{2}}^t \left( \frac{\xi_{j+\frac{1}{2}} + \hat{\xi}_{j+\frac{1}{2}}}{2} \right) - u_{j-\frac{1}{2}}^* (L\mathbf{P})_{j-\frac{1}{2}}^t \left( \frac{\xi_{j-\frac{1}{2}} + \hat{\xi}_{j-\frac{1}{2}}}{2} \right) \right), \end{aligned} \quad (4.4) $$

while in the second one,

$$ \begin{aligned} \mathbf{Q}_j(t) &= (L\mathbf{Q})_j(t) + \frac{\xi_{j-\frac{1}{2}} - \hat{\xi}_{j-\frac{1}{2}}}{\Delta x} \sum_{k=1}^{m} \omega_k (L\mathbf{P})_{j-\frac{1}{2}}^t (\xi_{j-\frac{1}{2},k}) + \frac{\hat{\xi}_{j+\frac{1}{2}} - \xi_{j+\frac{1}{2}}}{\Delta x} \sum_{k=1}^{m} \omega_k (L\mathbf{P})_{j+\frac{1}{2}}^t (\xi_{j+\frac{1}{2},k}) \\ &= (L\mathbf{Q})_j(t) - \frac{\Delta t}{\Delta x} \left( u_{j+\frac{1}{2}}^* \sum_{k=1}^{m} \omega_k (L\mathbf{P})_{j+\frac{1}{2}}^t (\xi_{j+\frac{1}{2},k}) - u_{j-\frac{1}{2}}^* \sum_{k=1}^{m} \omega_k (L\mathbf{P})_{j-\frac{1}{2}}^t (\xi_{j-\frac{1}{2},k}) \right), \end{aligned} $$

where the other integral $\frac{1}{\Delta x} \int_{\hat{\xi}_{j+\frac{1}{2}}}^{\xi_{j+\frac{1}{2}}} (L\mathbf{P})_{j+\frac{1}{2}}^t(\xi) dx$ has been evaluated in a similar way.

Second-order approximation in time. The last step of the second-order method consists in obtaining the right accuracy in time. In order to do this we simply use the Runge-Kutta TVD scheme at second order [23]. However, we have to specify that it has to be applied to the overall scheme (Lagrangian and remap step together) in order to avoid diffusion due to the splitting.

## 4.2 Second-order approximation with varying parameters

As we have done before for the well-balanced first-order scheme, here we present two different well-balanced second-order methods, the first exploits the approximate Riemann solver of section 2.2 and the other one the hydrostatic reconstruction approach.

Again, it is sufficient to focus on the Lagrangian step and nothing has to be changed for the second-order projection step as described above and the Runge-Kutta TVD procedure as it will be easily seen that they preserve the "man at eternal rest" solution. However, regarding the projection step, we specify that, if we would like to preserve a different stationary solution, we would have to modify it. For more details refer to [11].

### 4.2.1 Based on the approximate Riemann solver

Here we describe the second-order extension of the Lagrangian step explained in section 3.2.1 which makes use of the approximate Riemann solver of section 2.2.2 in order to maintain the well-balancedness of the method.
Thus, in order to obtain a second-order approximation in space, the idea is to exploit the reconstructed values at cell interfaces and then apply the usual updating formulae. However, we cannot simply use the reconstructed polynomial (4.1) as it would prevent the scheme to preserve the stationary solutions. Thus, the idea is to compute the slopes in such a way that they become equal to zero when the "man at eternal rest" condition (2.7) is satisfied. With this in mind we suggest to make use of the so-called fluctuations $\mathbf{D}$, refer to [25, 6]. Let us enter the details. Given a time $t$ and the $j$-th cell, to determine the slope at second-order of accuracy we need the values $\mathbf{Q}_{j-1}(t), \mathbf{Q}_j(t), \mathbf{Q}_{j+1}(t)$; thus for $k = j - 1, j, j + 1$, we define the so-called fluctuations

$$ \mathbf{D}_{k,j}(t) = \mathbf{Q}_k(t) - \frac{1}{\Delta x} \int_{x_{k-\frac{1}{2}}}^{x_{k+\frac{1}{2}}} \mathbf{Q}_j^{t,e}(x) dx, \quad (4.5) $$

where $\mathbf{Q}_j^{t,e}(x)$ denotes a reconstructed stationary solution we want to preserve and which satisfies

$$ \frac{1}{\Delta x} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} \mathbf{Q}_j^{t,e}(x) dx = \mathbf{Q}_j(t). \quad (4.6) $$

Note that, since we are interested in a second-order accurate scheme, we solve the integral in (4.5) using the mid-point rule in space. Usually it is not straightforward to compute $\mathbf{Q}_j^{t,e}(x)$ with the constraint (4.6), however, since we only want to preserve the "man at eternal rest" solution, we can automatically define it such that

$$ (K\sqrt{A_j^{t,e}})(x) = K_j(\sqrt{A_j} - \sqrt{A_{0,j}}) + (K\sqrt{A_0})(x) \quad \text{and} \quad u_j^{t,e}(x) = u_j^t. \quad (4.7) $$

Consequently we denote the reconstruction operator as

$$ \mathbf{P}_j^t(x) = \mathbf{P}_j^t(x; \mathbf{Q}_j(t), \mathbf{D}_{j-1,j}(t), \mathbf{D}_{j,j}(t), \mathbf{D}_{j+1,j}(t)) = \mathbf{Q}_j(t) + \Delta_j^t(x - x_j) \quad (4.8) $$

where, for the sake of clarity, we specify that $\Delta_j^t = \Delta_j^t(\mathbf{D}_{j-1,j}(t), \mathbf{D}_{j,j}(t), \mathbf{D}_{j+1,j}(t))$, for whose definition we will use either the ENO or the MINMOD slope. Let us observe that in our specific case $\mathbf{D}_{j,j}(t) = 0$ always and that $\mathbf{D}_{j\pm 1,j}(t) = 0$ when $\mathbf{Q}_j(t)$ satisfies a stationary solution. Therefore it is clear that the slopes equal zero when the "man at eternal rest" condition is satisfied.

Equipped with the definition for the slope and the same definition (4.2) for the left and right traces, we can now compute $u_{j+\frac{1}{2}}^*$ and $\Pi_{j+\frac{1}{2}}^*$ as in (3.19), namely

$$ u_{j+\frac{1}{2}}^*(t) = u^*(\mathbf{Q}_{j+\frac{1}{2},L}(t), \mathbf{Q}_{j+\frac{1}{2},R}(t)) \quad \text{and} \quad \Pi_{j+\frac{1}{2}}^*(t) = \Pi^*(\mathbf{Q}_{j+\frac{1}{2},L}(t), \mathbf{Q}_{j+\frac{1}{2},R}(t)) \quad (4.9) $$

with

$$ \mathbf{Q}_{j+\frac{1}{2}L}(t) = \mathbf{P}_j^t(x_{j+\frac{1}{2}}) \quad \text{and} \quad \mathbf{Q}_{j+\frac{1}{2}R}(t) = \mathbf{P}_{j+1}^t(x_{j+\frac{1}{2}}). $$

More specifically, we state

$$ \left\{ \begin{aligned} u_{j+\frac{1}{2}}^*(t) &= \frac{1}{2}(u_{j+\frac{1}{2}R}^t + u_{j+\frac{1}{2}L}^t) - \frac{1}{2a_{j+\frac{1}{2}}^t}(\Pi_{j+\frac{1}{2}R}^t - \Pi_{j+\frac{1}{2}L}^t) + \frac{\Delta m_{j+1/2}^t}{2a_{j+1/2}^t}s_{j+1/2}^t \\ \Pi_{j+\frac{1}{2}}^*(t) &= \frac{1}{2}(\Pi_{j+\frac{1}{2}R}^t + \Pi_{j+\frac{1}{2}L}^t) - \frac{a_{j+\frac{1}{2}}^t}{2}(u_{j+\frac{1}{2}R}^t - u_{j+\frac{1}{2}L}^t). \end{aligned} \right. $$

At last, one needs to specify the value of $K$ used in the definition $\Pi_{j+\frac{1}{2}R,L}^t$. It turns out that using the natural value of $K$ in $x_{j+\frac{1}{2}}$ leads to the loss of the well-balanced property. Therefore we propose to reconstruct $K$ exploiting the usual reconstruction polynomial and corresponding slope $\Delta_j^t$. However, since the equilibrium part of $K$ reads $K_j^{t,e}(x) = K(x)$ as $K$ is known and does not depend on time, its fluctuations are null and, thus, we can simply state that

$$ K_{j+\frac{1}{2}L} = K_j \quad \text{and} \quad K_{j+\frac{1}{2}R} = K_{j+1}, $$

which is very convenient as with this choice of $K$ we do not have to do further modifications to preserve the well-balancedness of the scheme. At last, to second-order of accuracy, for the source term we simply consider (3.20)-(3.21) with $M_{j+1/2}^t = M((A_{j+\frac{1}{2},L}(t); A_{0,j}, K_j); (A_{j+\frac{1}{2},R}(t); A_{0,j+1}, K_{j+1}))$ defined as in (2.33), and the updating formulae (3.18).
**Theorem 3.** The Lagrangian step (4.5)-(4.9) with updating formulae (3.18) here described is well-balanced under the "man at eternal rest" condition (2.7).

*Proof.* Once again, let us suppose that $Q_j^n$ satisfies the "man at eternal rest" condition (2.7) and, as such, $u_j^n = 0 \forall j$. Due to definitions (4.5)-(4.7) for the fluctuations, we have that the slope satisfies $\Delta_j^t = 0$ and, thus, as for the first-order scheme, we have $u_{j+\frac{1}{2}}^* = 0$ and $A_j^{n+1-} = A_j^n$. As for the proof of the corresponding first-order scheme, simple algebraic computations show that $q_j^{n+1-} = q_j^n$. □

Lastly, for the benefit of the reader, in the following we insert a summary of this second order scheme.

<table><thead><tr><th>Algorithm 1.</th><th>Second order scheme based on the approximate Riemann solver.</th></tr></thead><tbody><tr><td>Step 1.</td><td>Start with the Lagrangian step; look for the reconstructed stationary solution Q<sub>j</sub><sup>t,e</sup>(x) that satisfies (4.6).</td></tr><tr><td>Step 2.</td><td>Exploit it to compute the fluctuations D<sub>k,j</sub>(t) as in (4.5).</td></tr><tr><td>Step 3.</td><td>Define the reconstruction operator P<sub>j</sub><sup>t</sup>(x) as in (4.8).</td></tr><tr><td>Step 4.</td><td>Find the reconstructed values Q<sub>j+1/2,L</sub>(t) = P<sub>j</sub><sup>t</sup>(x<sub>j+1/2</sub>), Q<sub>j+1/2,R</sub>(t) = P<sub>j+1</sub><sup>t</sup>(x<sub>j+1/2</sub>).</td></tr><tr><td>Step 5.</td><td>Compute u<sub>j+1/2</sub><sup>*</sup>, Π<sub>j+1/2</sub><sup>*</sup> as in (3.19)-(4.9) and the source term as in (3.20)-(3.21) exploiting the reconstructed values.</td></tr><tr><td>Step 6.</td><td>Solve system (3.18) written in Lagrangian coordinates obtaining (LQ)(t).</td></tr><tr><td>Step 7.</td><td>Continue with the remap step; define the polynomial P<sub>j</sub><sup>t</sup>(x) as in (4.3).</td></tr><tr><td>Step 8.</td><td>Update Q(t) using formula (4.4).</td></tr><tr><td>Step 9.</td><td>Apply the Runge-Kutta scheme in order to reach the second order of accuracy in time.</td></tr></tbody></table>

### 4.2.2 Based on the hydrostatic reconstruction

Let us now see how to modify the second-order accurate scheme exploiting the hydrostatic reconstruction, already introduced in section 3.2.2. Once again, in order to have a well-balanced scheme, we only have to modify the Lagrangian step.

Hence, in this case, we would like to combine two different kinds of reconstruction, one for the well-balancedness and one for the high-order accuracy. Thus, at the end we will have a unique reconstruction function for $Q_j^t(x)$ which will consist of two parts: the fluctuation $P_j^t(x)$ and the equilibrium $Q_j^{t,e}(x)$ one,

$$Q_j^t(x) = Q_j^{t,e}(x) + P_j^t(x), \qquad (4.10)$$

as suggested in [25]. Subsequently, we will use the values $Q_{j+\frac{1}{2}L}^t = Q_j^t(x_{j+\frac{1}{2}})$ and $Q_{j+\frac{1}{2}R}^t = Q_{j+1}^t(x_{j+\frac{1}{2}})$ to compute $u_{j+\frac{1}{2}}^*$ and $\Pi_{j+\frac{1}{2}}^*$ in (3.5) and to find $(LQ)_j(t)$ according to (3.18).

At this stage we have to define $P_j^t(x)$ and $Q_j^{t,e}(x)$. For the latter, since we are interested in preserving only the "man at eternal rest solution", we simply use the well-balanced reconstructed values (3.24), in particular

$$Q_j^{t,e}(x_{j+\frac{1}{2}}) = Q_{j+\frac{1}{2},L}^{t,e} = \begin{pmatrix} A_{j+\frac{1}{2},L}^{t,e} \\ A_{j+\frac{1}{2},L}^{t,e} u_i^t \end{pmatrix} \quad \text{and} \quad Q_{j+1,x_{j+\frac{1}{2}}}^{t,e} = Q_{j+\frac{1}{2},R}^{t,e} = \begin{pmatrix} A_{j+\frac{1}{2},R}^{t,e} \\ A_{j+\frac{1}{2},R}^{t,e} u_i^t \end{pmatrix},$$

with $A_{j+\frac{1}{2},L}^{t,e}$, $A_{j+\frac{1}{2},R}^{t,e}$ computed as in (3.22)-(3.23). See either section 3.2.2 or [20] for more details. As far as the reconstruction polynomial $P_j^t(x)$ is concerned, we use a similar but not equal strategy to the one explained in the previous section. Indeed, here we write $P_j^t(x)$ only depending on the fluctuations $D_{j,j}(t)$, namely

$$P_j^t(x) = P_j^t(x; D_{j-1,j}(t), D_{j,j}(t), D_{j+1,j}(t)) = D_{j,j}(t) + \Delta_j^t(x - x_j), \qquad (4.11)$$

with $\Delta_j^t = \Delta_j^t(D_{j-1,j}(t), D_{j,j}(t), D_{j+1,j}(t))$ and $D_{k,j}(t)$, with $k = j - 1, j, j + 1$, defined as in (4.5). Let us remark that here $D_{j,j}(t) = 0$ always and, thus, $Q_j^t(x)$ reads

$$Q_j^t(x) = Q_j^{t,e}(x) + \Delta_j^t(x - x_j) = Q_j^{t,e}(x) + Q_j^{t,f}(x), \qquad (4.12)$$
where we have renamed $Q_j^{t,f}(x) = \Delta_j^t(x - x_j)$ for the fluctuations part. Then, let us compare (4.8) and (4.12). Indeed the term $Q_j(t)$ that appears in the right hand side of (4.8) can be understood as $Q_j^{t,e}(x)$ in (4.12) but replaced by $Q_j(t)$ using (4.6) and the mid-point rule.

Finally, we only need to specify the definition of the source term, which in general is defined by

$$s_j^t = \frac{1}{\Delta x} \int_{\xi_{j-\frac{1}{2}}}^{\xi_{j+\frac{1}{2}}} s(\xi, t) d\xi$$

with $s(A; A_0, K)$ given by (2.4). Since we aim to reach the second order of accuracy, in this case the mid-point rule in space suffices. Thus, using the equilibrium and fluctuation decomposition (4.12) for the cross-sectional area $A$, we can write

$$ \begin{aligned} s_j^t &= \frac{1}{\Delta x} \int_{\xi_{j-\frac{1}{2}}}^{\xi_{j+\frac{1}{2}}} s(A; A_0, K)(\xi, t) d\xi = \frac{1}{\Delta x} \int_{\xi_{j-\frac{1}{2}}}^{\xi_{j+\frac{1}{2}}} s(A^e + A^f; A_0, K)(\xi, t) d\xi \\ &= \frac{1}{\Delta x} \int_{\xi_{j-\frac{1}{2}}}^{\xi_{j+\frac{1}{2}}} (s(A^e + A^f; A_0, K)(\xi, t) - s(A^e; A_0, K)(\xi, t)) d\xi + \frac{1}{\Delta x} \int_{\xi_{j-\frac{1}{2}}}^{\xi_{j+\frac{1}{2}}} s(A^e; A_0, K)(\xi, t) d\xi \\ &= \frac{1}{\Delta x} \int_{\xi_{j-\frac{1}{2}}}^{\xi_{j+\frac{1}{2}}} s(A^e; A_0, K)(\xi, t) d\xi, \end{aligned} $$

where the last quality holds as, when applying the mid-point rule, the fluctuations part $A^f$ disappears leaving only $\frac{\xi_{j+\frac{1}{2}} - \xi_{j-\frac{1}{2}}}{\Delta x} (s(A^e + A^f; A_0, K)(x_j, t) - s(A^e; A_0, K)(x_j, t)) = \frac{\xi_{j+\frac{1}{2}} - \xi_{j-\frac{1}{2}}}{\Delta x} (s(A^e + 0; A_0, K)(x_j, t) - s(A^e; A_0, K)(x_j, t)) = 0$. Hence, similarly to the first-order scheme, the source term finally reads

$$s_j^t = \frac{1}{\Delta x} \int_{\xi_{j-\frac{1}{2}}}^{\xi_{j+\frac{1}{2}}} s(A^e; A_0, K)(\xi, t) d\xi = \frac{\tilde{p}_{j+\frac{1}{2},L}^e(t) - \tilde{p}_{j-\frac{1}{2},R}^e(t)}{\Delta x}. \quad (4.13)$$

**Theorem 4.** The Lagrangian step (4.10)-(4.13) is well-balanced under the "man at eternal rest" condition (2.7).

*Proof.* It is straightforward to see that, under the "man at eternal rest" condition (2.7), the fluctuations part in (4.12) satisfies $Q_j^{t,f}(x) = 0$. Consequently, the Lagrangian step is reduced to the first-order one, which we already proved to be well-balanced. □

As we have done in the previous section, for the sake of clarity we insert an algorithm to summarize this second order method. Let us note that the projection and Runge-Kutta steps are the same of the previous scheme.

<table><thead><tr><th colspan="2">Algorithm 2. Second order scheme based on the hydrostatic reconstruction.</th></tr></thead><tbody><tr><td>Step 1.</td><td>Start with the Lagrangian step; compute the well-balanced reconstructed values Q<sub>j</sub><sup>t,e</sup>(x) as in (3.24).</td></tr><tr><td>Step 2.</td><td>Find the fluctuations D<sub>k,j</sub>(t) by (4.5).</td></tr><tr><td>Step 3.</td><td>Define the reconstruction polynomial P<sub>j</sub><sup>t</sup>(x) only depending on the fluctuations as in (4.11).</td></tr><tr><td>Step 4.</td><td>Use Q<sub>j</sub><sup>t,e</sup>(x) and P<sub>j</sub><sup>t</sup>(x) to determine the reconstruction function Q<sub>j</sub><sup>t</sup>(x) by (4.10).</td></tr><tr><td>Step 5.</td><td>Find the reconstructed values Q<sub>j+1/2</sub><sup>t</sup> = Q<sub>j</sub><sup>t</sup>(x<sub>j+1/2</sub>) and Q<sub>j+1/2</sub><sup>t</sup> = Q<sub>j+1</sub><sup>t</sup>(x<sub>j+1/2</sub>).</td></tr><tr><td>Step 6.</td><td>Compute u<sub>j+1/2</sub><sup>*</sup> and Π<sub>j+1/2</sub><sup>*</sup> as in (3.5) exploiting the reconstructed values.</td></tr><tr><td>Step 7.</td><td>Define the source term s<sub>j</sub><sup>t</sup> according to (4.13).</td></tr><tr><td>Step 8.</td><td>Solve system (3.18) written in Lagrangian coordinates obtaining (LQ)(t).</td></tr><tr><td>Step 9.</td><td>Continue with the remap step; define the polynomial P<sub>j</sub><sup>t</sup>(x) as in (4.3).</td></tr><tr><td>Step 10.</td><td>Update Q(t) using formula (4.4).</td></tr><tr><td>Step 11.</td><td>Apply the Runge-Kutta scheme in order to reach the second order of accuracy in time.</td></tr></tbody></table>
# 5 Numerical simulations

In this section we carry out different numerical simulations that aim to show the good behaviour of the proposed numerical schemes. First of all, we numerically prove that, both in the case of system of conservation and balance laws, the numerical schemes reach the required order of convergence. A Riemann problem is also presented in the case of constant parameters. Then, different test cases are introduced in order to assess the well-balancedness and the wave-capturing properties of the numerical methods when applied to the non-conservative system.

*Time step and CFL condition.* Since the Lagrange-projection approach leads to a splitting of the original system (2.3) into the acoustic (2.13) and the transport (2.14) one, the time step is computed as the minimum between the two time steps obtained from (2.13) and (2.14). As far as the Lagrangian system is concerned, the Courant-Friedrichs-Lewy (CFL) condition reads

$$ \Delta t \leq \text{CFL}_t \frac{\Delta x}{\max_j \{ \max(\tau_j^n, \tau_{j+1}^n) a_{j+\frac{1}{2}} \}}, \quad (5.1) $$

while, for the transport system

$$ \Delta t \leq \text{CFL}_t \frac{\Delta x}{\max_j \{ u_{j-\frac{1}{2}}^{+} - u_{j+\frac{1}{2}}^{-} \}}, \quad (5.2) $$

where CFL$_t$ and CFL$_t$ are respectively the CFL number for the Lagrangian and the transport system, and finally

$$ u_{j-\frac{1}{2}}^{+} = \max(u_{j-\frac{1}{2}}^{*}, 0) \quad \text{and} \quad u_{j+\frac{1}{2}}^{-} = \min(u_{j+\frac{1}{2}}^{*}, 0). $$

When considering a first-order scheme, we can take $CFL_t \leq 0.5$ and $CFL_t < 1$. For more details refer to [10, 11, 12].

*Remark.* It is not difficult to prove that the first-order approximation (3.14) we presented preserves the strict positivity of the cross-sectional area A under the CFL conditions (5.1)-(5.2) with $CFL_t \leq \frac{1}{2}$ and $CFL_t < 1$, see also [5]. This statement remains true even if the parameters K and $A_0$ are not constant in space.

## 5.1 System of conservation laws

*Convergence study.*

In this section we start assessing that the second-order scheme of section 4.1 reaches the right order of accuracy. For this purpose, we need to compare the numerical solution with the exact one, specifying that, in order to obtain the correct order of accuracy, the smoothness of the exact solution is required. Since in general it is not known for system (2.5), we have to exploit the method of the manufactured solution, for which we refer to [28]. In a nutshell, given an acceptable smooth function $\hat{Q}$, we have to modify the sought system in such a way that $\hat{Q}$ is actually one of its solutions. This is achieved by adding a source term to the starting system, namely passing from the homogeneous version of system (2.5), $\partial_t \hat{Q} + \partial_x \hat{F}(\hat{Q}) = 0$, to

$$ \partial_t \hat{Q} + \partial_x \hat{F}(\hat{Q}) = \hat{S}(\hat{Q}). $$

$\hat{S}(\hat{Q})$ is usually found through an algebraic manipulator, and thus we will not report here the modified source term.

In particular, referring to [27], we have considered the following solution,

$$ \hat{Q} = \begin{pmatrix} \hat{A} \\ \hat{q} \end{pmatrix} = \begin{pmatrix} \tilde{A} + \tilde{a} \sin(\frac{2\pi}{L}x) \cos(\frac{2\pi}{T_0}t) \\ \tilde{q} - \tilde{a} \frac{L}{T_0} \cos(\frac{2\pi}{L}x) \sin(\frac{2\pi}{T_0}t) \end{pmatrix} \quad (5.3) $$

with the cross-sectional area at equilibrium and the wall stiffness given by

$$ \hat{a}_0 = \tilde{A}_0 + \tilde{a}_0 \sin \left( \frac{2\pi}{L} x \right) \quad \text{and} \quad \hat{K} = \left( 1 + 0.01 \sin \left( \frac{2\pi}{L} x \right) \right) \tilde{K}. $$

In particular we take $\tilde{A} = \tilde{A}_0 = 4.0 \times 10^{-4} \text{m}^2$, $\tilde{a} = \tilde{a}_0 = 4.0 \times 10^{-5} \text{m}^2$, $\tilde{q} = 0 \text{m}^3 \text{s}^{-1}$, $\tilde{K} = 2500 \text{kPa}$, $T_0 = 1.0 \text{s}$, the length of the vessel $L = 1.0 \text{m}$ and $\rho = 1050.0 \text{kg m}^{-3}$. We take CFL = 0.25 and exploit the MINMOD slope. Being (5.3) a periodic solution, as boundary conditions we use periodic ones. As initial condition we take $\hat{Q}$ at initial time $t = 0$.
At this stage, let us give the definition of the p-norm of the global error $E^n$,

$$ ||E(t^n, \Delta x)||_p = (\Delta x \sum_{i=-\infty}^{\infty} |v_i^n - v(x_i, t^n)|^p)^{\frac{1}{p}}, $$

where $v_i^n$ is the numerical solution and $v(x_i, t^n)$ is the exact solution computed in $(x_i, t^n)$. Note that we will use $p = 1$, $p = 2$ and $p = +\infty$. In our case we take either $v = A$ or $v = Au$. Given an increasing sequence of mesh $M_k$, with $k = 1, \dots, N$, and respective dimension $\Delta x_k$, we can now define the empirical order of accuracy $p_{k+1}$ as:

$$ p_{k+1} = \frac{\ln\left(\frac{E_{k+1}(T_{\text{out}}, \Delta x_{k+1})}{E_k(T_{\text{out}}, \Delta x_k)}\right)}{\ln\left(\frac{\Delta x_{k+1}}{\Delta x_k}\right)}, $$

with ending time $t = T_{\text{out}}$. $p_k$ should tend to the theoretical order of accuracy $p$, for sufficiently large $k$.

At last, in table 1, we show the errors and order of convergence of the variables $A$ and $q = Au$ at the ending time $T_{\text{out}} = 0.08\text{s}$. Indeed, we observe that the non-well-balanced numerical scheme described in section 4.1 reaches the second order of accuracy.

<table><caption>Table 1: Errors and empirical convergence rates for norms L<sup>1</sup>, L<sup>2</sup> and L<sup>∞</sup>. Mesh of size M = (16, 32, 64, 128, 256). Second-order scheme of section 4.1.</caption><thead><tr><th>Variable</th><th>Mesh M</th><th>err L<sup>1</sup></th><th>err L<sup>2</sup></th><th>err L<sup>∞</sup></th><th>O(L<sup>1</sup>)</th><th>O(L<sup>2</sup>)</th><th>O(L<sup>∞</sup>)</th></tr></thead><tbody><tr><td rowspan="5">Area A</td><td>16</td><td>0.1553 × 10<sup>-5</sup></td><td>0.1722 × 10<sup>-5</sup></td><td>0.2233 × 10<sup>-5</sup></td><td>—</td><td>—</td><td>—</td></tr><tr><td>32</td><td>0.0392 × 10<sup>-5</sup></td><td>0.0463 × 10<sup>-5</sup></td><td>0.0761 × 10<sup>-5</sup></td><td>1.9872</td><td>1.8948</td><td>1.5535</td></tr><tr><td>64</td><td>0.0094 × 10<sup>-5</sup></td><td>0.0122 × 10<sup>-5</sup></td><td>0.0222 × 10<sup>-5</sup></td><td>2.0663</td><td>1.9205</td><td>1.7761</td></tr><tr><td>128</td><td>0.0024 × 10<sup>-5</sup></td><td>0.0032 × 10<sup>-5</sup></td><td>0.0057 × 10<sup>-5</sup></td><td>1.9819</td><td>1.9284</td><td>1.9585</td></tr><tr><td>256</td><td>0.0006 × 10<sup>-5</sup></td><td>0.0008 × 10<sup>-5</sup></td><td>0.0014 × 10<sup>-5</sup></td><td>1.9606</td><td>1.9728</td><td>2.0130</td></tr><tr><td rowspan="5">Flow q</td><td>16</td><td>0.6980 × 10<sup>-5</sup></td><td>0.6994 × 10<sup>-5</sup></td><td>0.7699 × 10<sup>-5</sup></td><td>—</td><td>—</td><td>—</td></tr><tr><td>32</td><td>0.1581 × 10<sup>-5</sup></td><td>0.1610 × 10<sup>-5</sup></td><td>0.1958 × 10<sup>-5</sup></td><td>2.1420</td><td>2.1186</td><td>1.9755</td></tr><tr><td>64</td><td>0.0332 × 10<sup>-5</sup></td><td>0.0340 × 10<sup>-5</sup></td><td>0.0422 × 10<sup>-5</sup></td><td>2.2513</td><td>2.2431</td><td>2.2156</td></tr><tr><td>128</td><td>0.0075 × 10<sup>-5</sup></td><td>0.0079 × 10<sup>-5</sup></td><td>0.0112 × 10<sup>-5</sup></td><td>2.1383</td><td>2.1023</td><td>1.9169</td></tr><tr><td>256</td><td>0.0018 × 10<sup>-5</sup></td><td>0.0019 × 10<sup>-5</sup></td><td>0.0030 × 10<sup>-5</sup></td><td>2.1066</td><td>2.0739</td><td>1.9024</td></tr></tbody></table>

### Riemann problem: the ideal tourniquet.

Since Riemann problems are simple and idealized test cases but still useful to give a better understanding of the numerical schemes, here we present the following problem, the ideal tourniquet, for which we refer to [15, 36]. A tourniquet is placed and then immediately removed. As such, as initial data we consider

$$ Q(x, t = 0) = \begin{cases} Q_L & \text{if } x < L/2 \\ Q_R & \text{if } x > L/2 \end{cases} $$

with initial velocity $u_L = u_R = 0\text{m/s}$, initial radius $R_L = 5 \times 10^{-3}\text{m}$, $R_R = 4 \times 10^{-3}\text{m}$ and initial area computed as $A = \pi R^2$. Regarding the other parameters we take $K = \frac{1}{\sqrt{\pi}} \times 10^7\text{Pa/m}$, the length of the vessel $L = 0.08\text{m}$ and $\rho = 1060.0\text{kg m}^{-3}$. For the first and second-order schemes we respectively use CFL<sub>l</sub> = 0.45 and CFL<sub>r</sub> = 0.25. Finally the ending time is given by $T_{\text{out}} = 0.005\text{s}$ and once again we exploit the MINMODE slope. In figure 2 we compare the result for the first and second-order schemes against the exact solution. Respectively on the left and the right we used $M = 100$ and $M = 500$ cells, where $\Delta x = \frac{L}{M}$. Both the schemes approximate the exact solution well; obviously the second-order scheme results to be less diffusive than the first-order one. On the right, we can see that the numerical solution converges to the exact one.

## 5.2 System of balance laws

For the numerical simulations in this section we refer to [20, 36]. We distinguish the well-balanced schemes based on the approximate Riemann solver and the hydrostatic reconstruction respectively as WB-ARS and WB-HR. If second-order
<figure><img src="image_3.png" /><figcaption>Figure 2: Ideal tourniquet problem, radius (left) and flow (right). First-order (-o symbol), second-order (-* symbol) and exact (red line) solution.</figcaption></figure>

<table><caption>Table 2: Parameters values for test cases in section 5.2.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>L</td><td>10 [cm]</td></tr><tr><td>R<sub>in</sub></td><td>0.5 [cm]</td></tr><tr><td>K<sub>in</sub></td><td>1 &times; 10<sup>5</sup> [gcm<sup>-2</sup>s<sup>-2</sup>]</td></tr><tr><td>&rho;</td><td>1 [gcm<sup>-3</sup>]</td></tr></tbody></table>

accurate, we call them WB-ARS2 and WB-HR2. Let us give some numerical details we use for the subsequent test problems, unless it is specified otherwise. In general as initial condition we take

$$ \begin{cases} A(x, t = 0) = A_0(x) \\ q(x, t = 0) = 0. \end{cases} $$

As far as the cross-sectional radius at equilibrium $R_0$ and the wall rigidity $K$ are concerned, we use the following relations

$$ R_0(x) = \begin{cases} R_{in} & \text{if } x < x_s \text{ or } x > x_f \\ R_{in}(1 - \frac{\Delta G}{2}(1 + \cos(\pi + 2\pi \frac{x-x_s}{x_f-x_s}))) & \text{if } x_s \le x \le x_f, \end{cases} \quad (5.4) \\ K(x) = \begin{cases} K_{in} & \text{if } x < x_s \text{ or } x > x_f \\ K_{in}(1 + \frac{\Delta G}{2}(1 + \cos(\pi + 2\pi \frac{x-x_s}{x_f-x_s}))) & \text{if } x_s \le x \le x_f, \end{cases} $$

where $A_0 = \pi R_0^2$, $x_s = \frac{3L}{10}$, $x_f = \frac{7L}{10}$ and $\Delta G \in \{1\%, 10\%, 30\%, 60\%\}$. The other parameters values can be found in table 2. Regarding the boundary conditions, we impose the following flow at the inlet of the domain,

$$ q_{in} = ShapA_{in}c_{in} $$

where a value for $A_{in}$ consistent with the Shapiro number (Shap) has been estimated in [20] to be

$$ A_{in} = A_0(x = 0)(1 + Shap)^2. $$
<table><caption>Table 3: Errors and empirical convergence rates for norms $L^1$, $L^2$ and $L^\infty$ for the energy discharge $E$ and flow $q$. Mesh of size $M = (16, 32, 64, 128)$. WB-ARS2 (top) and WB-HR2 (bottom) methods.</caption><thead><tr><th>Numerical method</th><th>Variable</th><th>Mesh M</th><th>err L<sup>1</sup></th><th>err L<sup>2</sup></th><th>err L<sup>∞</sup></th><th>O(L<sup>1</sup>)</th><th>O(L<sup>2</sup>)</th><th>O(L<sup>∞</sup>)</th></tr></thead><tbody><tr><td rowspan="8">WB-ARS2</td><td rowspan="4">Energy</td><td>16</td><td>45.9648</td><td>20.2167</td><td>15.2702</td><td>—</td><td>—</td><td>—</td></tr><tr><td>32</td><td>12.1161</td><td>5.3498</td><td>4.5128</td><td>1.9236</td><td>1.9180</td><td>1.7586</td></tr><tr><td>64</td><td>2.9137</td><td>1.3085</td><td>1.1747</td><td>2.0560</td><td>2.0315</td><td>1.9418</td></tr><tr><td>128</td><td>0.7133</td><td>0.3244</td><td>0.2970</td><td>2.0302</td><td>2.0121</td><td>1.9835</td></tr><tr><td rowspan="4">Flow</td><td>16</td><td>0.0978</td><td>0.0525</td><td>0.0403</td><td>—</td><td>—</td><td>—</td></tr><tr><td>32</td><td>0.0280</td><td>0.0147</td><td>0.0112</td><td>1.8034</td><td>1.8334</td><td>1.8475</td></tr><tr><td>64</td><td>0.0070</td><td>0.0037</td><td>0.0029</td><td>2.0016</td><td>1.9787</td><td>1.9542</td></tr><tr><td>128</td><td>0.0018</td><td>0.0010</td><td>0.0007</td><td>1.9621</td><td>1.9439</td><td>1.9875</td></tr><tr><td rowspan="8">WB-HR2</td><td rowspan="4">Energy</td><td>16</td><td>44.4265</td><td>18.8607</td><td>13.4593</td><td>—</td><td>—</td><td>—</td></tr><tr><td>32</td><td>11.4229</td><td>5.0448</td><td>4.3839</td><td>1.9595</td><td>1.9025</td><td>1.6183</td></tr><tr><td>64</td><td>2.7538</td><td>1.2443</td><td>1.1764</td><td>2.05254</td><td>2.0194</td><td>1.8978</td></tr><tr><td>128</td><td>0.7072</td><td>0.3185</td><td>0.3001</td><td>1.9613</td><td>1.9660</td><td>1.9711</td></tr><tr><td rowspan="4">Flow</td><td>16</td><td>0.0973</td><td>0.0523</td><td>0.0377</td><td>—</td><td>—</td><td>—</td></tr><tr><td>32</td><td>0.0281</td><td>0.0148</td><td>0.0111</td><td>1.7910</td><td>1.8237</td><td>1.7686</td></tr><tr><td>64</td><td>0.0070</td><td>0.0037</td><td>0.0029</td><td>2.0057</td><td>1.9883</td><td>1.9403</td></tr><tr><td>128</td><td>0.0018</td><td>0.0010</td><td>0.0007</td><td>1.9634</td><td>1.9448</td><td>1.9854</td></tr></tbody></table>

The Shapiro number is the analogous to the Mach number for the compressible Euler equations and it is defined as $Shap = \frac{u}{c}$.

In particular we take $Shap = 0, 10^{-2}, 10^{-3}$. Let us note that with $Shap = 10^{-2}$ we are already in the subsonic regime, thus it could be interesting to use an implicit scheme in this case. Moreover, usually in the arteries the average value for the Shapiro number is indeed of order $Shap = 10^{-2}$. We refer again to [20] for more details. Then, we can find the boundary value's for the cross-sectional area $A_{in}$ exploiting the right Riemann invariant, i.e. $I^+ = u - 4c$, and imposing $I^+(Q_1) = I^+(Q_{in})$, where the subscript 1 indicates the values in the first cell of the computational domain. Whereas, regarding the right boundary condition, we impose

$$A_{out} = A_{0}(x = L)(1 + Shap)^{2},$$

and then the flow value $q_{out}$ is found exploiting the left Riemann invariant, namely $I^-(Q_{end}) = I^-(Q_{out})$, where with $Q_{end}$ we mean the value in the last cell. For more details about this test cases we refer again to [20].

### Stationary solution.

First of all, to assess the well-balanced property, we take $Shap = 0$, and check that the numerical schemes preserve the stationary solution $A = A_0$, $q = 0$. Indeed, we observe that they maintain it up to an error of order $10^{-12}$.

### Convergence study.

Then, in order to check that also the well-balanced schemes reach the right order of accuracy, we take $Shap = 10^{-2}$, $\Delta G = 10\%$ and final time $tEnd = 100.0s$. Note that here the exact solution is a steady state with non-zero velocity, namely it is given by

$$ \begin{cases} q_{ex} = q_{in} \\ E_{ex} = \frac{q_{in}^2}{2A_{out}^2} + \frac{K(end)}{\rho}(\sqrt{A_{out}} - \sqrt{A_{0(end)}}) \end{cases} \qquad (5.5) $$

with $E$ the energy and $q_{ex} \neq 0$. Indeed our schemes are able to preserve stationary solutions with only zero-velocity, thus the numerical solutions should converge to (5.5) when refining the mesh.

In table 3 the numerical errors and orders of convergence are exhibited in norms $L^1$, $L^2$ and $L^\infty$ for both WB-ARS2 and WB-HR2; the results seem to be satisfying.

### Wave propagation test case.
<table><caption>Table 4: Parameters values for test cases in section 5.2.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>L</td><td>0.16 [m]</td></tr><tr><td>R<sub>in</sub></td><td>4.0 × 10<sup>-3</sup> [m]</td></tr><tr><td>ΔR</td><td>1.0 × 10<sup>-3</sup> [m]</td></tr><tr><td>K</td><td><sup>1</sup>⁄<sub>√π</sub> × 10<sup>8</sup> [Pam]</td></tr><tr><td>ρ</td><td>1060 [Kgm<sup>-3</sup>]</td></tr></tbody></table>

In this test we assess the wave-capturing properties of the well-balanced schemes. We suppose that a single wave propagates in the vessel, with parameters defined by (5.4) and table 2. Namely, we impose the following unsteady inlet flow

$$q_{\text{in}}(t) = \begin{cases} q_{\text{pulse}}\sin(2\pi\frac{t}{t_{\text{pulse}}}) & \text{if } t \le \frac{t_{\text{pulse}}}{2} \\ 0 & \text{otherwise} \end{cases}$$

where once again we define $q_{\text{pulse}}$ as $q_{\text{pulse}} = \text{ShapA}_{\text{in}}c_{\text{in}}$. For the right boundary condition we simply use the transmissive one. Finally $t_{\text{pulse}} = 0.04\text{s}$ and as ending time we take $t_{\text{Out}} = 0.045\text{s}$. We compute a reference solution with the WB-HR2 and $M = 2048$ cells. In figure 3 we insert the results only for WB-ARS and WB-ARS2 as the ones obtained with WB-HR/2 are similar. Of course, solutions obtained with first-order schemes are more diffusive than the ones found exploiting higher order methods, but both outputs tend to the reference one. We also observe that there are no spurious oscillations in the results.

*Propagation of a pulse to/from an expansion.*

Here we want to consider two different cases, a pulse which propagates to and from an expansion. In the former case as initial radius we take

$$R(x, t = 0) = \begin{cases} R_0(x)(1 + \varepsilon \sin(\frac{100}{20L}\pi(x - \frac{65L}{100}))) & \text{if } \frac{65L}{100} \le x \le \frac{85L}{100} \\ R_0 & \text{if } \text{else,} \end{cases} \quad (5.6)$$

while in the second one

$$R(x, t = 0) = \begin{cases} R_0(x)(1 + \varepsilon \sin(\frac{100}{20L}\pi(x - \frac{15L}{100}))) & \text{if } \frac{15L}{100} \le x \le \frac{35L}{100} \\ R_0 & \text{if } \text{else} \end{cases} \quad (5.7)$$

where $\varepsilon = 5.0 \times 10^{-3}$. In this last numerical problem we assume the wall rigidity $K$ to be constant, while the radius at equilibrium is given by

$$R_0(x) = \begin{cases} R_{in} + \Delta R & \text{if } x < x_s \\ R_{in} + \frac{\Delta R}{2}(1 + \cos(\pi \frac{x-x_s}{x_f-x_s})) & \text{if } x_s \le x \le x_f \\ R_{in} & \text{if } \text{else} \end{cases}$$

where $A_0 = \pi R_0^2$, $x_s = \frac{19L}{40}$ and $x_f = \frac{L}{2}$. Note that $A_1 > A_{\text{end}}$. The other parameters values can be found in table 4. We use transmissive boundary conditions.

In figure 4 we present the outputs for this two problems; for the results on the left and right we respectively use initial conditions (5.6) and (5.7). We compare the outputs of WB-HR2 for $M = 200$ cells against a reference solution attained with WB-HR2 and $M = 2048$ cells. Once again we show the results of only one of the numerical schemes as their outputs are very similar. Indeed, the numerical solutions are satisfying and comparable with the ones of [15, 36].

# 6 Concluding remarks

In this paper we presented two different second-order well-balanced numerical schemes for the 1D blood flow equations, where the source term is due to varying mechanical and geometrical properties. By well-balanced we mean that the numerical method is able to preserve the zero-velocity "man at eternal rest" stationary solution. The first scheme is based on an
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 3: Wave propagation problem. We used $\Delta G = 10\%$ (top), $\Delta G = 30\%$ (middle), $\Delta G = 60\%$ (bottom) and $Shap = 10^{-2}$ (left), $Shap = 10^{-3}$ (right). Comparison among WB-ARS (-o blue symbol), WB-ARS2 (-* green symbol) and reference (red line) solution. $M = 100$ cells.</figcaption></figure>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 4: Propagation of a pulse to (left) and from (right) an expansion. Radius minus radius at equilibrium at time $t = 0.0s$ (green), $t = 0.002s$ (blue), $t = 0.006s$(red). Comparison between WB-HR2 ($^*$ symbol) and reference (line) solution. $M = 200$ cells.</figcaption></figure>

approximate Riemann solver and the numerical source is defined in such a way that it is consistent in the integral sense with the source term, while the second method exploits the hydrostatic reconstruction approach. Only the Lagrangian step has to be modified in order to satisfy the well-balanced property. Both the numerical schemes proved to be satisfying and their results are almost identical. On one hand, future works are expected to deal with an implicit formulation of the Lagrangian step as it would speed up the numerical methods. Indeed, the arteries are known to have an average Shapiro number of order $10^{-2}$ and this could imply a restriction on the time step value. On the other hand, efforts could be done to obtain fully well-balanced schemes able to preserve also stationary solutions for non-zero velocity. In this regard, we refer to the works [5]-[25] for fully well-balanced Lagrange-projection schemes applied to the shallow water system at first and high order of accuracy respectively. Another interesting reference in this framework is [20], in which the authors considered the low-Shapiro number steady states, which could be more easily preserved than the classical steady states with non-zero velocity.

**Acknowledgments.** This work was partially supported by a grant from Région Île-de-France. We also would like to thank Professor E. F. Toro for his suggestions about this manuscript.

## References

[1] Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein and Benoit Perthame. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. Siam Journal on Scientific Computing, 25, 2050-2065, 2004. https://doi.org/10.1137/S1064827503431090

[2] Michaël Baudin, Christophe Berthon, Frédéric Coquel, Roland Masson, Quang HuyTran. A relaxation method for two-phase flow models with hydrodynamic closure law. Numerische Mathematik. 99, 411-440, 2005. 10.1007/s00211-004-0558-1.

[3] Christophe Berthon, Françoise Foucher. Efficient well-balanced hydrostatic upwind schemes for shallow-water equations. Journal of Computational Physics. 231, 4993-5015, 2012. 10.1016/j.jcp.2012.02.031.
[4] François Bouchut. Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources. Frontiers in mathematics, 2004. 10.1007/b93802.

[5] Manuel J. Castro Díaz, Christophe Chalons and Tomás Morales De Luna. A fully well-balanced Lagrange-Projection type scheme for the Shallow-water equations. SIAM J. Numer. Anal., 56(5), 3071–3098, 2018. 10.1137/17M1156101.

[6] Manuel J. Castro, Tomás Morales de Luna and Carlos Parés. Well-Balanced Schemes and Path-Conservative Numerical Methods. Handbook of Numerical Analysis, Vol. 18, Pages 131-175, 2017. https://doi.org/10.1016/bs.hna.2016.10.002

[7] Manuel J. Castro Díaz, Alberto Pardo Milanés and Carlos Parés. Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique. Mathematical Models and Methods in Applied Sciences. 5, 2055-2113, 2007. 10.1142/S021820250700256X.

[8] Manuel J. Castro and Carlos Parés. Well-Balanced High-Order Finite Volume Methods for Systems of Balance Laws. Journal of Scientific Computing. 82, 48, 2020. 10.1007/s10915-020-01149-5.

[9] Christophe Chalons, Mathieu Girardin, Samuel Kokh. An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes. Journal of Computational Physics, Elsevier, 335, pp.885-904. 10.1016/j.jcp.2017.01.017.

[10] Christophe Chalons, Mathieu Girardin, Samuel Kokh. An All-Regime Lagrange-Projection Like Scheme for the Gas Dynamics Equations on Unstructured Meshes. Communications in Computational Physics. 20, 2014. 10.4208/cicp.260614.061115a.

[11] Christophe Chalons, Pierre Kestener, Samuel Kokh, and Maxime Stauffert. A large time-step and well-balanced Lagrange-Projection type scheme for the Shallow-water equations. Communications in Mathematical Sciences. 15, 2016. 10.4310/CMS.2017.v15.n3.a9.

[12] Christophe Chalons, Samuel Kokh, Mathieu Girardin. Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms. SIAM Journal on Scientific Computing. 35, 2013. 10.1137/130908671.

[13] Frederic Coquel, Edwige Godlewski, Benoit Perthame, Arun In, and Paul Rascle. Some new Godunov and relaxation methods for two-phase flow problems. Godunov methods (Oxford, 1999), pages 179–188, 2001. https://doi.org/10.1007/978-1-4615-0663-8_18

[14] Frederic Coquel and Benoit Perthame. Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics. SIAM Journal on Numerical Analysis, 35(6):2223–2249, 1998. https://doi.org/10.1137/S0036142997318528.

[15] Olivier Delestre, Pierre-Yves Lagrée. A "well-balanced" finite volume scheme for blood flow simulation. International Journal for Numerical Methods in Fluids, Wiley, 72 (2), pp.177-205, 2013. 10.1002/fld.3736.

[16] Frédéric Duboc, Cédric Enaux, Stéphane Jaouen, Hervé Jourdren, Marc Wolff. High-order dimensionally split Lagrange-remap schemes for compressible hydrodynamics. Comptes Rendus Mathematique Volume 348, Issues 1-2, Pages 105-110, January 2010. 10.1016/j.crma.2009.12.008.

[17] Luca Formaggia, Alfio Quarteroni, Alessandro Veneziani. Cardiovascular Mathematics: Modeling and Simulation of the Circulatory System. Springer-Verlag: Italia, Milano, 2009. 10.1007/978-88-470-1152-6.

[18] Gérard Gallice. Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source. C. R. Math. Acad. Sci. Paris 334, no. 8, pp. 713-716, 2002. 10.1016/S1631-073X(02)02307-5.

[19] Gérard Gallice. Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates. Numer. Math. 94 no. 4, pp. 673-713, 2003. 10.1007/s00211-002-0430-0
[20] Arthur R. Ghigo, Olivier Delestre, Jose-Maria Fullana, Pierre-Yves Lagrée. Low-Shapiro hydrostatic reconstruction technique for blood flow simulation in large arteries with varying geometrical and mechanical properties. Journal of Computational Physics, Elsevier, 331, pp.108 - 136, 2017. 10.1016/j.jcp.2016.11.032.

[21] Victor Michel-Dansac, Christophe Berthon, Stéphane Clain, Françoise Foucher. A well-balanced scheme for the shallow-water equations with topography. Computers and Mathematics with Applications, Elsevier, 72, pp.568 - 593, 2016. 10.1016/j.camwa.2016.05.015.

[22] Victor Michel-Dansac, Christophe Berthon, Stéphane Clain, Françoise Foucher. A well-balanced scheme for the shallow-water equations with topography or Manning friction. Journal of Computational Physics, Elsevier, 335, pp.115-154, 2017. 10.1016/j.jcp.2017.01.009.

[23] Sigal Gottlieb and Chi-Wang Shu. Total variation diminishing RUNGE-KUTTA schemes. Mathematics of Computation. 67, 1996. 10.1090/S0025-5718-98-00913-2.

[24] Amiram Harten, Peter D. Lax, and Bram Van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Review, 25:pp. 35-61, 1983.

[25] Tomás Morales De Luna, Manuel J. Castro Díaz and Christophe Chalons. High order fully well-balanced Lagrange-Projection scheme for Shallow-water. Commun. Math. Sci., Vol. 18, No. 3, pp. 781-807, 2020. https://dx.doi.org/10.4310/CMS.2020.v18.n3.a9

[26] Lucas O. Müller, Carlos Parés and Eleuterio F. Toro. Well-balanced high-order numerical schemes for one-dimensional blood flow in vessels with varying mechanical properties. Journal of Computational Physics 242, 53-85, 2013. 10.1016/j.jcp.2013.01.050.

[27] Lucas O. Müller, and Eleuterio F. Toro. Well-balanced high-order solver for blood flow in networks of vessels with variable properties. Int. J. Numer. Meth. Biomed. Engng. 29:1388-1411, 2013. 10.1002/cnm.2580.

[28] Kambiz Salari, Patrick Knupp. Code Verification by the Method of Manufactured Solutions. Report, 2000. https://digital.library.unt.edu/ark:/67531/metadc702130/

[29] Ion Suliciu. On the thermodynamics of fluids with relaxation and phase transitions. Fluids with relaxation. Int. J. Engag. Sci. 36, pp. 921-947, 1998.

[30] Eleuterio F. Toro and Annunziato Siviglia. PRICE: Primitive Centred Schemes for Hyperbolic Systems. Int. J. Numer. Meth. in Fluids, 42:1263-1291, 2003. 10.1002/fld.491.

[31] Eleuterio F. Toro and Annunziato Siviglia. Simplified blood flow model with discontinuous vessel properties: analysis and exact solutions. Modeling, Simulation and Applications, 5, 2011. 10.1007/978-88-470-1935-5_2.

[32] Eleuterio F. Toro and Annunziato Siviglia. Flow in collapsible tubes with discontinuous mechanical properties: mathematical model and exact solutions. Communications in Computational Physics; 13(2):361-385, 2013. 10.4208/cicp.210611.240212a.

[33] Eleuterio F. Toro. Lecture notes on computational haemodynamics. Mathematics Department, University of Trento, Italy, 2017.

[34] Eleuterio F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics, Third Edition. Springer-Verlag, 2009. 10.1007/b79761_5.

[35] Eleuterio F. Toro. Brain venous haemodynamics, neurological diseases and mathematical modelling. A review. Applied Mathematics and Computation, 272, 542-579, 2016. 10.1016/j.amc.2015.06.066.

[36] Zhenzhen Wang, Gang Li and Olivier Delestre. Well-balanced finite difference weighted essentially non-oscillatory schemes for the blood flow model. Int. J. Numer. Meth. Fluids, 82, 2016. 10.1002/fld.4232.