# Numerical approximation and uncertainty quantification for arterial blood flow models with viscoelasticity

Christophe Chalons, Alessia Del Grosso, Eleuterio Francisco Toro

► To cite this version:

Christophe Chalons, Alessia Del Grosso, Eleuterio Francisco Toro. Numerical approximation and uncertainty quantification for arterial blood flow models with viscoelasticity. Journal of Computational Physics, 2022, 457, 10.1016/j.jcp.2022.111071. hal-03993016

HAL Id: hal-03993016

https://hal.science/hal-03993016

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical approximation and uncertainty quantification for arterial blood flow models with viscoelasticity

Christophe Chalons, Alessia Del Grosso†and Eleuterio F. Toro ‡

December 18, 2019

**Abstract.** The importance of the study of the blood flow equations is widely recognized as it is a tool to understand the circulatory system. Arteries and veins result to have both elastic and viscous behaviour. Models for the first case are much more studied as they result to be simpler and still satisfying if compared to experimental data. In this paper, we consider a model which encompasses both the elastic and viscoelastic response in arterial walls, respectively leading to a conservative and a non-conservative system. We present a second-order scheme based on the first-order Price-T scheme and the MUSCL-Hancock strategy. This new approach automatically adapts to the above conservative and non-conservative cases.

Then, we perform a Sensitivity Analysis (SA) based on the Continuous Sensitivity Equation Method (CSEM), whose aim is the study of how changes in the inputs of a model can affect its outputs. In particular, the sensitivity is defined as the derivative (with respect to an uncertain parameter a) of the solution of the system taken into consideration. Since the CSEM cannot be directly applied to discontinuous solutions, we add a source term to compensate the spikes associated to the Dirac delta functions that can arise in the sensitivity variables.

One of the main applications of SA is uncertainty quantification, which is investigated for Riemann problems as well as for a network of 37 arteries. Details on junctions for coupling two or more vessels are also given.

## 1 Introduction and governing equations

The simulation of blood flows in human arteries has been a very active topic of research over the last decades. Indeed, recovering experimental data about the cardiovascular system proves to be extremely complicated and this motivates the numerical studies that have been done so far. The numerical outputs can help to expand the knowledge in this field, for instance on many pathologies such as hypertension or presence of occlusions.

Models in one-dimensional (1D) space dimension are often used as an alternative to the more complex three-dimensional (3D) fluid-structure interaction models. Indeed, they are able to provide a good description of the propagation of pressure waves in arteries with a low computational cost compared to 3D models, which is especially interesting when considering global human circulatory systems represented by a network of vessels. We refer for instance the reader to [3], [12], [19], [27] and [39] for more details.

In this paper, we are interested in the numerical approximation of the solutions of the following 1D system composed of two partial differential equations,

$$ \begin{cases} \partial_t A + \partial_x (Au) = 0 \\ \partial_t (Au) + \partial_x (Au^2) + \frac{A}{\rho} \partial_x p = f, \end{cases} \qquad (1) $$

consisting of the mass conservation and momentum balance equations. Here $x$ represents the axial coordinate along the longitudinal axis of the vessel and $t > 0$ is the time. $A(x, t) > 0$ is the cross-sectional area of the vessel, $u(x, t)$ the averaged velocity of blood at cross section and consequently $q = Au$ is the blood flow, $p$ is the pressure and $\rho$ represents the blood density. At last, $f = -Ru$ accounts for the friction forces where $R > 0$ is the viscous resistance given by

$$ R = 22\pi \frac{\mu}{\rho}, $$

<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>

<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>‡Laboratory of Applied Mathematics, Department of Civil and Environmental Engineering, University of Trento, Trento, Italy. E-mail: eleuterio.toro@unitn.it</footnote>
where $\mu$ is the dynamic viscosity. We assume that the initial area $A(x, t = 0)$ and velocity $u(x, t = 0)$ are given at time $t = 0$. For more details about the derivation of system (1), we refer to [18] and [37].

In order to close this system, a relationship between the pressure and the cross-sectional area is needed. This is the so-called tube law. Depending on the underlying mechanical considerations and elastic or viscoelastic assumptions on the response of the artery, the form of the closure relation and consequently the nature of system (1) can significantly differ. In this paper, we will more precisely consider the following two relations, for which we refer to [37] and [25].

*Elastic tube law.* The first closure relation has been proposed for example in [29] and assumes that the arterial walls are purely elastic, leading to the following relation

$$p = p_{ext} + \psi \qquad (2)$$

where $p_{ext}$ is the external pressure acting on the vessel walls and $\psi$ is the transmural pressure, namely the difference between the internal and external pressures. Since we are interested in blood flow in arteries, the transmural pressure reads

$$\psi(A, A_0, K) = K\phi(A, A_0)$$

with

$$\phi(A, A_0) = \sqrt{\frac{A}{A_0}} - 1 \quad \text{and} \quad K(x) = \frac{4}{3}\sqrt{\frac{\pi}{A_0}}Eh_0,$$

where $h_0$ is the wall thickness, $E$ the Young modulus and $A_0$ is the cross-sectional area of the vessel at equilibrium. Note that for the sake of simplicity, we will assume that $h_0, A_0$ and $E$ are constant along a single vessel. Hence, substituting (2) into (1), we obtain the ultimate form of the model,

$$\left\{ \begin{array}{l} \partial_t A + \partial_x(Au) = 0 \\ \partial_t(Au) + \partial_x(Au^2 + \gamma A^{\frac{3}{2}}) = f, \end{array} \right. \qquad (3)$$

with

$$\gamma = \frac{K}{3\rho\sqrt{A_0}},$$

and where the left-hand side of the momentum balance equation is now in conservation form. Easy calculations show that the two eigenvalues of the system are $\lambda_e^{\pm} = u \pm c$, where the wave speed is defined by

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

Consequently the convective part of (3) is strictly hyperbolic as $\lambda_e^{\pm}$ are real and distinct as long as $c \neq 0$. Finally, both the two characteristic fields are genuinely non-linear and the Riemann invariants associated with $\lambda_e^{\pm}$ are respectively given by $I_e^{-} = u + 4c$ and $I_e^{+} = u - 4c$. For more details refer to [36].

*Viscoelastic tube law.* The second closure relation has been proposed for instance in [25] and considers that the viscoelastic response in arterial walls cannot be neglected. In this case, a viscoelastic term is added in the tube law, leading to the more general relation

$$p(x, t) = p_{ext} + \psi(A, A_0, K) + \varphi(A, A_0)\partial_t A, \qquad (4)$$

where we take

$$\varphi(A, A_0) = \frac{\Gamma}{A_0\sqrt{A}}.$$

The coefficient $\Gamma$ depends on the viscoelastic properties of the vessel's wall and we note in particular that taking $\Gamma = 0$ boils down to considering an elastic closure. In practice, we will refer to [1] and take

$$\Gamma = \frac{2}{3}sh_0\sqrt{\pi},$$

where $s$ is the viscosity of the (silicone) vessel.
At this stage, it is important to note that using the mass conservation equation of (1), (4) is equivalent to

$$p = p_{ext} + \psi(A, A_0, K) - \varphi(A, A_0)\partial_x q,$$

so that

$$\partial_x p = K \partial_A \phi \partial_x A - \partial_A \varphi \partial_x A \partial_x q - \varphi \partial_x^2 q.$$

Therefore, it is observed that a second-order spatial derivative now appears in the momentum balance equation of (1), making the viscoelastic system of advection-diffusion-reaction type unlike the elastic model (3) which is an advection-reaction system. However, one of the two goals of the present contribution is to propose an appropriate second-order accurate numerical method to approximate the elastic and viscoelastic blood flow models in the same formalism. Therefore, we suggest to proceed as in the paper of Montecinos et al. [25], and to approximate the original parabolic (advection-diffusion-reaction) viscoelastic problem by a hyperbolic (advection-reaction) system with stiff source term using a relaxation technique. More precisely, we will relax system (1) by following the constitutive Cattaneo's law [7]

$$\partial_t \Psi - \frac{1}{\varepsilon} \partial_x q = - \frac{1}{\varepsilon} \Psi,$$

where $\varepsilon > 0$ is a parameter such that in the asymptotic limit $\varepsilon \to 0$, the new variable $\Psi$ converges towards the original gradient $\partial_x q$. Therefore, with this new variable we can rewrite the pressure gradient as

$$\partial_x p = K \partial_A \phi \partial_x A - \partial_A \varphi \partial_x A \Psi - \varphi \partial_x \Psi$$

and reformulate system (1) as

$$\left\{
\begin{array}{l}
\partial_t A + \partial_x q = 0 \\
\partial_t q + \partial_x \left(\frac{q^2}{A}\right) + \frac{A}{\rho}(K \partial_A \phi - \Psi \partial_A \varphi) \partial_x A - \frac{A}{\rho} \varphi \partial_x \Psi = f \\
\partial_t \Psi - \frac{1}{\varepsilon} \partial_x q = - \frac{1}{\varepsilon} \Psi.
\end{array}
\right.
\qquad (5)$$

It is thus clear that this relaxation approach allows us to reformulate the viscoelastic system while maintaining the physical meaning of the system. However, the striking difference with the elastic model (3) is that the advection part of (5) is not in conservative form. This property will be a key feature of the proposed numerical method.

As for the elastic system, it is not difficult to check that the convective part of (5) admits three eigenvalues given by $\lambda_v^{\pm} = u \pm \tilde{c}$ and $\lambda_v^0 = 0$ where the wave speed is now defined by

$$\tilde{c} = \sqrt{c^2 + \omega} \quad \text{with} \quad \omega = \frac{\varphi A}{\rho \varepsilon} + \frac{a_{\Gamma}}{2},$$

with $a_{\Gamma} = \frac{\Psi \varphi}{\rho}$. In particular, the system is strictly hyperbolic if $c^2 + \omega \ge 0$, so if

$$\frac{1}{\varepsilon} \geq -\frac{\Psi}{2A} - \frac{c^2 \rho}{\varphi A}.$$

For more details about the analysis of the system, we refer to [25]. In addition, we observe that the characteristic fields associated to $\lambda_v^{\pm}$ are genuinely non linear, while the field associated to $\lambda_v^0$ is linearly degenerate. Finally, the Riemann invariants associated with $\lambda_v^-$ and $\lambda_v^+$ are respectively

$$I_{v}^{-,1} = \frac{A}{\varepsilon} + \Psi \quad \text{and} \quad I_{v}^{-,2} = u + \int \frac{\tilde{c}}{A} dA$$

and

$$I_{v}^{+,1} = \frac{A}{\varepsilon} + \Psi \quad \text{and} \quad I_{v}^{+,2} = u - \int \frac{\tilde{c}}{A} dA$$

while for $\lambda_v^0$ we obtain the Riemann Invariants

$$I_{v}^{0,1} = p_{ext} + \psi - \varphi \Psi + \frac{1}{2} \rho u^2 \quad \text{and} \quad I_{v}^{0,1} = q,$$
see for instance [25] for more details.

For the sake of conciseness, elastic and viscoelastic systems (3) and (5) will be given the following compact form

$$ \partial_t \mathbf{Q} + \partial_x \mathbf{F}(\mathbf{Q}) + \mathbf{A}(\mathbf{Q}) \partial_x \mathbf{Q} = \mathbf{S}(\mathbf{Q}), \qquad (6) $$

where $\mathbf{Q}$ denotes the vector of unknowns, $\mathbf{F}(\mathbf{Q})$ the flux of the conservative part, $\mathbf{A}(\mathbf{Q})$ the coefficient matrix of the non conservative part and $\mathbf{S}(\mathbf{Q})$ the source term. In particular, we have

$$ \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} = \begin{pmatrix} 0 \\ f \end{pmatrix}, \qquad (7) $$

for the elastic model, and

$$ \mathbf{Q} = \begin{pmatrix} A \\ q \\ \Psi \end{pmatrix}, \quad \mathbf{F}(\mathbf{Q}) = \begin{pmatrix} q \\ \frac{q^2}{A} + \gamma A^{\frac{3}{2}} \\ -\frac{1}{\epsilon}q \end{pmatrix}, \quad \mathbf{A} = \begin{pmatrix} 0 & 0 & 0 \\ \frac{a\Gamma}{2} & 0 & -\frac{\varphi A}{\rho} \\ 0 & 0 & 0 \end{pmatrix} \quad \text{and} \quad \mathbf{S} = \begin{pmatrix} 0 \\ f \\ -\frac{1}{\epsilon}\Psi \end{pmatrix} \qquad (8) $$

for the viscoelastic one. Note that the non conservative term $\mathbf{A}(\mathbf{Q})$ and the third row of $\mathbf{F}(\mathbf{Q})$ and $\mathbf{S}(\mathbf{Q})$ do not appear in the elastic system (3). Our first objective will be to propose a second-order accurate numerical scheme to approximate the solution of (6) and with the consistency property that the scheme applied to the viscoelastic system with $\Gamma = 0$ should be equivalent to the scheme applied directly to the elastic system. As it will be seen, this is not trivial considering the finite volume framework of first-order Price-T schemes for non conservative systems coupled with a second-order MUSCL-Hancock methodology taken into account in the present paper. We also refer the reader to the following recent contributions [17, 24, 25, 26] and the references therein for related finite volume methods applied to the elastic and viscoelastic blood flow models separately.

As usual when dealing with mathematical models, the 1D systems used in the present paper to describe the human circulatory system are based on modelling assumptions and parameters to describe the flow and geometrical or physical properties of the vessels, for instance the vessel thickness $h_0$, the cross-sectional area of the vessel at equilibrium $A_0$, the Young modulus $E$ or the initial values $A(x, t = 0)$ of the area at time $t = 0$. In practice, these parameters are estimated or calibrated and it is natural to wonder how the numerical simulations are sensitive to these approximations and uncertainties. Sensitivity Analysis (SA), that is the study of how changes in the inputs of a model affect the outputs, is a powerful tool to perform Uncertainty Quantification (UQ). It consists in computing numerically the derivatives of the solutions of the system of equations with respect to one or several chosen input parameters. There exist essentially two main approaches to SA with own advantages and disadvantages depending on the target applications. A detailed comparison between the two for optimization problems is done in [20]. The first one, known as discrete sensitivity equation method, starts with discretizing the system of equations and then differentiating it to obtain an approximation of the sensitivities. Whereas, the second one consists in differentiating the equations first, obtaining in this way the system of sensitivity equations which will be then discretised and solved, see for instance [4, 10, 15, 16]. The latter is known as the differentiate-then-discretise approach or Continuous Sensitivity Equation Method (CSEM) and it is the one that will be considered here. In this context, our second objective is to introduce the use of sensitivity analysis to study the uncertainties in the solutions of blood flow equations. This will be done both for the elementary solutions of Riemann problems (which, although being the simplest non trivial initial value problem, are at the basis of many numerical methods), and a complex network of 37 arteries already simulated in [1, 24].

*Outline of the paper.* The paper is organized as follows. In the next section we introduce the Price-T scheme and our version of MUSCL-Hancock method, respectively of first and second-order of accuracy. Outputs for Riemann problems are also shown. In sections 3 and 4 we present respectively the sensitivity equations and the modified numerical methods in order to take into account the discontinuities in the state variables. Results for Riemann problems and uncertainty quantification are carried out at the end of section 4. Subsequently, in section 5 we explain how to take care of the presence of junctions among two or three vessels, a necessary step to obtain the results for the network of 37 vessels, which is presented in section 5 as well. Finally, in section 6 conclusions are drawn.
<figure><img src="image_2.png" /><figcaption>Figure 1: Staggered grid for Price-T and MUSCL-Hancock schemes, refer to [34].</figcaption></figure>

## 2 Centred numerical scheme

The aim of this section is to propose a finite volume method which can be used to resolve non homogeneous hyperbolic system (6) both in conservative ($A(Q) = 0$) and non conservative ($A(Q) \neq 0$) case. We are especially interested in second-order centred strategies based on MUSCL reconstructions. We focus on the homogeneous system

$$ \partial_t Q + \partial_x F(Q) + A(Q) \partial_x Q = 0 \qquad (9) $$

as the source term $S(Q)$ will be taken into account thanks to a classical fractional step method at first-order accuracy, and an usual Strang splitting at second-order accuracy. Note that the source term can be solved exactly or using an ODE solver. As far as the homogeneous system is concerned, observe that the proposed method should automatically reduce to a conservative one if applied to a conservative system. This is a very important request as a non conservative method applied to a conservative system could converge to a wrong solution [21]. For this purpose, we will present the first-order Price-T scheme, and a second-order extension based on the MUSCL-Hancock strategy, see [25] for more details about these schemes.

Let us now introduce some notations. When necessary, we will first use $Q = (q_1, \dots, q_m)^t$ to denote the component of the vector of unknowns. Then, we introduce a constant space step $\Delta x$ and constant time step $\Delta t$. The mesh interfaces are defined by $x_{i+1/2} = i\Delta x$ for $i \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 $Q_i^n$ of the solution in the interval $[x_{i-1/2}, x_{i+1/2})$, $i \in \mathbb{Z}$. Therefore, a piecewise constant approximate solution $x \to Q_{\Delta t, \Delta x}(x, t^n)$ of the solution $Q$ is given by

$$ Q_{\Delta t, \Delta x}(x, t^n) = Q_i^n \quad \text{for all } x \in C_i = [x_{i-1/2}; x_{i+1/2}), \quad j \in \mathbb{Z}, \quad n \in \mathbb{N}. $$

When $n = 0$, we set

$$ Q_i^0 = \frac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} Q_0(x) dx, \quad \text{for all } i \in \mathbb{Z}. $$

### 2.1 First-order Price-T scheme

In this subsection, we briefly recall the first-order Price-T scheme, which is a variant of the Force scheme for non-conservative systems. We refer the reader to [34] and [38] for more details. Force and Price-T methods are centred and therefore do not utilize wave propagation informations coming from exact or approximate Riemann solvers, which in turn can be difficult or even impossible to define. Instead, the idea at the base of these methods is to use integral relations that can be obtained from (9). Assuming that the sequence $(Q_i^n)_{i \in \mathbb{Z}}$ is known, the main steps to define the sequence $(Q_i^{n+1})_{i \in \mathbb{Z}}$ are as follows. Note that we will not show all the computations but again refer to [34] for more details. Figure 1 is also proposed to help visualizing the reasoning and the different steps.

**Step 1.** First of all and for a given $i$, the idea is to virtually solve the Riemann problems set at interfaces $x_{i-1/2}$ and $x_{i+1/2}$ and respectively associated with the left and right states $(Q_{i-1}^n, Q_i^n)$ and $(Q_i^n, Q_{i+1}^n)$, see figure 1. Up to a space translation
of $x_{i-1/2}$ and $x_{i+1/2}$ and a time translation of $t^n$, the respective solutions are denoted $(x, t) \to \hat{\mathbf{Q}}_{i-\frac{1}{2}}^{n+\frac{1}{2}}(x, t)$ and $(x, t) \to \hat{\mathbf{Q}}_{i+\frac{1}{2}}^{n+\frac{1}{2}}(x, t)$. Note that these solutions are self-similar and composed of simple waves (shocks, contact discontinuities and rarefaction waves) emanating from the interfaces $x_{i-1/2}$ and $x_{i+1/2}$, see [34] for more details.

**Step 2.** The aim of this step is to evaluate the following integral averages of the two Riemann problem solutions at time $\Delta t/2$, namely

$$ \mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}} = \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \hat{\mathbf{Q}}_{i-\frac{1}{2}}^{n+\frac{1}{2}}(x, \frac{\Delta t}{2})dx \quad \text{and} \quad \mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}} = \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \hat{\mathbf{Q}}_{i+\frac{1}{2}}^{n+\frac{1}{2}}(x, \frac{\Delta t}{2})dx. $$

These intermediate states are simply evaluated by integrating (9) on the control volume $(-\frac{\Delta x}{2}, \frac{\Delta x}{2}) \times (0, \frac{\Delta t}{2})$, and doing a linearization in which we substitute the matrix $\mathbf{A}$ with a constant matrix $\hat{\mathbf{A}}_{i\pm 1/2}$. More precisely, we obtain

$$ \mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}} = \frac{1}{2}(\mathbf{Q}_{i-1}^n + \mathbf{Q}_i^n) - \frac{1}{2}\frac{\Delta t}{\Delta x}(\mathbf{F}(\mathbf{Q}_i^n) - \mathbf{F}(\mathbf{Q}_{i-1}^n)) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{\mathbf{A}}_{i-\frac{1}{2}}(\mathbf{Q}_i^n - \mathbf{Q}_{i-1}^n), \quad (10) $$

$$ \mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}} = \frac{1}{2}(\mathbf{Q}_i^n + \mathbf{Q}_{i+1}^n) - \frac{1}{2}\frac{\Delta t}{\Delta x}(\mathbf{F}(\mathbf{Q}_{i+1}^n) - \mathbf{F}(\mathbf{Q}_i^n)) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{\mathbf{A}}_{i+\frac{1}{2}}(\mathbf{Q}_{i+1}^n - \mathbf{Q}_i^n), \quad (11) $$

where we set

$$ \hat{\mathbf{A}}_{i-\frac{1}{2}} = \mathbf{A}(\frac{1}{2}(\mathbf{Q}_{i-1}^n + \mathbf{Q}_i^n)) \quad \text{and} \quad \hat{\mathbf{A}}_{i+\frac{1}{2}} = \mathbf{A}(\frac{1}{2}(\mathbf{Q}_i^n + \mathbf{Q}_{i+1}^n)). $$

**Steps 3 and 4.** Equipped with these updated values at intermediate time and on a staggered grid, we now follow the same steps as before to evaluate the approximate solution at the final time $t^{n+1}$ and on the original grid. Namely, one first virtually solves the Riemann problem set at interface $x_i$ and associated with the left and right states $(\mathbf{Q}_{i-\frac{1}{2}}^{n+1/2}, \mathbf{Q}_{i+\frac{1}{2}}^{n+1/2})$, see again figure 1. We denote this solution by $(x, t) \to \hat{\mathbf{Q}}_i^{n+1}(x, t)$ up to the same space and time translations. And we evaluate the integral

$$ \mathbf{Q}_i^{n+1} = \frac{1}{\Delta x} \int_{-\frac{\Delta x}{2}}^{\frac{\Delta x}{2}} \hat{\mathbf{Q}}_i^{n+1}(x, \frac{\Delta t}{2})dx $$

using another integration of (9) on the control volume $(-\frac{\Delta x}{2}, \frac{\Delta x}{2}) \times (0, \frac{\Delta t}{2})$. In this way we arrive at a new solution at time level $n+1$

$$ \mathbf{Q}_i^{n+1} = \frac{1}{2}(\mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}} + \mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}}) - \frac{1}{2}\frac{\Delta t}{\Delta x}(\mathbf{F}(\mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}}) - \mathbf{F}(\mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}})) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{\mathbf{A}}_i(\mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}} - \mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}}) \quad (12) $$

where we set

$$ \hat{\mathbf{A}}_i = \mathbf{A}(\frac{1}{2}(\mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}} + \mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}})). $$

This concludes the description of the Price-T method as proposed in [34]. It is clear that this scheme is a centred non conservative scheme which is consistent with the non conservative system (9). However, when applied to a conservative system the coefficient matrices $\hat{\mathbf{A}}$ will disappear, leaving only the terms with the numerical flux. Thus, the Price-T scheme easily boils down to the conservative Force scheme (see Toro and Billett [31])

$$ \mathbf{Q}_i^{n+1} = \mathbf{Q}_i^n - \frac{\Delta t}{\Delta x}(\mathbf{F}_{i+1/2}^{\text{force}} - \mathbf{F}_{i-1/2}^{\text{force}}), $$

where the numerical flux $\mathbf{F}_{i+1/2}^{\text{force}} = \mathbf{F}_{i+1/2}^{\text{force}}(\mathbf{Q}_i^n, \mathbf{Q}_{i+1}^n)$ is usually written as the arithmetic average of the Lax-Friedrichs and Lax-Wendroff fluxes, namely with clear notations

$$ \mathbf{F}_{i+1/2}^{\text{force}} = \frac{1}{2}(\mathbf{F}_{i+1/2}^{\text{lf}} + \mathbf{F}_{i+1/2}^{\text{lw}}) $$

with

$$ \mathbf{F}_{i+1/2}^{\text{lf}} = \frac{1}{2}(\mathbf{F}(\mathbf{Q}_i^n) + \mathbf{F}(\mathbf{Q}_{i+1}^n)) - \frac{1}{2}\frac{\Delta x}{\Delta t}(\mathbf{Q}_{i+1}^n - \mathbf{Q}_i^n) $$

and

$$ \mathbf{F}_{i+1/2}^{\text{lw}} = \mathbf{F}(\mathbf{Q}_{i+1/2}^{n+1/2}). $$
## 2.2 MUSCL-Hancock Price-T scheme

In this subsection, we present a second-order extension of the Price-T method using a MUSCL-Hancock strategy. As already pointed out before, we will pay a particular attention to the conservation property of the scheme when applied to a conservative system ($A(Q) = 0$). Actually, our objective in this case will be to recover the Slope Limiter Centred (SLIC) method with Force flux, see for instance [38] for more details. As we will see, a correction term in a "naive" (but natural) extension of the Price-T method will be necessary.

Before going further, let us first briefly recall that the original MUSCL-Hancock strategy uses reconstruction of data together with upwind Godunov or Godunov-type methods. Instead, the aim of the SLIC scheme is to avoid the use of Riemann problems by exploiting different low-order schemes, for instance the centred Price-T flux (or Force flux in the conservative case). For more details, we refer again the reader to [38].

The different steps of the method are as follows.

**Step 1: reconstruction and cell-boundary values.** As it is customary, this step aims at defining boundary extrapolated values by means of polynomial data reconstructions. More precisely, using for each cell $I_i$ a reconstructed polynomial vector $P_i(x) = (p_{1,i}(x), \dots, p_{m,i}(x))^t$ with components

$$p_{k,i}(x) = Q_{k,i}^n + (x - x_i)\Delta_{k,i}, \quad k = 1, \dots, m$$

where $\Delta_{k,i}$ are the slopes (recall that $m$ is the number of variables), we define the boundary values by

$$ \left\{ \begin{array}{l} Q_i^L = P_i(x_{i-\frac{1}{2}}) = Q_i^n - \Delta_i \frac{\Delta x}{2} \\ Q_i^R = P_i(x_{i+\frac{1}{2}}) = Q_i^n + \Delta_i \frac{\Delta x}{2}. \end{array} \right. \qquad (13) $$

Regarding the definition of the slope, a first choice is given by the ENO strategy [38]

$$ \Delta_{k,i} = \begin{cases} \frac{\Delta_{k,i-1/2}}{\Delta x} = \frac{q_{k,i-1} - q_{k,i-1}}{\Delta x} & \text{if } |\Delta_{k,i-1/2}| < |\Delta_{k,i+1/2}| \\ \frac{\Delta_{k,i+1/2}}{\Delta x} = \frac{q_{k,i+1} - q_{k,i}}{\Delta x} & \text{otherwise} \end{cases} $$

where $k = 1, \dots, m$ refers to the unknown variables of $Q$. Other possibilities are given by the well-known Van Leer or Super Bee limiters, namely

$$ \Delta_{k,i} = \phi\left(\frac{q_{k,i} - q_{k,i-1}}{q_{k,i+1} - q_{k,i}}\right) \frac{q_{k,i+1} - q_{k,i}}{\Delta x} $$

with

$$ \phi(\theta) = \frac{\theta + |\theta|}{1 + |\theta|} $$

for the Van Leer limiter and

$$ \phi(\theta) = \max(\min(2\theta, 1), \min(\theta, 2))^{+} $$

for the Super Bee limiter. See [34] for more details.

**Step 2: time evolution of the states $Q_i^L$ and $Q_i^R$ for a time $\Delta t/2$**. In this step, we propose to evolve in time the boundary values $Q_i^L$ and $Q_i^R$ for a time $\Delta t/2$ according to

$$ \left\{ \begin{array}{l} \bar{Q}_i^L = Q_i^L - \frac{1}{2} \frac{\Delta t}{\Delta x} (F(Q_i^R) - F(Q_i^L)) - \frac{1}{2} \frac{\Delta t}{\Delta x} \hat{A}_i (Q_i^R - Q_i^L), \\ \bar{Q}_i^R = Q_i^R - \frac{1}{2} \frac{\Delta t}{\Delta x} (F(Q_i^R) - F(Q_i^L)) - \frac{1}{2} \frac{\Delta t}{\Delta x} \hat{A}_i (Q_i^R - Q_i^L). \end{array} \right. \qquad (14) $$

where $\hat{A}_i = A(Q_i^n)$. These formulas mimic both the classical MUSCL-Hancock method for the conservative part and the Price-T scheme recalled in the previous section for the non conservative one.

**Step 3: intermediate states at interfaces.** In order to compute the intermediate states $Q_{i \pm 1/2}^{n+1/2}$ at interfaces, we now exploit the Price-T formulas (10) and (11) with the evoluted values $(\bar{Q}_{i-1}^R, \bar{Q}_i^L)$ and $(\bar{Q}_i^R, \bar{Q}_{i+1}^L)$. Therefore we obtain

$$ \bar{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}} = \frac{1}{2}(\bar{Q}_{i-1}^R + \bar{Q}_i^L) - \frac{1}{2} \frac{\Delta t}{\Delta x} (F(\bar{Q}_i^L) - F(\bar{Q}_{i-1}^R)) - \frac{1}{2} \frac{\Delta t}{\Delta x} \hat{A}_{i-\frac{1}{2}} (\bar{Q}_i^L - \bar{Q}_{i-1}^R) \qquad (15) $$
and

$$ \mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}} = \frac{1}{2}(\mathbf{Q}_i^R + \bar{\mathbf{Q}}_{i+1}) - \frac{1}{2}\frac{\Delta t}{\Delta x}(\mathbf{F}(\mathbf{Q}_{i+1}) - \mathbf{F}(\bar{\mathbf{Q}}_i)) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{\mathbf{A}}_{i+\frac{1}{2}}(\bar{\mathbf{Q}}_{i+1} - \bar{\mathbf{Q}}_i^R), \quad (16) $$

where we set

$$ \hat{\mathbf{A}}_{i-\frac{1}{2}} = \mathbf{A}(\frac{1}{2}(\bar{\mathbf{Q}}_{i-1}^R + \bar{\mathbf{Q}}_i^L)) \quad \text{and} \quad \hat{\mathbf{A}}_{i+\frac{1}{2}} = \mathbf{A}(\frac{1}{2}(\bar{\mathbf{Q}}_{i}^R + \bar{\mathbf{Q}}_{i+1}^L)). $$

Step 4: new solution. At this stage, it would be tempting to update the solution using formula (12) with the new intermediate states, leading to

$$ \mathbf{Q}_{i}^{n+1} = \frac{1}{2}(\mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}} + \mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}}) - \frac{1}{2}\frac{\Delta t}{\Delta x}(\mathbf{F}(\mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}}) - \mathbf{F}(\mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}})) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{\mathbf{A}}_{i}(\mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}} - \mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}}), \quad (17) $$

with

$$ \hat{\mathbf{A}}_i = \mathbf{A}(\frac{1}{2}(\mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}} + \mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}})). $$

However, it was observed that in the conservative case $\mathbf{A}(\mathbf{Q}) = 0$ such a formula does not provide us with a finite volume conservative scheme, and in particular it does not reduce to the expected SLIC scheme associated with the Price-T flux, namely to

$$ \mathbf{Q}_{i}^{n+1} = \mathbf{Q}_{i}^{n} - \frac{\Delta t}{\Delta x}(\mathbf{F}_{i+1/2}^{slic} - \mathbf{F}_{i-1/2}^{slic}), $$

with $\mathbf{F}_{i+1/2}^{slic} = \mathbf{F}_{i+1/2}^{force}(\bar{\mathbf{Q}}_i^R, \bar{\mathbf{Q}}_{i+1}^L)$ for all $i$. Refer to [38] for more details. This is of course not satisfactory since it is well-known that non conservative schemes used to approximate conservative systems generally lead to wrong solutions. This will be illustrated in the numerical experiments below.

In order to overcome this issue, easy calculations show that a correction term $\mathbf{L}_i$ given by

$$ \mathbf{L}_i = \frac{\mathbf{F}(\bar{\mathbf{Q}}_i^R) - \mathbf{F}(\bar{\mathbf{Q}}_i^L)}{\Delta x} - \frac{\mathbf{F}(\mathbf{Q}_i^R) - \mathbf{F}(\mathbf{Q}_i^L)}{\Delta x}, $$

and clearly approximating zero can be added in the update formula to recover the SLIC scheme in the conservative case. More precisely, we simply suggest to replace (17) with

$$ \mathbf{Q}_{i}^{n+1} = \frac{1}{2}(\mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}} + \mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}}) - \frac{1}{2}\frac{\Delta t}{\Delta x}(\mathbf{F}(\mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}}) - \mathbf{F}(\mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}})) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{\mathbf{A}}_{i}(\mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}} - \mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}}) - \frac{\Delta t}{2}\mathbf{L}_{i}, \quad (18) $$

which concludes the description of our second-order centred scheme.

## 2.3 Numerical results

This section shows the numerical results given by the proposed numerical schemes in the elastic and viscoelastic cases. At this stage, we consider academic test cases, namely three different Riemann problems with initial data

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

where we recall that $\mathbf{Q} = (A, Au)^t$ in the elastic case, $\mathbf{Q} = (A, Au, \Psi)^t$ in the viscoelastic case and $L$ is the length of the blood vessel. Note that we will always use $\Psi_L = \Psi_R = 0$. No friction forces are accounted for. A more complex test case of a network with 37 arteries will be considered in section 5.

The parameters and initial states used for the simulations are given in the following tables 1 and 2. As we will see and as far as the elastic case is concerned, test 1 presents a left shock and a right rarefaction, test 2 has two rarefactions and test 3 two shocks.

The time step $\Delta t$ is defined at each time iteration by

$$ \Delta t = \text{CFL} \frac{\Delta x}{\max_i (|u_i| + \tilde{c}_i)} $$
<table><caption>Table 1: Parameters for Riemann problems in sections 2.3 and 4.3.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td><i>p</i><sub>ext</sub></td><td>0 [Pa]</td></tr><tr><td><i>L</i></td><td>0.4 [m]</td></tr><tr><td><i>h</i><sub>0</sub></td><td>0.5 × 10<sup>-3</sup> [m]</td></tr><tr><td><i>r</i><sub>0</sub></td><td>0.01 [m]</td></tr><tr><td><i>E</i></td><td>1.2 × 10<sup>6</sup> [N/m<sup>2</sup>]</td></tr><tr><td><i>ρ</i></td><td>1050 [kg/m<sup>3</sup>]</td></tr></tbody></table>

<table><caption>Table 2: Values for area A, velocity u and final time for the proposed Riemann problems.</caption><thead><tr><th>Test</th><th><i>A</i><sub>L</sub></th><th><i>u</i><sub>L</sub></th><th><i>A</i><sub>R</sub></th><th><i>u</i><sub>R</sub></th><th>Timeout</th></tr></thead><tbody><tr><td>1</td><td><i>A</i><sub>0</sub></td><td>0</td><td>2<i>A</i><sub>0</sub></td><td>0</td><td>0.013</td></tr><tr><td>2</td><td>2<i>A</i><sub>0</sub></td><td>-1</td><td>2<i>A</i><sub>0</sub></td><td>1</td><td>0.013</td></tr><tr><td>3</td><td>2<i>A</i><sub>0</sub></td><td>1</td><td>2<i>A</i><sub>0</sub></td><td>-1</td><td>0.013</td></tr></tbody></table>

and we take CFL = 0.9. The stiff source term associated with the relaxation variable Ψ is solved exactly.

Before showing the results, we underline that the proposed relaxation approach to reformulate the viscoelastic equations as a hyperbolic system with stiff source term introduces an error of order $O(ε)$. Considering that the numerical scheme approximates the solution with a consistency error $O(Δx^r)$ where $r = 1$ for the Price-T scheme and $r = 2$ for the MUSCL-Hancock Price-T scheme, it is natural to choose $Δx$ in such a way that $ε ≈ Δx^r$. In practice, we follow the optimal choice proposed in [25] and given by

$$ε = \frac{O(1)Δx^r}{K(r)} \quad \text{with} \quad K(r) = \frac{1 - 2^{-\frac{1}{2}}}{2^{r-\frac{1}{2}} - 1} \qquad (19)$$

It is important to note that, fixed the relation time and the maximum acceptable value for $O(1)$, if we refine excessively the mesh, namely we choose a small $Δx$, then the formulation error overcomes the numerical one and consequently

$$\frac{εK(r)}{Δx^r} > O(1).$$

Since the error related to the choice of $ε$ cannot be completely removed, it represents a limit to the accuracy we can obtain when using a numerical method of order $r$ on a mesh of size $Δx$. Referring to [25], in general for the numerical computations we take $O(1) = 15$, and for each test we choose the values for the relaxation time and the mesh size in agreement with (19).

**Elastic case.** Figure 2 shows the numerical solutions for the reference area $A/A_0$ and the flux $q = Au$ in the elastic case ($Γ = 0$). We compared both first and second order method with mesh of $M = 100$ points, where the mesh size is $Δx = L/M$, against the exact solution, which was computed using an exact Riemann solver [36]. For the MUSCL-Hancock scheme, the Van Leer limiter was used. We note that both schemes are satisfying, in particular MUSCL-Hancock scheme approximates the exact solution very well as expected.

In figure 3 we exhibit the MUSCL-Hancock solution for test 1 using different values for the mesh size. In particular we take $M = [62, 125, 250, 500, 1000, 2000]$ cells and illustrate that the numerical output converges to the exact solution.

Finally and for the sake of completeness, we point out in figure 4 the importance of the correction term $L_i$. As we already observed, the MUSCL-Hancock Price-T scheme applied to a conservative system does not reduce to a conservative method when we neglect $L_i$. As a consequence, the numerical outcome could converge to a wrong solution. Indeed, when we do not add the correction term, we clearly see that the approximate solution departs from the exact one.

**Viscoelastic case.** Figure 5 shows the results obtained in the viscoelastic case with $Γ = 1$ Pa s m and $M = 200$. We used both first and second order scheme, with again Van Leer limiter for MUSCL-Hancock method. We insert also a reference solution computed with the second-order MUSCL-Hancock scheme and a mesh with $M = 2000$ cells. Consequently, we
<figure><img src="image_3.png" /><figcaption>(a) Test 1 - Left shock, right rarefaction</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Test 2 - Two rarefactions</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Test 3 - Two shocks</figcaption></figure>

Figure 2: Comparison among Price-T (blue symbol), MUSCL-Hancock (green symbol) and exact solution (red line). Reference area $A/A_0$ (left) and flux $q$ (right) for the Riemann Problems in table 2. Mesh with $M = 100$ cells, $\Gamma = 0$.
<figure><img src="image_6.png" /><figcaption>Figure 3: Converged mesh solution for the reference area $A/A_0$ (left) and the flow $q$ (right) computed with MUSCL-Hancock method (dashed line) against exact solution (red line). Test 1 of table 2, mesh with $M = [62, 125, 250, 500, 1000, 2000]$ cells, $\Gamma = 0$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: Comparison among MUSCL-Hancock (green line), MUSCL-Hancock without correction term $L_i$ (black line) and exact solution (red dashed line). Reference area $A/A_0$ (left) and flux $q$ (right) for Riemann Problems in table 2. Mesh with $M = 500$ cells, $\Gamma = 0$.</figcaption></figure>
had to use a relaxation time $\varepsilon$ of order $O(10^{-6})$ according to condition (19). It is important to note that to take a finer mesh size implies the use of a smaller relaxation time according to this relation, with the consequence of a decreased time step $\Delta t$ and an increased number of iterations. Moreover, for a fixed relaxation time, the higher is the order of accuracy of the method, the larger is the mesh size we have to take to satisfy condition (19). Thus, we could have used a finer mesh if we had only considered the Price-T scheme, as it is first-order accurate. Clearly, with a fixed mesh size, Price-T solution results to be more diffusive than the second-order one.

Moreover, we note that the first-order outcome has already reached the boundaries at time $t = 0.013s$. For this reason, we exhibit the solution for the first test case with an enlarged domain, in order to show that the transmissive boundary conditions are not interfering with the results. We take the vessel length $L = 0.6m$ and consequently we consider the initial discontinuity at gate $= 0.3m$. Thus, in figure 6 satisfying results are reported.

Finally, also for the viscoelastic case we insert a convergence test using a sequence of mesh with $M = [62, 125, 250, 500, 1000, 2000]$ cells. As initial conditions, we consider the first test case of table 2, using once again the enlarged domain and relaxation time of order $O(10^{-6})$.

### 3 Continuous sensitivity equations

In this section, we consider the continuous sensitivity equation method for computing the first-order sensitivities of our elastic and viscoelastic models (6). We refer for instance the reader to [10] and [16] for more details about this method. These equations will be exploited for uncertainty quantification in sections 4.3 and 5.

Let us briefly recall that when speaking of blood flow in arteries and veins, sensitivity analysis is a crucial topic. Indeed, acquiring experimental results can be very challenging and the values of the parameters considered in the numerical simulations can present uncertainties, not to mention the underlying modeling assumptions. It is therefore very important to understand how the solution (flow and pressure in particular) could transform when varying the parameters values.

In the following, we will denote by $a$ the uncertain parameter and system (6) will be called the state system. Our main objective is to derive the so-called sensitivity system associated with the derivative of $Q$ with respect to $a$ that will be denoted by $Q_a$. To compute the sensitivity system we first take as hypothesis the smoothness of the solution and simply differentiate the state equations with respect to $a$. It amounts to formally exchange the derivative with respect to $a$ with the ones with respect to space and time. As we will see, the obtained system is proved to be weakly hyperbolic. Then, we will consider the possibility of having discontinuous state solutions $Q$ and correct the sensitivity system accordingly by adding a Dirac source term.

**Elastic case.** Let us differentiate system (3) with respect to $a$. After some calculations and exchanging the derivatives in time and space with the ones with respect to $a$, we easily obtain the sensitivity equations

$$ \begin{cases} \partial_t A_a + \partial_x q_a = 0 \\ \partial_t q_a + \partial_x ((c^2 - u^2)A_a + 2uq_a + \gamma_a A_a^{\frac{3}{2}}) = -R_a u - R u_a. \end{cases} \quad (20) $$

Also note that differentiating (2) gives

$$ p_a(A, A_a) = p_{ext,a} + 3\rho \left( \gamma_a(\sqrt{A} - \sqrt{A_0}) + \gamma\left(\frac{A_a}{2\sqrt{A}} - \frac{A_{0a}}{2\sqrt{A_0}}\right) \right) + 3\rho_a\gamma(\sqrt{A} - \sqrt{A_0}). $$

Considering both the state equations (3) and the sensitivity equations (20) together gives a new system with four equations and four unknowns. Its Jacobian matrix is

$$ J = \begin{pmatrix} J_{el} & \mathbf{0} \\ \mathbf{R} & J_{el} \end{pmatrix} = \begin{pmatrix} 0 & 1 & 0 & 0 \\ c^2 - u^2 & 2u & 0 & 0 \\ 0 & 0 & 0 & 1 \\ \frac{1}{A} \left( \left( \frac{c^2}{2} + 2u^2 \right) A_a - 2uq_a + \frac{A c^2 \gamma_a}{\gamma} \right) & \frac{1}{A} (2q_a - 2uA_a) & c^2 - u^2 & 2u \end{pmatrix}, $$

with $J_{el}$ the Jacobian matrix of the elastic model (3). Hence, the eigenvalues of the global system are still given by $\lambda_e^{\pm}$ but now with multiplicity two. Therefore, the strict hyperbolicity is immediately lost. Actually, it is easy to see that hyperbolicity
<figure><img src="image_8.png" /><figcaption>(a) Test 1</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Test 2</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(c) Test 3</figcaption></figure>

Figure 5: Price-T (blue symbol '*' ) and MUSCL-Hancock (green symbol 'o') solution computed using a mesh with M = 200 cells, reference solution (red line) with M = 2000 cells. Reference area A/A₀ (left) and flux q (right) for Riemann Problems in table 2, Γ = 1 Pa s m.
<figure><img src="image_11.png" /><figcaption>Figure 6: Price-T (blue symbol '*' ) and MUSCL-Hancock (green symbol 'o') solution computed using a mesh with M = 200 cells, reference solution (red line) with M = 2000 cells. Reference area A/A₀ (left) and flux q (right) for Riemann Problems in table 2, Γ = 1 Pa s m.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 7: Converged mesh solution for the reference area A/A₀ (left) and flow q (right) for test 1 of table 2. Mesh with M = [62, 125, 250, 500, 1000, 2000] cells, MUSCL-Hancock method, Γ = 1 Pa s m.</figcaption></figure>

is also lost in general. Indeed, the global system turns out to be only weakly hyperbolic, which means that all eigenvalues are real but no complete set of linearly independent eigenvectors exists, or equivalently the Jacobian matrix is not R - diagonalizable. In this case the characteristic polynomial reads

$$
p(x) = (x - \lambda_e^{-})^2(x - \lambda_e^{+})^2,
$$

hence the minimal polynomial should be at most of degree 2 in order to have distinct roots.
Thus, to prove that **J** is diagonalizable, we should have

$$
(\mathbf{J} - \lambda_e^- I_4)(\mathbf{J} - \lambda_e^+ I_4) = 0
$$

where $I_4$ is the 4 × 4 identity matrix. This is equivalent to write

$$
\begin{align*}
& \begin{pmatrix} \mathbf{J}_{el} - \lambda_e^- I_2 & \mathbf{0} \\ \mathbf{R} & \mathbf{J}_{el} - \lambda_e^- I_2 \end{pmatrix} \begin{pmatrix} \mathbf{J}_{el} - \lambda_e^+ I_2 & \mathbf{0} \\ \mathbf{R} & \mathbf{J}_{el} - \lambda_e^+ I_2 \end{pmatrix} = \\
&= \begin{pmatrix} (\mathbf{J}_{el} - \lambda_e^- I_2)(\mathbf{J}_{el} - \lambda_e^+ I_2) & \mathbf{0} \\ \mathbf{R}(\mathbf{J}_{el} - \lambda_e^+ I_2) + (\mathbf{J}_{el} - \lambda_e^- I_2)\mathbf{R} & (\mathbf{J}_{el} - \lambda_e^- I_2)(\mathbf{J}_{el} - \lambda_e^+ I_2) \end{pmatrix} = 0
\end{align*}
$$
and in order to have this condition satisfied, we should impose $\mathbf{R}(\mathbf{J}_{el} - \lambda_e^+ I_2) + (\mathbf{J}_{el} - \lambda_e^- I_2)\mathbf{R} = 0$. However, computing its coefficients shows that in general they are different from zero. More precisely, we have

$$ (\mathbf{R}(\mathbf{J}_{el} - \lambda_e^+ I_2) + (\mathbf{J}_{el} - \lambda_e^- I_2)\mathbf{R})_{1,1} = \frac{1}{A}((\frac{c^2}{2} + 2u^2)A_a - 2uq_a + \frac{Ac^2\gamma_a}{\gamma}) \quad (21) $$

and

$$ (\mathbf{R}(\mathbf{J}_{el} - \lambda_e^+ I_2) + (\mathbf{J}_{el} - \lambda_e^- I_2)\mathbf{R})_{1,2} = \frac{1}{A}(2q_a - 2uA_a) = 2u_a \quad (22) $$

so that we would need $u_a = 0$ by (22) and then $A_a = 2A\gamma_a/\gamma$ by (21). Both conditions are very specific and in general are not true, which proves the weak hyperbolicity.

**Viscoelastic case.** Let us now turn to the viscoelastic case and differentiate (5) with respect to the uncertain parameter $a$. After some computations, we obtain

$$ \left\{ \begin{aligned} & \partial_t A_a + \partial_x q_a = 0 \\ & \partial_t q_a + \partial_x \left( (c^2 - u^2)A_a + 2uq_a + \gamma_a A^{\frac{3}{2}} \right) + \frac{a_{\Gamma_a}}{2} \partial_x A + \frac{a_{\Gamma}}{2} \partial_x A_a - \\ & \qquad - \zeta \left( \frac{1}{2\sqrt{A}} A_a \partial_x \Psi + \sqrt{A} \partial_x \Psi_a \right) - \zeta_a \sqrt{A} \partial_x \Psi = -R_a u - R u_a \\ & \partial_t \Psi_a - \frac{1}{\epsilon} \partial_x q_a + \frac{1}{\epsilon^2} \epsilon_a \partial_x q = - \frac{1}{\epsilon} \Psi_a + \frac{1}{\epsilon^2} \epsilon_a \Psi, \end{aligned} \right. \quad (23) $$

with $\zeta = \frac{\Gamma}{\rho A_0}$. Moreover, the derivative of the viscoelastic pressure reads

$$ \begin{aligned} p_a(A, \Psi, A_a, \Psi_a) = & p_{ext,a} + 3\rho \left( \gamma_a(\sqrt{A} - \sqrt{A_0}) + \gamma_a \left( \frac{A_a}{2\sqrt{A}} - \frac{A_{0a}}{2\sqrt{A_0}} \right) \right) + 3\rho_a \gamma(\sqrt{A} - \sqrt{A_0}) - \\ & - \Gamma_a \frac{\Psi}{A_0 \sqrt{A}} + \Gamma \frac{\Psi A_{0a}}{A_0^2 \sqrt{A}} + \Gamma \frac{\Psi A_a}{2A_0 \sqrt{A^3}} - \Gamma \frac{\Psi_a}{A_0 \sqrt{A}}. \end{aligned} \quad (24) $$

In the same way we did for the elastic system, it can be easily shown that the global system (5)-(23) is weakly hyperbolic with real eigenvalues given by $\lambda_v^\pm$ and $\lambda_v^0$ and multiplicity two.

For the sake of conciseness, we rewrite the global systems (3)-(20) for the elastic case and (5)-(23) for the viscoelastic case in the compact form

$$ \left\{ \begin{aligned} & \partial_t \mathbf{Q} + \partial_x \mathbf{F}(\mathbf{Q}) + \mathbf{A}(\mathbf{Q}) \partial_x \mathbf{Q} = \mathbf{S}(\mathbf{Q}) \\ & \partial_t \mathbf{Q}_a + \partial_x \mathbf{F}_a(\mathbf{Q}, \mathbf{Q}_a) + \mathbf{B}(\mathbf{Q}, \mathbf{Q}_a) \partial_x \mathbf{Q} + \mathbf{A}(\mathbf{Q}) \partial_x \mathbf{Q}_a = \mathbf{S}_a(\mathbf{Q}, \mathbf{Q}_a) \end{aligned} \right. \quad (25) $$

where in the new equation on $\mathbf{Q}_a$, we have $\mathbf{B}(\mathbf{Q}, \mathbf{Q}_a) = 0$ and

$$ \mathbf{Q}_a = \begin{pmatrix} A_a \\ q_a \\ \Psi_a \end{pmatrix}, \quad \mathbf{F}_a(\mathbf{Q}) = \begin{pmatrix} q_a \\ (c^2 - u^2)A_a + 2uq_a + \gamma_a A^{\frac{3}{2}} \\ -\frac{1}{\epsilon}q_a + \frac{\epsilon_a}{\epsilon^2}q \end{pmatrix}, \quad \mathbf{S}_a(\mathbf{Q}, \mathbf{Q}_a) = \begin{pmatrix} 0 \\ -R_a u - R u_a \\ -\frac{1}{\epsilon} \Psi_a + \frac{1}{\epsilon^2} \epsilon_a \Psi \end{pmatrix} \quad (26) $$

for the elastic model, and

$$ \begin{aligned} \mathbf{Q}_a &= \begin{pmatrix} A_a \\ q_a \\ \Psi_a \end{pmatrix}, \quad \mathbf{F}_a(\mathbf{Q}) = \begin{pmatrix} q_a \\ (c^2 - u^2)A_a + 2uq_a + \gamma_a A^{\frac{3}{2}} \\ -\frac{1}{\epsilon}q_a + \frac{\epsilon_a}{\epsilon^2}q \end{pmatrix}, \\ \mathbf{B} &= \partial_a \mathbf{A}(\mathbf{Q}) = \begin{pmatrix} 0 & 0 & 0 \\ \frac{a_{\Gamma_a}}{2} & 0 & -\zeta \frac{A_a}{2\sqrt{A}} - \zeta_a \sqrt{A} \\ 0 & 0 & 0 \end{pmatrix}, \quad \mathbf{S}_a(\mathbf{Q}, \mathbf{Q}_a) = \begin{pmatrix} 0 \\ -R_a u - R u_a \\ -\frac{1}{\epsilon} \Psi_a + \frac{1}{\epsilon^2} \epsilon_a \Psi \end{pmatrix} \end{aligned} \quad (27) $$
for the viscoelastic one. Note that the global viscoelastic system has two additional terms due to the non-conservative part of the viscoelastic state system.

**The case of discontinuous solutions.** So far, it was implicitly assumed that the solutions had sufficient regularity to differentiate with respect to $a$ and exchange the derivatives. However and at least in the elastic case, $Q$ can be discontinuous in practice (see for instance the numerical simulations above). In order to avoid the appearance of Dirac delta functions in the sensitivity solution $Q_a$ whenever $Q$ is not continuous, it was proposed in [8] and [10] to add a compensation term $\hat{S}(Q)$ of the form

$$ \hat{S}(Q) = \sum_{k=1}^{N_s} \delta_k \rho_k \qquad (28) $$

where $N_s$ is the number of discontinuities,

$$ \delta_k = \delta(x - x_{k,s}) $$

is the Dirac delta function centered in the position $x_{s,k}$ of the $k$-th discontinuity separating the states $Q^-$ and $Q^+$, and $\rho_k$ is the amplitude of the correction for the $k$-th discontinuity given by

$$ \rho_k(t) = \sigma_{k,a}(Q^+ - Q^-). \qquad (29) $$

We refer for instance the reader to [17] and [10] and the references therein for more details. Finally, the compact forms of the global systems are now given with clear notations by

$$ \left\{ \begin{aligned} \partial_t Q + \partial_x F(Q) + A(Q)\partial_x Q &= S(Q) \\ \partial_t Q + \partial_x F_a(Q, Q_a) + B(Q, Q_a)\partial_x Q + A(Q)\partial_x Q_a &= S_a(Q, Q_a) + \hat{S}(Q). \end{aligned} \right. \qquad (30) $$

# 4 Centred numerical scheme for the continuous sensitivity equations

The aim of this section is to briefly describe the Price-T and MUSCL-Hancock schemes applied to the sensitivity variables $Q_a$. The main difference with the general setting proposed in section 2 lies in the presence of the source term $\hat{S}(Q)$ in the sensitivity equations. Note that the friction forces or the stiff source term associated with the relaxation variable $\Psi$ are considered using the same operator splitting method as for the state equations on $Q$. On the contrary, we handle $\hat{S}$ in a different way and take it into account at the same level as for the convective terms.

## 4.1 First-order Price-T scheme

Since we consider also the source term $\hat{S}(Q)$ for the sensitivity variables, we modify the updating formulae (10), (11) and (12) in the following way:

$$ \begin{aligned} Q_{a,i-\frac{1}{2}}^{n+\frac{1}{2}} = & \frac{1}{2}(Q_{a,i-1}^n + Q_{a,i}^n) - \frac{1}{2}\frac{\Delta t}{\Delta x}(F_a(Q_i^n, Q_{a,i}^n) - F_a(Q_{i-1}^n, Q_{a,i-1}^n)) - \\ & - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{B}_{i-\frac{1}{2}}(Q_i^n - Q_{i-1}^n) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{A}_{i-\frac{1}{2}}(Q_{a,i}^n - Q_{a,i-1}^n) + \frac{\Delta t}{2}\hat{S}_{i-\frac{1}{2}}^n, \end{aligned} $$

$$ \begin{aligned} Q_{a,i+\frac{1}{2}}^{n+\frac{1}{2}} = & \frac{1}{2}(Q_{a,i}^n + Q_{a,i+1}^n) - \frac{1}{2}\frac{\Delta t}{\Delta x}(F_a(Q_{i+1}^n, Q_{a,i+1}^n) - F_a(Q_i^n, Q_{a,i}^n)) - \\ & - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{B}_{i+\frac{1}{2}}(Q_{i+1}^n - Q_i^n) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{A}_{i+\frac{1}{2}}(Q_{a,i+1}^n - Q_{a,i}^n) + \frac{\Delta t}{2}\hat{S}_{i+\frac{1}{2}}^n, \end{aligned} $$

$$ \begin{aligned} Q_{a,i}^{n+1} = & \frac{1}{2}(Q_{a,i-\frac{1}{2}}^{n+\frac{1}{2}} + Q_{a,i+\frac{1}{2}}^{n+\frac{1}{2}}) - \frac{1}{2}\frac{\Delta t}{\Delta x}(F_a(Q_{i+\frac{1}{2}}^{n+\frac{1}{2}}, Q_{a,i+\frac{1}{2}}^{n+\frac{1}{2}}) - F_a(Q_{i-\frac{1}{2}}^{n+\frac{1}{2}}, Q_{a,i-\frac{1}{2}}^{n+\frac{1}{2}})) - \\ & - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{B}_i(Q_{i+\frac{1}{2}}^{n+\frac{1}{2}} - Q_{i-\frac{1}{2}}^{n+\frac{1}{2}}) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{A}_i(Q_{a,i+\frac{1}{2}}^{n+\frac{1}{2}} - Q_{a,i-\frac{1}{2}}^{n+\frac{1}{2}}) + \frac{\Delta t}{2}\hat{S}_i^{n+\frac{1}{2}}, \end{aligned} $$
where $\hat{\mathbf{S}}_{i-\frac{1}{2}}^n$ is the source term referred to the cells $I_{i-1}$ and $I_i$, while $\hat{\mathbf{S}}_i^{n+\frac{1}{2}}$ to $I_i$. In agreement with formulae (28) and (29), we suggest the following formulas,

$$ 
\begin{aligned}
\hat{\mathbf{S}}_{i-\frac{1}{2}}^n &= \frac{1}{\Delta x} \left( \sigma_{1a,i-\frac{1}{2}} (\mathbf{Q}_{i-\frac{1}{2}}^* - \mathbf{Q}_{i-1}^n) \delta_{1,i-\frac{1}{2}} + \sigma_{2a,i-\frac{1}{2}} (\mathbf{Q}_i^n - \mathbf{Q}_{i-\frac{1}{2}}^*) \delta_{2,i-\frac{1}{2}} \right), \\
\hat{\mathbf{S}}_i^{n+\frac{1}{2}} &= \frac{1}{\Delta x} \left( \sigma_{1a,i} (\mathbf{Q}_i^{*,n+\frac{1}{2}} - \mathbf{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}}) \delta_{1,i} + \sigma_{2a,i} (\mathbf{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}} - \mathbf{Q}_i^{*,n+\frac{1}{2}}) \delta_{2,i} \right),
\end{aligned}
\quad (31) $$

where for $k = 1, 2$, $\delta_{k,i+1/2} = \delta_{k,i+1/2}(\mathbf{Q}_i^n, \mathbf{Q}_{i+1}^n)$, $\delta_{k,i} = \delta_{k,i}(\mathbf{Q}_{i-1/2}^{n+1/2}, \mathbf{Q}_{i+1/2}^{n+1/2})$ are shock detectors, while $\sigma_{ka}(\mathbf{Q}_i^n, \mathbf{Q}_{i+1}^n)$, $\sigma_{ka,i} = \sigma_{ka}(\mathbf{Q}_{i-1/2}^{n+1/2}, \mathbf{Q}_{i+1/2}^{n+1/2})$ and $\mathbf{Q}_{i+1/2}^* = \mathbf{Q}^*(\mathbf{Q}_i^n, \mathbf{Q}_{i+1}^n)$, $\mathbf{Q}_i^* = \mathbf{Q}^*(\mathbf{Q}_{i-1/2}^{n+1/2}, \mathbf{Q}_{i+1/2}^{n+1/2})$ respectively approximate the derivative with respect to $a$ of the speed of propagation $\sigma_k$ of the corresponding shock, and its left or right state.

In order to evaluate these quantities and since in general the exact solution of a Riemann problem is not explicitly known (at least for the viscoelastic system), we use the Harten-Lax-van Leer formalism applied to non conservative systems (see [38] and [2]) to define a HLL approximate Riemann solver. More precisely, if we denote $\mathbf{Q}_L$ and $\mathbf{Q}_R$ the left and right Riemann initial states, we define $\sigma_k = \sigma_k(\mathbf{Q}_L, \mathbf{Q}_R)$ and $\mathbf{Q}^* = \mathbf{Q}^*(\mathbf{Q}_L, \mathbf{Q}_R)$ using the consistency relation

$$ \Delta \mathbf{F} + \hat{\mathbf{A}} \Delta \mathbf{Q} = \sigma_1(\mathbf{Q}^* - \mathbf{Q}_L) + \sigma_2(\mathbf{Q}_R - \mathbf{Q}^*), $$

leading to

$$ \mathbf{Q}^* = \frac{\Delta \mathbf{F} + \hat{\mathbf{A}} \Delta \mathbf{Q} + \mathbf{Q}_L \sigma_1 - \mathbf{Q}_R \sigma_2}{\sigma_1 - \sigma_2}, $$

where we have used the classical notations $\Delta \mathbf{F} = \mathbf{F}(\mathbf{Q}_R) - \mathbf{F}(\mathbf{Q}_L)$, $\Delta \mathbf{Q} = \mathbf{Q}_R - \mathbf{Q}_L$ and $\hat{\mathbf{A}}$ is such that the consistency condition

$$ \lim_{\mathbf{Q}_L, \mathbf{Q}_R \to \mathbf{Q}} \hat{\mathbf{A}}(\mathbf{Q}_L, \mathbf{Q}_R) = \mathbf{A}(\mathbf{Q}) $$

holds true. We simply suggest

$$ \hat{\mathbf{A}} = \mathbf{A}(\frac{1}{2}(\mathbf{Q}_L + \mathbf{Q}_R)). $$

Before defining the speeds of propagation and the shock detectors, observe that using a HLL approximate Riemann solver is the more natural way for the elastic system since the exact Riemann solution itself contains exactly two waves. Therefore, we are not excluding any of them and the approach results to be complete (but approximate). On the other hand, for systems with three or more waves we would neglect the middle ones. To avoid this problem, one could determine the intermediate values exploiting an HLLC approach, refer for instance to [35, 10].

At last, for the elastic model, we propose to define the speed of propagation as follows,

$$ \sigma_1 = \min(u_L - c_L, u_R - c_R), \quad \sigma_2 = \max(u_L + c_L, u_R + c_R), $$

refer to [13]. Therefore, the derivatives $\sigma_{1a} = \partial_a \sigma_1$ and $\sigma_{2a} = \partial_a \sigma_2$ will be

$$ \sigma_{1a} = \begin{cases} u_{aL} - c_{aL} & \text{if } \sigma_1 = u_L - c_L \\ u_{aR} - c_{aR} & \text{if } \sigma_1 = u_R - c_R \end{cases}, \quad \sigma_{2a} = \begin{cases} u_{aL} + c_{aL} & \text{if } \sigma_2 = u_L + c_L \\ u_{aR} + c_{aR} & \text{if } \sigma_2 = u_R + c_R \end{cases}, $$

with

$$ c_a = \frac{1}{2c} \left( \frac{3}{2} \gamma_a \sqrt{A} + \frac{3\gamma}{4\sqrt{A}} A_a \right). $$

Regarding the shock detectors, we refer for instance to [11] and we set

$$ \delta_{1,i-\frac{1}{2}} = \begin{cases} 1 & \text{if } A_{i-\frac{1}{2}}^* > A_{i-1}^n \\ 0 & \text{if } A_{i-\frac{1}{2}}^* \le A_{i-1}^n \end{cases}, \quad \delta_{2,i-\frac{1}{2}} = \begin{cases} 1 & \text{if } A_{i-\frac{1}{2}}^* > A_i^n \\ 0 & \text{if } A_{i-\frac{1}{2}}^* \le A_i^n \end{cases}, $$

$$ \delta_{1,i} = \begin{cases} 1 & \text{if } A_i^{*,n+\frac{1}{2}} > A_{i-\frac{1}{2}}^{n+\frac{1}{2}}, \\ 0 & \text{if } A_i^{*,n+\frac{1}{2}} \le A_{i-\frac{1}{2}}^{n+\frac{1}{2}}, \end{cases}, \quad \delta_{2,i} = \begin{cases} 1 & \text{if } A_i^{*,n+\frac{1}{2}} > A_{i+\frac{1}{2}}^{n+\frac{1}{2}}, \\ 0 & \text{if } A_i^{*,n+\frac{1}{2}} \le A_{i+\frac{1}{2}}^{n+\frac{1}{2}}. \end{cases} $$
## 4.2 MUSCL-Hancock Price-T scheme

Proceeding in a similar way, we revise the MUSCL-Hancock scheme proposed in section 2.2 to include the term $\hat{S}$ for the sensitivity equations. First of all we have to modify the evolved boundary values (14), namely

$$ \bar{Q}_{a,i}^{L} = Q_{a,i}^{L} - \frac{1}{2} \frac{\Delta t}{\Delta x} (F_a(Q_i^R, Q_{a,i}^R) - F_a(Q_i^L, Q_{a,i}^L)) - \frac{1}{2} \frac{\Delta t}{\Delta x} \hat{A}_i^n (Q_{a,i}^R - Q_{a,i}^L) - \frac{1}{2} \frac{\Delta t}{\Delta x} \hat{B}_i^n (Q_i^R - Q_i^L) + \frac{\Delta t}{2} \hat{S}_i, $$

$$ \bar{Q}_{a,i}^{R} = Q_{a,i}^{R} - \frac{1}{2} \frac{\Delta t}{\Delta x} (F_a(Q_i^R, Q_{a,i}^R) - F_a(Q_i^L, Q_{a,i}^L)) - \frac{1}{2} \frac{\Delta t}{\Delta x} \hat{A}_i^n (Q_{a,i}^R - Q_{a,i}^L) - \frac{1}{2} \frac{\Delta t}{\Delta x} \hat{B}_i^n (Q_i^R - Q_i^L) + \frac{\Delta t}{2} \hat{S}_i, $$

with $\hat{S}_i$ evaluated as in (31) with left and right states $Q_i^L$ and $Q_i^R$. Then, as we did for the Price-T scheme, we have to consider the source term $\hat{S}$ for both the intermediate states $Q_{a,i \pm \frac{1}{2}}$ and the new solution $Q_i^{n+1}$, leading to

$$ \begin{aligned} \bar{Q}_{a,i-\frac{1}{2}}^{n+\frac{1}{2}} = & \frac{1}{2}(\bar{Q}_{a,i-1}^R + \bar{Q}_{a,i}^L) - \frac{1}{2}\frac{\Delta t}{\Delta x}(F_a(\bar{Q}_i^L, \bar{Q}_{a,i}^L) - F_a(\bar{Q}_{i-1}^R, \bar{Q}_{a,i-1}^R)) - \\ & - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{B}_{i-\frac{1}{2}}(\bar{Q}_i^L - \bar{Q}_{i-1}^R) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{A}_{i-\frac{1}{2}}(\bar{Q}_{a,i}^L - \bar{Q}_{a,i-1}^R) + \frac{\Delta t}{2}\hat{S}_{i-\frac{1}{2}}^n, \end{aligned} $$

$$ \begin{aligned} \bar{Q}_{a,i+\frac{1}{2}}^{n+\frac{1}{2}} = & \frac{1}{2}(\bar{Q}_{a,i}^R + \bar{Q}_{a,i+1}^L) - \frac{1}{2}\frac{\Delta t}{\Delta x}(F_a(\bar{Q}_{i+1}^L, \bar{Q}_{a,i+1}^L) - F_a(\bar{Q}_i^R, \bar{Q}_{a,i}^R)) - \\ & - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{B}_{i+\frac{1}{2}}(\bar{Q}_{i+1}^L - \bar{Q}_i^R) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{A}_{i+\frac{1}{2}}(\bar{Q}_{a,i+1}^L - \bar{Q}_{a,i}^R) + \frac{\Delta t}{2}\hat{S}_{i+\frac{1}{2}}^n, \end{aligned} $$

$$ \begin{aligned} \bar{Q}_{a,i}^{n+1} = & \frac{1}{2}(\bar{Q}_{a,i-\frac{1}{2}}^{n+\frac{1}{2}} + \bar{Q}_{a,i+\frac{1}{2}}^{n+\frac{1}{2}}) - \frac{1}{2}\frac{\Delta t}{\Delta x}(F_a(\bar{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}}, \bar{Q}_{a,i+\frac{1}{2}}^{n+\frac{1}{2}}) - F_a(\bar{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}}, \bar{Q}_{a,i-\frac{1}{2}}^{n+\frac{1}{2}})) - \\ & - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{B}_i(\bar{Q}_{i+\frac{1}{2}}^{n+\frac{1}{2}} - \bar{Q}_{i-\frac{1}{2}}^{n+\frac{1}{2}}) - \frac{1}{2}\frac{\Delta t}{\Delta x}\hat{A}_i(\bar{Q}_{a,i+\frac{1}{2}}^{n+\frac{1}{2}} - \bar{Q}_{a,i-\frac{1}{2}}^{n+\frac{1}{2}}) - \frac{\Delta t}{2}\mathbf{L}_i + \frac{1}{2}\Delta t \hat{S}_i^{n+\frac{1}{2}}, \end{aligned} $$

with

$$ \mathbf{L}_i = \frac{F_a(\bar{Q}_i^R, \bar{Q}_{a,i}^R) - F_a(\bar{Q}_i^L, \bar{Q}_{a,i}^L)}{\Delta x} - \frac{F_a(\bar{Q}_i^R, \bar{Q}_{a,i}^R) - F_a(\bar{Q}_i^L, \bar{Q}_{a,i}^L)}{\Delta x} $$

and $\hat{S}$ computed as in (31).

## 4.3 Numerical results and uncertainty quantification

This section shows the outcomes of the proposed numerical schemes applied to the sensitivity equations. In particular, our aim is to illustrate the behavior of the solutions with respect to the emergence or disappearance of Dirac peaks when the source term $\hat{S}$ is active or not. Despite the test cases considered here are academic, we believe that they provide a better understanding of the influence of the source term $\hat{S}$ and its numerical approximation. As already said, a more realistic test case of a complex network will be considered in section 5.

As far as the state variables are concerned, we consider the same Riemann problems as in section 2.3, with values listed in table 2. Regarding the sensitivity variables and choosing the uncertain parameter $a$, we will simply have

$$ \mathbf{Q}_{aL} = \partial_a \mathbf{Q}_L, \quad \mathbf{Q}_{aR} = \partial_a \mathbf{Q}_R. $$

In particular, we study the sensitivity of the solution with respect to the initial values of the cross-sectional area, namely we take $a = A_L$ and $a = A_R$. Furthermore, as we mentioned in the introduction, an important parameter in this field is the arterial stiffness, as it associated to many cardiovascular illnesses. Thus, we consider as uncertain parameters also the vessel thickness $h_0$, the radius $r_0$ of the cross-sectional area at equilibrium and the Young modulus $E$, as they are all related to the arterial stiffness.

Therefore, using for instance $a = A_R$, we obtain

$$ \mathbf{Q}_{aL} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} \quad \text{and} \quad \mathbf{Q}_{aR} = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}. $$
<figure><img src="image_13.png" /><figcaption>Figure 8: Comparison among Price-T (blue symbol '*'), MUSCL-Hancock (green symbol 'o') and exact solution (red line). Derivative of the area $A_α$ and flux $q_α$ with $a = A_L$ for test 1 of table 2. Mesh with $M = 200$ cells, $Γ = 0$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 9: Comparison among Price-T (blue symbol '*'), MUSCL-Hancock (green symbol 'o') and exact solution (red line). Derivative of the area $A_α$ and flux $q_α$ with $a = A_R$ for test 1 of table 2. Mesh with $M = 200$ cells, $Γ = 0$.</figcaption></figure>

Observe that if we take $a \neq A_R$ and $a \neq A_L$, $Q_{aL}$ and $Q_{aR}$ will be zero vectors.

If not specified, we underline that we used a CFL condition of $CFL = 0.9$.

**Elastic Case.** In all the graphics we insert the exact solution, which was computed using an exact Riemann solver, refer to Chifari's thesis [11]. Figures 8 - 12 correspond to test 1 of table 2 and show the outputs with and without the numerical source term $\hat{S}$ for the different choices of $a$. As expected, we first observe that neglecting the correction term $\hat{S}$ implies the presence of spikes at the position of the discontinuities of the state variables. For $a = A_L$ and $a = A_R$, the presence of $\hat{S}$ removes completely the spikes and gives the sought solution. Instead, for $a = h_0$, $a = r_0$ and $a = E$, we note very little spikes on the sensitivity variable $A_α$ in correspondence of the discontinuity on $A$. In that case, our discretization of $\hat{S}$ appears to over-compensate the initial spikes, in the sense that it is now in the opposite direction with respect to the figures in which we neglected $\hat{S}$.

Next, we aim to illustrate empirically that the MUSCL-Hancock solution converges to the exact one when we consider the correction term $\hat{S}$ and we refine the mesh at the same time. Therefore, we use a mesh with $M = 2000$ cells to compare the two different outcomes. Indeed, in figure 13 we can see that the MUSCL-Hancock approximation is very satisfying, only in graphics 13.(c), 13.(d) and 13.(e) small spikes remain in correspondence of the shock of the cross-sectional area $A$. Moreover, there is a minor clipping of extreme values which is expected to further reduce refining the mesh.

Figure 14 shows the results for test 2 of table 2. Since in this Riemann problem there are two rarefactions, there is no need to insert the source term $\hat{S}$ in the sensitivity equations, as it would give the same solution. Let us also observe that in figures 14.c, 14.d and 14.e, there is a clipping of extreme values.
<figure><img src="image_15.png" /><figcaption>Figure 10: Comparison among Price-T (blue symbol '*'), MUSCL-Hancock (green symbol 'o') and exact solution (red line). Derivative of the area $A_a$ and flux $q_a$ with $a = h_0$ for test 1 of table 2. Mesh with $M = 200$ cells, $\Gamma = 0$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 11: Comparison among Price-T (blue symbol '*'), MUSCL-Hancock (green symbol 'o') and exact solution (red line). Derivative of the area $A_a$ and flux $q_a$ with $a = r_0$ for test 1 of table 2. Mesh with $M = 200$ cells, $\Gamma = 0$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 12: Comparison among Price-T (blue symbol '*'), MUSCL-Hancock (green symbol 'o') and exact solution (red line). Derivative of the area $A_a$ and flux $q_a$ with $a = E$ for test 1 of table 2. Mesh with $M = 200$ cells, $\Gamma = 0$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 13: Comparison between MUSCL-Hancock coupled with correction term $\hat{S}$ (green symbol 'o') and exact solution (red line). Derivative of the area $A_a$ and flux $q_a$ for test 1 of table 2. Mesh with $M = 2000$ cells, $\Gamma = 0$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 14: Comparison between MUSCL-Hancock (green symbol) and exact solution (red line). Derivative of area $A_a$ and flux $q_a$ for test 2 of table 2. Mesh with $M = 200$ cells, $\Gamma = 0$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 15: Comparison between MUSCL-Hancock (green symbol 'o') and exact solution (red line). Derivative of the area $A_a$ and flux $q_a$ with $a = A_L$ for test 3 of table 2. Mesh with $M = 200$ cells, $\Gamma = 0$, CFL = 0.5.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 16: Comparison between MUSCL-Hancock (green symbol 'o') and exact solution (red line). Derivative of the area $A_a$ and flux $q_a$ with $a = A_R$ for test 3 of table 2. Mesh with $M = 200$ cells, $\Gamma = 0$, CFL = 0.5.</figcaption></figure>

Finally, figures 15 - 19 show the outcomes for test 3 of table 2. For this case we used CFL = 0.5 in order to be sure that the numerical solution would reach the exact plateau in the star region. Regarding the cases with $a = A_L$ and $a = A_R$, we observe that there are no more spikes with the source term $\hat{S}$. In 17.(b), 18.(b) and 19.(b), spikes are much smaller and compared to test 1, they go in the same direction as without $\hat{S}$.

*Viscoelastic Case.* Let us first recall that in this case there is no need to add the source term $\hat{S}$ to the sensitivity equations as the state variables do not present discontinuities. In all the graphics we insert a reference solution obtained with the MUSCL-Hancock method and a mesh of $M = 2000$ cells.

Outcomes of test 1 of table 2 for the derivative of the cross-sectional area $A$ and the flow $q$ are given in figure 20. For this Riemann problem we compare the Price-T and MUSCL-Hancock solution using a mesh with $M = 200$ cells. Clearly, being first and second-order accurate respectively, the Price-T outcomes are more diffusive than the MUSCL-Hancock ones. Taking $a = A_L$ and $a = A_R$ and comparing the elastic and viscoelastic outcomes (figures 8 - 9 and 20.(a) - 20.(b)), we observe that the sensitivity variables present analogous behaviours, even if in the latter case they are more diffusive. Whereas, to take into consideration $a = h_0$, $a = r_0$ and $a = E$ leads to more obvious dissimilarities, see figures 10 - 12 for the elastic case, 20.(c) - 20.(e) for the viscoelastic one. Indeed, in the viscoelastic sensitivity solutions we note low and wide "spikes" in correspondence of the shocks of the elastic state variables. They go in the same direction of the spikes of the elastic sensitivities, which are present when we neglect the source term $\hat{S}$. Whereas, concerning rarefactions, elastic and viscoelastic sensitivity variables exhibit similar behaviours.

Analogous considerations apply to the Riemann problems 2 and 3 of table 2. For these two cases we consider only
<figure><img src="image_22.png" /><figcaption>Figure 17: Comparison between MUSCL-Hancock (green symbol 'o') and exact solution (red line). Derivative of the area $A_a$ and flux $q_a$ with $a = h_0$ for test 3 of table 2. Mesh with $M = 200$ cells, $\Gamma = 0$, CFL = 0.5.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 18: Comparison between MUSCL-Hancock (green symbol 'o') and exact solution (red line). Derivative of the area $A_a$ and flux $q_a$ with $a = r_0$ for test 3 of table 2. Mesh with $M = 200$ cells, $\Gamma = 0$, CFL = 0.5.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 19: Comparison between MUSCL-Hancock (green symbol 'o') and exact solution (red line). Derivative of the area $A_a$ and flux $q_a$ with $a = E$ for test 3 of table 2. Mesh with $M = 200$ cells, $\Gamma = 0$, CFL = 0.5.</figcaption></figure>
MUSCL-Hancock solution with $M = 200$ cells, as Price-T outcome converges to it decreasing the mesh size. In particular, we insert the outputs for test 2 in figure 21. The viscoelastic sensitivity outcome shows a similar behaviour to the one of the elastic sensitivities, even if much more diffused. This was expected as the elastic state variables present only rarefactions in this test case, and thus the solution is regular enough.

Finally, in figure 22 the results for test 3 are reported. Being this test a two-shocks Riemann Problem, once again we note that for $a = h_0$, $a = r_0$ and $a = E$ the elastic and viscoelastic sensitivity outcomes are quite different. Indeed, when we do not consider $\bar{S}$, the big spikes present in the elastic sensitivities are in correspondence to the low and wide "spikes" we note in the viscoelastic outcomes. Whereas, for $a = A_L$ and $a = A_R$, the solutions of the two sensitivity systems are similar, even if the viscoelastic output is more diffused.

*Uncertainty quantification.* Now we want to start exploiting our work setting based on sensitivity analysis for uncertainty quantification. Recall that the main aim of uncertainty quantification is to compute confidence intervals. This could be done of course using probabilistic methods like Monte-Carlo, but at the price of a much higher computational cost compared to deterministic approach proposed here. For more details about this kind of methods, see for instance [17] and the references therein.

Given a random variable $X$, its confidence interval $CI_X$ is defined by

$$CI_X = [\mu_X - k\sigma_X, \mu_X + k\sigma_X],$$

where $\mu_X$ and $\sigma_X$ respectively denote the mean and standard deviation of $X$. The value of $k$ regulates the amplitude of the interval. For instance for Gaussian random variables, $k = 1.96$ means a 95% confidence interval. Sensitivity variables comes into play to easily compute an estimate of $\mu_X$ and $\sigma_X^2$, as we show now. Note that in practice, the random variable $X$ will be either the area $A$, or the flow $q$, or the pressure $p$, while the randomness of $X$ comes from the uncertain feature of one or several parameters $a_i$. We will denote by $M$ the number of uncertain parameters and they will be represented by the vector $\mathbf{a} = (a_1, ..., a_M)^t$.

Let $\mu_a$ and $\sigma_a$ respectively be the average and the covariance matrix of the uncertain vector $\mathbf{a}$,

$$\mu_a = \begin{pmatrix} \mu_{a_1} \\ \vdots \\ \mu_{a_M} \end{pmatrix} \quad \text{and} \quad \sigma_a = \begin{pmatrix} \sigma_{a_1}^2 & \operatorname{cov}(a_1, a_2) & \dots & \operatorname{cov}(a_1, a_M) \\ \operatorname{cov}(a_1, a_2) & \sigma_{a_2}^2 & \dots & \operatorname{cov}(a_2, a_M) \\ \vdots & \vdots & \ddots & \vdots \\ \operatorname{cov}(a_1, a_M) & \dots & \dots & \sigma_{a_M}^2 \end{pmatrix},$$

with $\mu_{a_i}$ the average of the $i$-th uncertain parameter and $\sigma_{a_i}^2$ its variance. In order to estimate the mean and standard deviation of $X$, we use a Taylor expansion for $X$ around the mean of $\mathbf{a}$, namely

$$X(\mathbf{a}) = X(\mu_{\mathbf{a}}) + \sum_{i=1}^{M} (a_i - \mu_{a_i})X_{a_i}(\mu_{a_i}) + o(||\mathbf{a}||^2),$$

where $X_{a_i}$ is the derivative of $X$ with respect to the $i$-th uncertain parameter, that is to say what we have called the sensitivity variable with respect to $a_i$. To first-order and taking the average, we get

$$\mu_X = \mathbb{E}[X(\mathbf{a})] = X(\mu_{\mathbf{a}}) + \sum_{i=1}^{M} X_{a_i}(\mu_{\mathbf{a}})\mathbb{E}[a_i - \mu_{a_i}] = X(\mu_{\mathbf{a}}).$$

Indeed, recall that $\mathbb{E}[a_i - \mu_{a_i}] = 0$ and observe that $X(\mu_{a_i})$ and $X_{a_i}(\mu_{a_i})$ are not random variables. As far as the variance is concerned, we find to first-order

$$\sigma_X^2 = \mathbb{E}[(X(\mathbf{a}) - \mu_X)^2] = \mathbb{E}\left[(\sum_{i=1}^{M} X_{a_i}(\mu_{\mathbf{a}})(a_i - \mu_{a_i}))^2\right] = \\ = \sum_{i=1}^{M} X_{a_i}^2(\mu_{\mathbf{a}})\mathbb{E}[(a_i - \mu_{a_i})^2] + \sum_{\substack{i,j=1 \\ i \neq j}}^{M} X_{a_i}(\mu_{\mathbf{a}})X_{a_j}(\mu_{\mathbf{a}})\mathbb{E}[(a_i - \mu_{a_i})(a_j - \mu_{a_j})].$$
<figure><img src="image_25.png" /><figcaption>Figure 20: Price-T (blue symbol '*' ) and MUSCL-Hancock (green symbol 'o') solution computed using a mesh with M = 200 cells, reference solution (red line) with M = 2000 cells. Derivative of the area A$_{a}$ and flux q$_{a}$ for test 1 of table 2. Γ = 1 Pa s m.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>(a) $a = a_L$</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(b) $a = a_R$</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(c) $a = h_0$</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(d) $a = r_0$</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(e) $a = E$</figcaption></figure>

Figure 21: MUSCL-Hancock solution computed using a mesh with $M = 200$ cells (green symbol 'o') and with $M = 2000$ cells (red line). Derivative of the area $A_a$ and flux $q_a$ for test 2 of table 2. $\Gamma = 1$ Pa s m.
<figure><img src="image_31.png" /><figcaption>Figure 22: MUSCL-Hancock solution computed using a mesh with $M = 200$ cells (green symbol 'o') and with $M = 2000$ cells (red line). Derivative of the area $A_a$ and flux $q_a$ for test 3 of table 2. $\Gamma = 1$ Pa s m.</figcaption></figure>
Finally, we have obtained the following first-order estimates for $\mu_X$ and $\sigma_X^2$,

$$ \mu_X = X(\mu_a), \quad \sigma_X^2 = \sum_{i=1}^{M} \sigma_{a_i}^2 X_{a_i}^2 + \sum_{\substack{i,j=1 \\ i \neq j}}^{M} X_{a_i} X_{a_j} \operatorname{cov}(a_i, a_j). \qquad (32) $$

Note that these formulas require to evaluate only once the state and sensitivity variables at the mean values of the uncertain parameters, which is really inexpensive compared to Monte-Carlo methods. We refer for instance to [10, 17] for more details.

Let us now present the numerical results. As uncertain parameters, we consider the ones that characterize the mechanical and geometrical properties of the vessels, i.e. the vessel thickness $h_0$, the radius at equilibrium $r_0$ and the Young modulus $E$, but also the the Riemann initial values for the area $A_L$ and $A_R$. We suppose all the parameters to be uncorrelated and, following the works of Petrella et al. [29] and Alastruey et al. [1], for each vessel we assume $h_0, r_0$ and $E$ to be respectively affected by an error of 2.5%, 3.5% and 5%. Finally, for $A_L$ and $A_R$ we take an error of 0.1%. Namely, we considered $\sigma_a = \operatorname{diag}((0.1\%A_L)^2, (0.1\%A_R)^2, (2.5\%h_0)^2, (3.5\%r_0)^2, (5\%E)^2)$.

In the following pictures, we show graphics for the average and the average plus/minus twice the standard deviation of the area $A$ and the flow $q$, as well as four samples. For the elastic model, outcomes with and without the correction term $\hat{S}$ are exhibited. Price-T and MUSCL-Hancock schemes are considered.

**Elastic Case.** Starting with test 1, in figures 23 and 24 we note that our confidence intervals are very satisfying in smooth regions. As expected, we observe large spikes near discontinuities when we do not take $\hat{S}$ into account. They disappear if we consider the correction term, however the confidence intervals are very narrow and the samples may fall outside. Note also the dependency of the shock speed with respect to the uncertain parameters has not been taken into account, which may cause errors, see for instance [10] for more details. Regarding the differences between the two numerical methods, we notice that the spikes of the Price-T are much smaller, which is in agreement with the sensitivities showed previously. Moreover, in correspondence of the head of the rarefaction, the Price-T solution exhibits little spikes which are emphasized by MUSCL-Hancock method. This could be related to the clipping of extreme values of the Price-T sensitivities. Observe also that the confidence intervals obtained with the Price-T scheme are more diffusive than MUSCL-Hancock, which is expected due to the numerical diffusion. However, let us also keep in mind that our approximations of the mean and variance for computing the confidence intervals are first-order, which may question the use of MUSCL-Hancock for the sensitivity equations.

Next in figure 25 we find the results of test 2 with 2-rarefaction waves. We see that the samples fall inside the predicted confidence intervals. On the other hand, we observe small spikes in correspondence of the heads of the rarefactions that reflect the behavior of the sensitivity variables.

Finally figure 26 show the outputs of Test 3. If not considering $\hat{S}$ entails the presence of large spikes in presence of shock, adding the correction term completely remove them. However, as we observed for test 1, the samples could fall outside of the confidence intervals when there are discontinuous solutions.

Finally, we can conclude this first batch of results saying that in general the sensitivity variables provide good confidence intervals, especially when the solutions are smooth. Moreover, recall that the computational cost is very low with respect to probabilistic method as Monte-Carlo. Again, while this scheme can require thousands of solutions of the state equations, our method needs only one solution of the state and as many solutions of the sensitivity equations as the number of uncertain parameters.

**Viscoelastic case.** In figures 27 and 28 we insert the results for test 1, for which we used respectively the Price-T and MUSCL-Hancock scheme. Let us recall that in this case there is no need to consider the correction term $\hat{S}$, as the viscoelastic state variables do not present shocks. Contrarily to the elastic system, now we do not have excessively large spikes. Solutions are smooth and the samples fall inside the predicted confidence intervals. Comparing the results obtained with the first and second-order scheme, we can do analogous observations to the ones for the elastic system. The confidence interval produced by the MUSCL-Hancock can be considered better as the ones of the Price-T scheme are more diffusive.

Similar considerations apply to tests 2 and 3, in graphics 29 and 30. In general for the viscoelastic system we do not have the problem of samples falling outside the confidence intervals, with the exception of the 2-shock case. Indeed, only one of the samples is not contained by the confidence interval for the cross-sectional area $A$. However, this happens only in a minor zone of the solution and with an error of order $O(10^{-8})$.
<figure><img src="image_32.png" /></figure>

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

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

<figure><img src="image_35.png" /><figcaption>Figure 23: Test 1 with Price-T scheme. Average and Average plus and minus twice the standard deviation in black. Four samples in magenta dashed line. Graphics for $A$ in (a) and (b), for $q$ in (c) and (d). $M = 500$ and $\Gamma = 0$.</figcaption></figure>
<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /><figcaption>Figure 24: Test 1 with MUSCL-Hancock scheme. Average and Average plus and minus twice the standard deviation in black. Four samples in magenta dashed line. Graphics for $A$ in (a) and (b), for $q$ in (c) and (d). $M = 500$ and $\Gamma = 0$.</figcaption></figure>

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

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

<figure><img src="image_40.png" /><figcaption>Figure 25: Test 2 MUSCL-Hancock schemes. Average and Average plus and minus twice the standard deviation in black. Four samples in magenta dashed line. Graphics for $A$ in (a) and (c), for $q$ in (b) and (d). $M = 500$ and $\Gamma = 0$.</figcaption></figure>

<figure><img src="image_41.png" /></figure>
<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /><figcaption>Figure 26: Test 3 with MUSCL-Hancock scheme. Average and Average plus and minus twice the standard deviation in black. Four samples in magenta dashed line. Graphics for $A$ in (a) and (b), for $q$ in (c) and (d). $M = 500$ and $\Gamma = 0$.</figcaption></figure>

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

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

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

<figure><img src="image_47.png" /><figcaption>Figure 27: Test 1 with Price-T scheme. Average and average plus and minus twice the standard deviation in black line for cross-sectional area $A$ (left) and flow $q$ (right). Four samples in magenta dashed line. $M = 500$ cells and $\Gamma = 1$ Pa m s.</figcaption></figure>
<figure><img src="image_48.png" /></figure>

Figure 28: Test 1 with MUSCL-Hancock scheme. Average and average plus and minus twice the standard deviation in black line for cross-sectional area A (left) and flow $q$ (right). Four samples in magenta dashed line. $M = 500$ cells and $\Gamma = 1$ Pa m s.

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

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

Figure 29: Test 2 with MUSCL-Hancock scheme. Average and average plus and minus twice the standard deviation (black line) for cross-sectional area A (left) and flow $q$ (right). Four samples in magenta dashed line. $M = 500$ cells and $\Gamma = 1$ Pa m s.

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

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

Figure 30: Test 3 with MUSCL-Hancock scheme. Average and average plus and minus twice the standard deviation (black line) for cross-sectional area A (left) and flow $q$ (right). Four samples in magenta dashed line. $M = 500$ cells and $\Gamma = 1$ Pa m s.

<figure><img src="image_53.png" /></figure>
<figure><img src="image_54.png" /><figcaption>Figure 31: Idealised model for vessel junction (N = 3).</figcaption></figure>

# 5 Numerical simulation of a human arterial network

The main objective of this section is to model and simulate a complex network made of 37 arteries. As a first step, it is necessary to propose a numerical treatment of the junctions between two or several arteries, considering both the state and the sensitivity equations.

## 5.1 Treatment of the junctions

As just motivated, we want to explain how to treat the presence of junctions points among vessels. For the sake of simplicity, we will illustrate in more details the cases of two and three vessels but the strategy can be extended to an arbitrary number, for more details we refer for instance to [24] and [26].

In the following, we will assume that our network of vessels is represented by a finite collection of incoming and outgoing directed arcs in one space dimension, which are connected by nodes and aligned in the same direction.

Considering N vessels sharing a node as represented in figure 31 for N = 3,

$$ \left\{ \begin{array}{l} \partial_t \mathbf{Q}_k + \partial_x \mathbf{F}(\mathbf{Q}_k) + \mathbf{A}(\mathbf{Q}_k) \partial_x \mathbf{Q}_k = \mathbf{S}(\mathbf{Q}_k) \\ \mathbf{Q}_k(x, 0) = \mathbf{Q}_k^{1D} \end{array} \right. \qquad (33) $$

with $k = 1, \dots, N$, where $\mathbf{Q}_k^{1D}$ represents the constant initial value in the k-th vessel. Waves with negative speed in incoming vessels and positive speed in outgoing vessels are expected to develop at the node, leading to the occurrence of N intermediate states denoted by $\mathbf{Q}_{*k}$. These states represent the traces of the solution at the node in each vessel, and should be connected to the initial states $\mathbf{Q}_k^{1D}$ by (approximate) nonlinear waves. In the following and for the sake of simplicity, we will always consider that

$$ \mathbf{Q}_k^{1D} = \begin{pmatrix} A_k^{1D} \\ q_k^{1D} \\ \Psi_k^{1D} \end{pmatrix} \quad \text{and} \quad \mathbf{Q}_{*k} = \begin{pmatrix} A_{*k} \\ q_{*k} \\ \Psi_{*k} \end{pmatrix}, $$

being implied that, unless otherwise stated, the last component $\Psi$ should not be considered when dealing with the elastic system.

In order to define the N unknown states $\mathbf{Q}_{*k}$, we suggest to impose

* the continuity of the mass flux $q = Au$ through the junction,

* the continuity of the total pressure $p + \frac{1}{2}\rho u^2$ through the junction,

* the constancy of the generalized Riemann invariants seen in section 1 in each vessel.
Note that the continuity of the total pressure is nothing but the Bernoulli relation which is easily recovered from the state equations assuming that the flow is incompressible (A is constant) and stationary at the junction. As far as the constancy of the generalized Riemann invariants is concerned, it means implicitly that the nonlinear waves created at the junctions and propagating in the vessels are rarefaction waves. Then, for $N = 2$ we obtain the following system

$$ 
\left\{
\begin{aligned}
& A_1^* u_1^* = A_2^* u_2^*, \\
& p(A_1^*, \Psi_1^*) + \frac{1}{2} \rho u_1^{*2} = p(A_2^*, \Psi_2^*) + \frac{1}{2} \rho u_2^{*2} \\
& \frac{A_k^*}{\epsilon} + \Psi_k^* = \frac{A_k^{1D}}{\epsilon} + \Psi_k^{1D}, && k = 1, 2, \\
& u_1^* + \int A_1^* \frac{\tilde{c}(A)}{A} dA = u_1^{1D} + \int A_1^{1D} \frac{\tilde{c}(A)}{A} dA, && \text{(left rarefaction)}, \\
& u_2^* - \int A_2^* \frac{\tilde{c}(A)}{A} dA = u_2^{1D} - \int A_2^{1D} \frac{\tilde{c}(A)}{A} dA, && \text{(right rarefaction)},
\end{aligned}
\right.
\tag{34}
$$

whereas, if we have a single incoming vessel that divides into two different outgoing ones ($N = 3$), we will have three more equations, namely

$$ 
\left\{
\begin{aligned}
& A_1^* u_1^* = \sum_{k=2}^{3} A_k^* u_k^*, \\
& p(A_1^*, \Psi_1^*) + \frac{1}{2} \rho u_1^{*2} = p(A_k^*, \Psi_k^*) + \frac{1}{2} \rho u_k^{*2}, && k = 2, 3, \\
& \frac{A_k^*}{\epsilon} + \Psi_k^* = \frac{A_k^{1D}}{\epsilon} + \Psi_k^{1D}, && k = 1, 2, 3, \\
& u_1^* + \int A_1^* \frac{\tilde{c}(A)}{A} dA = u_1^{1D} + \int A_1^{1D} \frac{\tilde{c}(A)}{A} dA, && \text{(left rarefaction)}, \\
& u_k^* - \int A_k^* \frac{\tilde{c}(A)}{A} dA = u_k^{1D} - \int A_k^{1D} \frac{\tilde{c}(A)}{A} dA, && k = 2, 3 \text{ (right rarefaction)}.
\end{aligned}
\right.
\tag{35}
$$

Note that in the first equation of system (35) the flux in the incoming vessel has to be equal to the sum of the fluxes in the two outgoing vessels. Instead the pressure is assumed to be the same in all three vessels, as we can see from the second equation.

**Linearization and sensitivity variables.** It is observed that the proposed coupling conditions are highly nonlinear due to the integral terms, which poses difficulties to compute explicitly the intermediate states, but also to define the coupling conditions for the sensitivity variables. With this in mind, we suggest to linearize (34) and (35). In the spirit of pressure relaxation methods, see for instance the pioneering work [22, 30], or more recently [9], we suggest to replace the real wave speed $\tilde{c}(A)$ by

$$ \tilde{c}(A) = \frac{\bar{a}}{A}, $$

where $\bar{a}$ is a constant defined by the so-called sub-characteristic condition

$$ \bar{a} = \max_k (A_k^{1D} \tilde{c}(Q_k^{1D})). $$

Therefore, relations (34) become

$$ 
\left\{
\begin{aligned}
& A_1^* u_1^* = A_2^* u_2^*, \\
& p(A_1^*, \Psi_1^*) + \frac{1}{2} \rho u_1^{*2} = p(A_2^*, \Psi_2^*) + \frac{1}{2} \rho u_2^{*2} \\
& \frac{A_k^*}{\epsilon} + \Psi_k^* = \frac{A_k^{1D}}{\epsilon} + \Psi_k^{1D}, && k = 1, 2, \\
& u_1^* - \frac{\bar{a}_k}{A_1^*} = u_1^{1D} - \frac{\bar{a}_k}{A_1^{1D}}, && \text{(left rarefaction)}, \\
& u_2^* + \frac{\bar{a}_k}{A_2^*} = u_2^{1D} + \frac{\bar{a}_k}{A_2^{1D}}, && \text{(right rarefaction)},
\end{aligned}
\right.
\tag{36}
$$

and analogously for system (35). Therefore and compared to [26], our procedure to couple two or more vessels is simpler. We highlight that it is important to reduce the computational cost for the junctions as they would be the most expensive part when computing the outputs of a network. In addition and as we will see in the numerical experiments, the simplification gives good results.
At last, it also allows to easily obtain the coupling conditions also for the sensitivity variables by differentiating the linearized conditions (36) with respect to the parameter a. Therefore for $N = 2$ we write

$$ egin{cases} 
A_1^* u_{a1}^* + A_{a1}^* u_1^* = A_2^* u_{a2}^* + A_{a2}^* u_2^*, \\ 
p_a(A_1^*, \Psi_1^*, A_{a1}^*, \Psi_{a1}^*) + \rho u_1^* u_{a1}^* = p_a(A_2^*, \Psi_2^*, A_{a2}^*, \Psi_{a2}^*) + \rho u_2^* u_{a2}^* \\ 
\frac{A_{ak}^*}{\epsilon} - \frac{A_k^*}{\epsilon^2} \epsilon_a + \Psi_{ak}^* = \frac{A_{ak}^{1D}}{\epsilon} - \frac{A_k^{1D}}{\epsilon^2} \epsilon_a + \Psi_{ak}^{1D} & k = 1, 2, \quad (37) \\ 
u_{a1}^* + \frac{\bar{a}}{A_1^{*2}} A_{1a}^* = u_{a1}^{1D} + \frac{\bar{a}}{A_1^{1D2}} A_{1a}^{1D}, & \text{(left rarefaction),} \\ 
u_{a2}^* - \frac{\bar{a}}{A_2^{*2}} A_{2a}^* = u_{a2}^{1D} - \frac{\bar{a}}{A_2^{1D2}} A_{2a}^{1D}, & \text{(right rarefaction)} 
\end{cases} $$

while, for $N = 3$:

$$ egin{cases} 
A_1^* u_{a1}^* + A_{a1}^* u_1^* = A_2^* u_{a2}^* + A_{a2}^* u_2^* + A_3^* u_{a3}^* + A_{a3}^* u_3^*, \\ 
p_a(A_1^*, \Psi_1^*, A_{a1}^*, \Psi_{a1}^*) + \rho u_1^* u_{a1}^* = p_a(A_k^*, \Psi_k^*, A_{ak}^*, \Psi_{ak}^*) + \rho u_k^* u_{ak}^* & k = 2, 3 \\ 
\frac{A_{ak}^*}{\epsilon} - \frac{A_k^*}{\epsilon^2} \epsilon_a + \Psi_{ak}^* = \frac{A_{ak}^{1D}}{\epsilon} - \frac{A_k^{1D}}{\epsilon^2} \epsilon_a + \Psi_{ak}^{1D}, & k = 1, 2, 3, \quad (38) \\ 
u_{a1}^* + \frac{\bar{a}}{A_1^{*2}} A_{1a}^* = u_{a1}^{1D} + \frac{\bar{a}}{A_1^{1D2}} A_{1a}^{1D}, & \text{(left rarefaction),} \\ 
u_{ak}^* - \frac{\bar{a}}{A_k^{*2}} A_{ka}^* = u_{ak}^{1D} - \frac{\bar{a}}{A_k^{1D2}} A_{ka}^{1D}, & k = 2, 3, \quad \text{(right rarefaction).} 
\end{cases} $$

Note that we have done a further simplification by assuming that $\bar{a}$ is independent from the uncertain parameter a.

**Remark.** For the numerical computations, the above-mentioned linearization procedure to simplify the coupling conditions was used for the elastic state and sensitivity variables as well. However, we highlight that, in this case, relations (36) are not reduced to the classical ones we would have if we had used the elastic Riemann invariants $I_e^- = u + 4c$ and $I_e^+ = u - 4c$. See [29] for further details.

**Numerical illustrations.** Here we illustrate the numerical outcomes obtained with our implementation of the junctions. We want to compare the results given by a numerical simulation considering two identical vessels linked by a junction and the one considering a unique equivalent vessel. This academic test case is useful to judge the efficiency of the strategy proposed to couple the vessels is efficient or not. We will also consider a test case with three vessels.

Let us briefly describe the boundary conditions. Regarding the left boundary condition, we impose an inlet flow given by

$$ q_{bc} = \hat{q} e^{-10000(t-0.025)^2} $$

where $\hat{q} = 100 \text{ ml s}^{-1}$, and we define the missing area and relaxation variable by exploiting the continuity of Riemann invariants, which leads to the following system,

$$ egin{cases} 
F_1(A^*, \Psi^*) = u_1 + \frac{\bar{a}}{A_1} - u^* - \frac{\bar{a}}{A^*} = 0 \\ 
F_2(A^*, \Psi^*) = \frac{A_1}{\epsilon} + \Psi_1 - \frac{A^*}{\epsilon} - \Psi^* = 0. 
\end{cases} 
\quad (39) $$

Here the subscript 1 refers to the value of the first cell. This system can be solved with the Newton-Rhapson Method. Regarding the sensitivity variables, we differentiate (39) to obtain:

$$ egin{cases} 
F_{a1}(A^*, \Psi^*, A_a^*, \Psi_a^*) = u_{a1} - \frac{\bar{a}}{A_1^2} A_{a1} - u_{a1}^* + \frac{\bar{a}}{A^*} A_a^* = 0 \\ 
F_{a2}(A^*, \Psi^*, A_a^*, \Psi_a^*) = \frac{A_{a1}}{\epsilon} - \frac{A_1}{\epsilon^2} \epsilon_a + \Psi_{a1} - \frac{A_{a1}^*}{\epsilon} + \frac{A^*}{\epsilon^2} \epsilon_a - \Psi_a^* = 0 
\end{cases} $$

The right boundary conditions are classical transmissive conditions.

The parameters of the vessels under consideration are given in the table below. More precisely, for the parameters of the two connected vessels we used the ones of vessels $v_2$ and $v_3$, and for the equivalent vessel we do refer to $v_1$. The missing parameters are taken from table 1. Here we are also adding the friction forces with $\mu = 2.5 \text{mPa S}$.

As initial condition for the state variables we used a Gaussian, namely

$$ Q = \begin{pmatrix} A \\ q \\ \Psi \end{pmatrix} = \begin{pmatrix} A_0(1 + \alpha e^{-\beta(x-gate)}) \\ 0 \\ 0 \end{pmatrix}. \quad (40) $$
<table><caption>Table 3: Parameters for different vessels in subsection (5.1).</caption><thead><tr><th>Vessel name</th><th>L [m]</th><th>r<sub>0</sub> [m]</th><th>h<sub>0</sub> [m]</th><th>Γ [Pa s m]</th></tr></thead><tbody><tr><td>v<sub>1</sub></td><td>0.04</td><td>0.01</td><td>5.0 × 10<sup>-4</sup></td><td>1.0</td></tr><tr><td>v<sub>2</sub></td><td>0.02</td><td>0.01</td><td>5.0 × 10<sup>-4</sup></td><td>1.0</td></tr><tr><td>v<sub>3</sub></td><td>0.02</td><td>0.01</td><td>5.0 × 10<sup>-4</sup></td><td>1.0</td></tr><tr><td>v<sub>4</sub></td><td>0.02</td><td>0.001</td><td>5.0 × 10<sup>-5</sup></td><td>3.0</td></tr></tbody></table>

with $\alpha = 0, \beta = 10000$. Therefore for the sensitivity variables we imposed

$$ \mathbf{Q}_a = \begin{pmatrix} A_a \\ q_a \\ \Psi_a \end{pmatrix} = \begin{pmatrix} A_{0a}(1 + \alpha e^{-\beta(x-gate)}) \\ 0 \\ 0 \end{pmatrix}. \qquad (41) $$

We show the results obtained with the MUSCL-Hancock scheme in the viscoelastic case. We inserted the results for the flow and pressure in figure 32. In 33, 34 and 35 we show the outputs for the sensitivity variables obtained with respectively $a = h, a = r_0$ and $a = E$. Note that when we differentiate the sensitivity equations with respect to the radius $r_0$, we are also considering uncertainty on the initial condition as it depends on $r_0$. We plotted the results for three different times with MUSCL-Hancock scheme. Even if linearized coupling conditions are used, one can be satisfied with the results in spite of using a coarse mesh. Regarding the sensitivities, the outputs seem to be slightly worse than the ones for flow and pressure, especially at time $t = 0.068$. Since the sensitivity equations are derived from the state ones, an imprecision in the latter will be reflected from the sensitivity. Thus, the error will be amplified as we are adding two of them: one from the state and one from our modeling of the junctions.

To be thorough, in 36, 37, 38 and 39 we also included the graphics obtained considering one mother vessel ($v_2$) which splits up into two daughter vessels ($v_3$ and $v_4$). Note that $v_3$ has the same characteristic of $v_2$, while the viscoelastic parameter $\Gamma$ of $v_4$ is three times $v_3$'s $\Gamma$. Also the radius $r_0$ and the vessel thickness $h$ are much smaller. This will entail that the blood flow in $v_4$ will be minor than the one of $v_3$ as we can see in (36). We remark that when considering three vessels, at the junction we are imposing the flow (or its derivative) of the mother vessel to be equal to the sum of the flows (or the sum of the derivatives of the fluxes) of the two daughter vessels. We can observe that these two facts seem to be fulfilled both from the state and the sensitivity variables.

## 5.2 A human arterial network

Following the works of Matthys et al. [24] and Alastruey et al. [1], in this subsection we aim to model an arterial network for both the state and the sensitivity equations. They constructed a 1 : 1 in-vitro model using silicone tubes for the 37 largest arteries of the systemic circulation, a pump to simulate the heart flow and terminal resistance tubes to represent the peripheral circulation. In particular they created a closed loop hydraulic system connecting these tubes to an overflow reservoir that represents a constant venous pressure. For the values of the parameters we refer the reader to [24]. We underline that they were measured directly from the in-vitro model, thus no data fitting is involved. Note that compared to [1], we will use the same parameters values with the exception of the radius which we kept constant, taking it as the average between the initial and final values for each vessel. Regarding the spatial discretization, in general we used $\Delta x = 2$ cm with the exception of the vessels shorter than 1.5 cm, for which we used only one cell. In this way we were legitimized to use $\varepsilon = 10^{-3}$. Moreover we took CFL = 0.9 and consequently we had a time step around $\Delta t = 700\mu s$ in the elastic case, 7 times larger than the one used in [1]. Instead in the viscoelastic case we had $\Delta t = 400\mu s$, 20 times larger. For the initial conditions, once again we used (40) and (41). Finally, since we want to obtain periodic solutions, we show graphics from the 14th cardiac cycle, with period $T = 0.827s$ (refer to [29] for more details).

**Boundary conditions.** In order to simulate the heart, the heart flow showed in figure (40) is imposed at the beginning of the ascending aorta. Then to find also the boundary values for the area and the relaxation variable, we exploit the right Riemann invariants as before. Regarding the right boundary conditions for terminal vessels, we follow [1] and we consider the flow

$$ q^* = \frac{p(A^*, \Psi^*) - p_{out}}{R_p} \qquad (42) $$
<figure><img src="image_55.png" /><figcaption>Figure 32: Flow and pressure obtained with MUSCL-Hancock method. Red symbol 'o' for the unique vessel, green and magenta symbol 'x' for the two coupled vessels. Mesh with $M = 120$ cells for the unique vessel.</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>Figure 33: Derivative of flow and pressure with respect to $a = h_0$ obtained MUSCL-Hancock method. Red symbol 'o' for the unique vessel, green and magenta symbol 'x' for the two coupled vessels. Mesh with $M = 120$ cells for the unique vessel.</figcaption></figure>
<figure><img src="image_57.png" /><figcaption>Figure 34: Derivative of flow and pressure with respect to $a = r_0$ obtained with MUSCL-Hancock method. Red symbol 'o' for the unique vessel, green and magenta symbol 'x' for two coupled vessels. Mesh with $M = 120$ cells for the unique vessel.</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>Figure 35: Derivative of flow and pressure with respect to $a = E$ obtained with MUSCL-Hancock method. Red symbol 'o' for the unique vessel, green and magenta symbol 'x' for two coupled vessels. Mesh with $M = 120$ cells for the unique vessel.</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>Figure 36: Flow and pressure obtained with MUSCL-Hancock method. Green, blue and magenta symbols represent $v_2$, $v_3$ and $v_4$ respectively. Mesh with $M = 60$ cells for each vessel.</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>Figure 37: Derivative of flow and pressure with respect to $a = h_0$ obtained MUSCL-Hancock method. Green, blue and magenta symbols represent $v_2$, $v_3$ and $v_4$ respectively. Mesh with $M = 60$ cells for each vessel.</figcaption></figure>
<figure><img src="image_61.png" /><figcaption>Figure 38: Derivative of flow and pressure with respect to $a = r_0$ obtained with MUSCL-Hancock method. Green, blue and magenta symbols represent $v_2$, $v_3$ and $v_4$ respectively. Mesh with $M = 60$ cells for each vessel.</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>Figure 39: Derivative of flow and pressure with respect to $a = E$ obtained with MUSCL-Hancock method. Green, blue and magenta symbols represent $v_2$, $v_3$ and $v_4$ respectively. Mesh with $M = 60$ cells for each vessel.</figcaption></figure>
<figure><img src="image_63.png" /><figcaption>Figure 40: Heart flow imposed at the beginning of the ascending aorta.</figcaption></figure>

with $p_{out} = 3.2$ mmHg induced by the overflow reservoir. Using left Riemann invariants and asking for their constancy, we find

$$ \left\{ \begin{aligned} F_1(A^*, \Psi^*) &= \mathrm{RI}_M A^* R_p + \bar{a} R_p - p(A^*, \Psi^*) + p_{out} = 0 \\ F_2(A^*, \Psi^*) &= \frac{A_M}{\bar{\epsilon}} + \Psi_M - \frac{A^*}{\bar{\epsilon}} - \Psi^* = 0, \end{aligned} \right. \qquad (43) $$

with $RI_M = u_M - \frac{\bar{a}}{A_M}$, where the subscript M indicates the last cell of the vector. Deriving (42) and (43), for the sensitivity variables we obtain

$$ q_a^* = \frac{p_a(A^*, \Psi^*, A_a^*, \Psi_a^*)}{R_p} = \frac{p_a^*}{R_p} $$

and

$$ \left\{ \begin{aligned} F_{a1}(A^*, \Psi^*, A_a^*, \Psi_a^*) &= \mathrm{RI}_{aM} A^* R_p + \mathrm{RI}_M A_a^* R_p - p_a^* = 0 \\ F_{a2}(A^*, \Psi^*, A_a^*, \Psi_a^*) &= \frac{A_{aM}}{\bar{\epsilon}} - \frac{A_M}{\bar{\epsilon}^2} \epsilon_a + \Psi_{aM} - \frac{A_a^*}{\bar{\epsilon}} + \frac{A_a^*}{\bar{\epsilon}^2} \epsilon_a - \Psi_a^* \end{aligned} \right. $$

with $RI_{aM} = u_{aM} + \frac{\bar{a}A_{aM}}{A_M^2}$.

*Numerical illustrations.* First of all we show the numerical results for the pressure and flow computed in four locations, namely in the midpoint of the thoracic aorta I, left renal, right iliac femoral III and right carotid. We compared the outputs obtained with the first and second order scheme, in both elastic (figure 41) and viscoelastic (figure 42) case. Then, using Price-T and MUSCL-Hancock method respectively, in graphics 43 and 44 we observe the dissimilarities between the conservative and non-conservative system. We use the ENO slope for the second-order method.

Both elastic and viscoelastic outcomes appear to have a similar profile with respect to the ones given by the experimental measurements. Indeed, this is an important achievement as it validates our numerical results. Then, we highlight that the use of the simple resistive boundary conditions produces non-physiological oscillations, as observed for instance in [1, 29]. As a consequence, high order methods overestimate the actual oscillations present in the experimental solution, leading to the contradictory situation in which the first-order schemes seem to approximate better the experimental measurements. Indeed, the numerical viscosity entails a damping of the oscillations when using a first-order method as the Price-T. This kind of problem is eluded by exploiting the viscoelastic formulation, which reduces the oscillations. Note that both the left renal and right carotid are terminal vessels and exhibit large oscillations if compared to the experimental ones.

We are now interested in plotting the average and variance of the flow and pressure, and then their confidence intervals. Regarding the uncertain parameters, we take the ones linked to the arterial stiffness, hence the radius at equilibrium $r_0$, the vessel thickness $h_0$ and the Young modulus $E$. We assume them to be affected by the same errors considered in subsection 4.3. Note that assuming $a = r_0$ means modifying the initial conditions for the SE as well, reminding that for each vessel we use conditions (40) and (41).

First of all in graphic 45 we insert the outputs for the sensitivity variables, that is to say the derivative of the flow and pressure with respect to $a = h_0$, $a = r_0$ and $a = E$. We computed them in the midpoint of the left renal with the MUSCL-Hancock scheme for both the elastic and viscoelastic system. As we said before, neglecting viscoelasticity produces overestimation of the oscillations of the experimental solution and this consequently affects the behavior of the elastic
<figure><img src="image_64.png" /></figure>

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

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

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

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

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

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

<figure><img src="image_71.png" /><figcaption>Figure 41: Comparison among experimental (exp, red line), Price-T (black line) and MUSCL-Hancock (MH, blue line) results for pressure (left) and flow (right) computed in the midpoint of the thoracic aorta I, left renal, r. iliac femoral III and right carotid. Elastic case ($Γ = 0$).</figcaption></figure>
<figure><img src="image_72.png" /></figure>

Thoracic Aorta I

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

Thoracic Aorta I

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

L. Renal

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

L. Renal

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

R. Iliac Femoral III

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

R. Iliac Femoral III

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

R. carotid

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

R. carotid

<caption>Figure 42: Comparison among experimental (exp, red line), Price-T (black line) and MUSCL-Hancock (MH, blue line) results for pressure (left) and flow (right) computed in the midpoint of the thoracic aorta I, left renal, r. iliac femoral III and right carotid. Viscoelastic case ($\Gamma \neq 0$).</caption>
<figure><img src="image_80.png" /></figure>

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

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

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

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

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

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

<figure><img src="image_87.png" /><figcaption>Figure 43: Comparison among experimental (exp, red line), numerical elastic (elas, black line) and viscoelastic (visc, dashed black line) pressure (left) and flow (right) computed in the midpoint of the thoracic aorta I, left renal, r. iliac femoral III and right carotid. Results obtained with Price-T method.</figcaption></figure>
<figure><img src="image_88.png" /></figure>

Thoracic Aorta I

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

L. Renal

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

R. Iliac Femoral III

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

R. carotid

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

Thoracic Aorta I

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

L. Renal

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

R. Iliac Femoral III

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

R. carotid

<caption>Figure 44: Comparison among experimental (exp, red line), numerical elastic (elas, black line) and viscoelastic (visc, dashed black line) pressure (left) and flow (right) computed in the midpoint of the thoracic aorta I, left renal, r. iliac femoral III and right carotid. Results obtained with MUSCL-Hancock method.</caption>
sensitivity variables. Indeed, especially for the derivative of the flow, we note large oscillations that are damped when using the viscoelastic formulation. Clearly, this leads to the same problem for the average and in particular for the variance of the flow, both of them shown in picture 46. In the same figure we insert the results for the average and variance of the pressure as well. We note that for them the dissimilarities between the two systems are less outstanding. However, we observe that the viscoelastic pressure variance results to be slightly bigger than the elastic one, a fact that could be due to the difference in values between the derivatives with respect to $a = r_0$ of the elastic and viscoelastic pressure, see figure 45(c).

Next, we insert the *maximum* in absolute value of the standard deviation for the pressure and flow for all the 37 vessels in graphic 47, in the elastic and viscoelastic case respectively. In these graphics we used different colors to discriminate the vessels with and without the resistive boundary conditions, in magenta and blue respectively. As we said, we note that in general the standard deviation of the viscoelastic pressure is slightly greater than the one of the elastic pressure. Then, for the elastic model, we observe that the maximum values of the standard deviation of the flow are much bigger than the viscoelastic flow ones. Once again, this could be due to the large oscillations we noted when using the elastic formulation. However, we highlight that in both cases the maximum values for the standard deviation of the flow are greater for non-terminal vessels. This could be related to the fact that these vessels compose the aortic path and therefore are larger and less resistive. Instead terminal vessels are usually smaller, and thus their flow is less pulsatile, where pulse wave speed is proportional to the stiffness, and as such to the uncertain parameters we have considered. This could explain why their flow is less sensitive to changes in the values of the parameters. Differences in values for the standard deviation of pressure between non-terminal and peripheral vessels are less outstanding.

Then, in order to understand which parameter the pressure and flow are most sensitive to, we show the maximum values of standard deviation considering only one value for $a$ each time, see graphic 48. For both variables we observe larger values when exploiting $a = r_0$, in particular when considering non-terminal vessel. We obtained similar results in both viscoelastic and elastic case, therefore we do report the outputs only in the first case.

Finally, we plot the average for the flow and pressure plus and minus twice the standard deviation. Indeed, figures 49 and 50 show the outputs for the elastic and viscoelastic case respectively. MUSCL-Hancock method was used. In general the pressure seems to be more sensitive than the flow to changes in the values of the parameters. We observe that the predicted interval of the flow computed with MUSCL-Hancock in elastic case is larger than the one in the viscoelastic case. As we already said, this could be due to the large oscillations we noted for both the state and sensitivity elastic variables.

# 6 Conclusion and perspectives

In this work we studied the blood flow equations with both elastic [37] and viscoelastic [25] formulations. In the latter case, an additional term is needed in order to express the viscoelastic behaviour of the arteries walls. This term includes a second-order spatial derivative of the flow, yielding to a parabolic advection-diffusion-reaction system. In order to transform it into a hyperbolic advection-reaction system as the elastic one, we followed the work of Montecinos et al. [25]. Thus, we used a relaxation approach which entails the introduction of a new variable $\Psi$ and relaxation time $\varepsilon$. This step was necessary in order to write these systems with the same formalism and to find a second-order method which could be applied to both of them.

Subsequently we introduced the centred first-order Price-T scheme, which is a variant of FORCE method for non-conservative systems [34]. In particular, we have reformulated the scheme so that it automatically reduces to the FORCE method when considering a conservative system; a crucial request as primitive schemes applied to a conservative system could give a solution which converges to the wrong one. Hence, we needed a second-order scheme written in non-conservative form with this property. We turned to the MUSCL-Hancock Price-T scheme presented in [34], which however does not reduce to a conservative method. Therefore we modified it, adding a correction term, in order to achieve the SLIC scheme [31] when considering a conservative system as the elastic one. Indeed, for this case we showed that the numerical output of the original MUSCL-Hancock Price-T scheme did not converge to the exact one, contrary to the solution of our version of the method. Considering the ADER approach [32, 34] could allow the construction of higher order numerical schemes with this property.

The second goal of this work was to apply the continuous sensitivity equations method [10] to the blood flow equations, with the aim to exploit it for uncertainty quantification. One of the reasons for choosing this method could be its low computational cost, especially if we compare it to probabilistic methods as Monte Carlo. Thus, differentiating the state equations with respect to an arbitrary uncertain parameter $a$, we found the so-called sensitivity system. However, an
<figure><img src="image_96.png" /><figcaption>(a) Derivative of pressure with $a = h_0$</figcaption></figure>

<figure><img src="image_97.png" /><figcaption>(b) Derivative of flow with $a = h_0$</figcaption></figure>

<figure><img src="image_98.png" /><figcaption>(c) Derivative of pressure with $a = r_0$</figcaption></figure>

<figure><img src="image_99.png" /><figcaption>(d) Derivative of flow with $a = r_0$</figcaption></figure>

<figure><img src="image_100.png" /><figcaption>(e) Derivative of pressure with $a = E$</figcaption></figure>

<figure><img src="image_101.png" /><figcaption>(f) Derivative of flow with $a = E$</figcaption></figure>

Figure 45: Derivative of pressure and flow in the midpoint of left renal in elastic (magenta line) and viscoelastic case (blue line). Results obtained with MUSCL-Hancock method.
<figure><img src="image_102.png" /><figcaption>(a) Average of flow</figcaption></figure>

<figure><img src="image_103.png" /><figcaption>(c) Average of pressure</figcaption></figure>

<figure><img src="image_104.png" /><figcaption>(b) Variance of flow</figcaption></figure>

<figure><img src="image_105.png" /><figcaption>(d) Variance of pressure</figcaption></figure>

Figure 46: Average and variance of flow and pressure in the midpoint of left renal in elastic (magenta line) and viscoelastic case (blue line). Results obtained with MUSCL-Hancock method.

<figure><img src="image_106.png" /><figcaption>(a) Elastic case</figcaption></figure>

<figure><img src="image_107.png" /><figcaption>(b) Viscoelastic case</figcaption></figure>

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

Figure 47: Maximum absolute value of standard deviation of pressure (top) and flow (bottom) in the midpoint of each vessels of the network (identified with numbers from 1 to 37). Magenta color for terminal vessels. Elastic (left) and viscoelastic (right) case, MUSCL-Hancock method.
<figure><img src="image_109.png" /><figcaption>Figure 48: Maximum absolute value of standard deviation of pressure (left) and flow (right) using $a = h_0$ (top), $a = r_0$ (middle) and $a = E$ (bottom) in the midpoint of each vessel of the network (identified with numbers from 1 to 37). Magenta color for terminal vessels. Viscoelastic case, MUSCL-Hancock method.</figcaption></figure>

assumption of the CSE method is the regularity of the solution, which could not hold in some cases (e.g. Riemann problems). The consequence is the appearance of Dirac delta functions, which can be compensated by introducing the source term $\hat{S}$ in the sensitivity system, refer to [8, 10, 17]. However special treatment is required to handle this term when applying the numerical schemes. Numerical results for the elastic system show that the addition of $\hat{S}$ reduces considerably or completely the large spikes that otherwise would be present in correspondence of the shocks of the state variables. Whereas, regarding the viscoelastic system, there is no need to consider the source term $\hat{S}$ for the sensitivities. Indeed, in this case the state solution does not present discontinuities due to the diffusive term present in the viscoelastic closure condition (4).

We specify that to define the term $\hat{S}$ we exploited an HLL approach as the elastic system has a two-waves structure and the viscoelastic equations do not require it. However, for systems with three or more equations this strategy would be incomplete and thus, the use of the HLLC approach could lead to an improvement of the results, refer for instance to [35, 10].

Finally we computed the confidence intervals for the cross-sectional area and the flow for three different Riemann problems. Regarding the uncertain parameters, we considered the initial values $A_L$ and $A_R$ for the cross-sectional area and then the parameters linked to the arterial stiffness, thus the radius at equilibrium $r_0$, the vessel thickness $h_0$ and the Young modulus $E$. In general, the numerical outputs were satisfying as the confidence intervals included the samples, with the exception of shock zones. This could be related to the fact that we neglected the dependency of the shock speed on the uncertain parameters. Moreover, in correspondence of the discontinuities of the state variables, we note large spikes that are removed when the source term $\hat{S}$ is active. On the other hand, in the majority of the results for the viscoelastic system the samples fall inside the confidence intervals, as in this case the variables do not present shocks and thus the predicted intervals are sufficiently wide. Ulterior enhancement on the computation of confidence intervals could be achieved using higher order sensitivities.

The last step of this work was to apply the proposed numerical models and schemes to an arterial network of 37 silicone vessels, for which we refer to [1]. This was an important request, as we would be able to validate our results. The first problem was to define the coupling conditions for the sensitivity unknowns. In order to be able to obtain them, we had to simplify the conditions for the state variables, namely doing a relevant linearization. Nevertheless the numerical outputs are satisfying, even if the results for the sensitivity variables are slightly worse as they sum the state and the modeling errors. It could be interesting to attempt to use more complex coupling conditions, with the remarks of paying attention to the computational cost, as the junctions are one of the most expensive part when considering a blood vessels network.
<figure><img src="image_110.png" /></figure>

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

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

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

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

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

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

<figure><img src="image_117.png" /><figcaption>Figure 49: Comparison among experimental results (exp, red line), average (black line) plus and minus twice the standard deviation (magenta line) of pressure (left) and flow (right) computed in the midpoint of the thoracic aorta I, left renal, r. iliac femoral III and right carotid. MUSCL-Hancock method, elastic case.</figcaption></figure>
<figure><img src="image_118.png" /></figure>

Thoracic Aorta I

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

L. Renal

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

R. Iliac Femoral III

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

R. carotid

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

Thoracic Aorta I

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

L. Renal

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

R. Iliac Femoral III

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

R. carotid

<caption>Figure 50: Comparison among experimental results (exp, red line), average (black line) plus and minus twice the standard deviation (magenta line) of pressure (left) and flow (right) computed in the midpoint of the thoracic aorta I, left renal, r. iliac femoral III and right carotid. MUSCL-Hancock method, viscoelastic case.</caption>
Afterwards, we exhibited our results for the state variables in the network case, and showed that the outputs had sim-
ilar profiles to the experimental ones. An interesting remark is that, when using the elastic formulation, we observe the
contradictory result that low-order methods approximate the experimental outcomes better than the high-order ones. This
is due to the fact that, to simulate the peripheral resistances, we used simple terminal resistive models which produce non-
physiological oscillations in the numerical results. While in the low-order solutions the numerical viscosity reduces these
oscillations, high order methods overestimate the ones present in the experimental results. This problem is overcome using
the viscoelastic formulation as it introduces a diffusive term which dumps the oscillations.

Finally, we computed the confidence intervals using the uncertain parameters linked to the arterial stiffness, as it is associated to many human illnesses. Then, we observed that the maximum standard deviation for terminal vessels' flow reaches lower values than that of the vessels which compose the aortic path. This could be due to the fact that their flow is less pulsatile, where pulse wave speed is linked to the arterial stiffness. An ulterior step could be the study of the sensitivity of the blood flow and pressure not only in relation to the arterial stiffness, but also with respect to the inlet and outlet flow, as they can depend on many parameters which are difficult to estimate as well.

References

[1] J. Alastruey, A. W. Khir, K. S. Matthys, P. Segers, S. J. Sherwin, P. R. Verdonck, K. H. Parker, J. Peiro. Pulse wave propa-
gation in a model human arterial network: Assessment of 1-D visco-elastic simulations against in vitro measurements.
*Journal of Biomechanics*, 44:2250-2258, 2011.

[2] A. Ambroso, C. Chalons, P. A. Raviar. A Godunov-type method for the seven-equation model of compressible two-phase
flow. 2010. hal-00517375

[3] P. J. Blanco, R. A. Feijoo, S. A. Urquiza. A unified variational approach for coupling 3D-1D models and its blood flow applications. *Comput. Methods Appl. Mech. Engrg.*, 196 (2007) 4391-4410.

[4] J. Borggaard and J. Burns. A PDE sensitivity equation method for optimal aerodynamic design. *Journal of Computational Physics*, 136(2):366–384, 1997.

[5] A. Brault, L. Dumas and D. Lucor. Uncertainty quantification of inflow boundary condition and proximal arterial stiffness coupled effect on pulse wave propagation in a vascular network. *International Journal for Numerical Methods in Biomedical Engineering*, 2016.

[6] A. Canestrelli, A. Siviglia, M. Dumbser, and E. F. Toro. Well-Balanced High-Order Centred Schemes for Non-Conservative Hyperbolic Systems. Application to Shallow Water with Fixed and Mobile Bed. *Advances in Water Resource*, 32(6):834–844, 2009.

[7] C. Cattaneo, Sulla conduzione del calore, Atti Semin. Mat. Fis. Univ. Modena 3 (1948) 83-101.

[8] C. Chalons, R. Duvigneau, C. Fiorini. Sensitivity analysis and numerical diffusion effects for hyperbolic PDE systems with discontinuous solutions. The case of barotropic Euler equations in Lagrangian coordinates. *SIAM Journal on Scientific Computing, Society for Industrial and Applied Mathematics*, 2018, 40 (6), pp.A3955-A3981.

[9] C. Chalons, J.-F. Coulombel. Relaxation approximation of the Euler equations. *J. Math. Anal. Appl.*, vol 348, no. 2, pp 872-893 (2008)

[10] C. Chalons, C. Fiorini, R. Duvigneau. Sensitivity equation method for Euler equations in presence of shocks applied to uncertainty quantification. 2018. hal-01817815

[11] F. Chifari. Sensitivity analysis for hyperbolic equations with application to blood flow in elastic arteries. University of Trento, Italy, July 2018.

[12] C. M. Colciago, S. Deparis, A. Quarteroni. Comparisons between reduced order models and full 3D models for fluid-structure interaction problems in haemodynamics. *Journal of Computational and Applied Mathematics*, 265(2014)120-138.
[13] S. F. Davis. Simplified Second-Order Godunov-Type Methods. *SIAM J. Sci. Stat. Comput.*, 9:445-473, 1988.

[14] L. Dumas, T. El Bouti and D. Lucor. A robust and subject-specific hemodynamic model of the lower limb based on non-invasive arterial measurement.

[15] R. Duvigneau and D. Pelletier. A sensitivity equation method for fast evaluation of nearby flows and uncertainty analysis for shape parameters. *Int. J. of CFD*, 20(7):497-512, August 2006.

[16] R. Duvigneau, D. Pelletier, and J. Borggaard. An improved continuous sensitivity equation method for optimal shape design in mixed convection. *Numerical Heat Transfer part B : Fundamentals*, 50(1):1-24, July 2006.

[17] C. Fiorini. Sensitivity analysis for nonlinear hyperbolic equations of conservation laws. Université Paris-Saclay, 2018, PhD.

[18] L. Formaggia, A. Quarteroni, A. Veneziani. Cardiovascular Mathematics: Modeling and Simulation of the Circulatory System. Springer-Verlag: Italia, Milano, 2009.

[19] L. Grinberg, E. Cheever, T. Anor, J. R. Madsen and G. E. Karniadakis. Modeling Blood Flow Circulation in Intracranial Arterial Networks: A Comparative 3D/1D Simulation Study. *Annals of Biomed Eng.*, 39(1):297-309, 2010.

[20] M. D. Gunzburger. Perspectives in flow control and optimization, volume 5. Siam, 2003.

[21] T. Y. Hou and P. G. Le Floch. Why Nonconservative Schemes Converge to Wrong Solutions: Error Analysis. *Mathematics of computation*, volume 62, number 206, 1994, pages 497-530.

[22] S. Jin and Z. Xin. The relaxation schemes for systems of conservation laws in arbitrary space dimension. *Comm. Pure. Appl. Math.*, 48:235-276, 1995.

[23] A. Manzoni, A. Quarteroni, C. Vergara. The Cardiovascular System: Mathematical Modeling, Numerical Algorithms, Clinical Applications. 2016.

[24] K. S. Matthys, J. Alastruey, A. W. Khir, P. Segers, S. J. Sherwin, P. R. Verdonck, K. H. Parker, J. Peiro. Pulse wave propagation in a model human arterial network: Assessment of 1-D numerical simulations against in vitro measurements. *Journal of Biomechanics*, 40,3476-3486, 2007.

[25] G. I. Montecinos, L. O. Müller, and E. F. Toro. Hyperbolic reformulation of a 1D visco-elastic blood flow model and ADER finite volume schemes. *Journal of Computational Physics*, 266:101-123, 2014.

[26] L. O. Müller, G. Leugering, and P. J. Blanco. Consistent treatment of visco-elastic effects at junctions in one-dimensional blood flow models. *Journal of Computational Physics*, 314:167-193, 2016.

[27] L. O. Müller and E. F. Toro. A global multiscale mathematical model for the human circulation with emphasis on the venous system. *International Journal for Numerical Methods in Biomedical Engineering*, 2014.

[28] C. Parés. Numerical Methods for Nonconservative Hyperbolic Systems: a Theoretical Framework. *SIAM Journal on Numerical Analysis*, 44:300-321, 2006.

[29] M. Petrella, S. Tokareva, and E. F. Toro. Uncertainty Quantification for Hyperbolic Systems with application to blood flow in arteries. *Journal of Computational Physics*, 386:405-427, 2019.

[30] I. Suliciu. On the thermodynamics of fluids with relaxation and phase transitions. Fluids with relaxation. *Int. J. Engag. Sci.*, 36:921-947, 1998.

[31] E. F. Toro and S. J. Billett. Centered TVD schemes for hyperbolic conservation laws. *IMA Journal of Numerical Analysis*, (2000) 20, 47-79.

[32] E. F. Toro, R. C. Millington and L. A. M. Nejad. Towards very high order to schemes. In Toro, E. F. (ed.), *Godunov Methods: Theory and Applications*, Kluwer Academic/ Plenum Publishers, New York, pp. 907-940. 2001.
[33] E. F. Toro and G. I. Montecinos. Advection-diffusion-reaction equations: hyperbolization and high order ADER discretizations. *SIAM Journal of Scientific Computing*, 36(5):A2423-A2457, 2014.

[34] E. F. Toro and A. Siviglia. PRICE: Primitive Centred Schemes for Hyperbolic Systems. *Int. J. Numer. Meth. in Fluids*, 42:1263-1291, 2003.

[35] E. F. Toro, M. Spruce, and W. Speares. Restoration of the Contact Surface in the HLL-Riemann Solver. Technical Report CoA-9204, Department of Aerospace Science, Colleague of Aeronautics, Cranfield Institute of Technology, UK, 1992.

[36] E. F. Toro. Brain venous haemodynamics, neurological diseases and mathematical modelling. A review. *Applied Mathematics and Computation*, 272 (2016) 542-579.

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

[38] E. F. Toro. *Riemann Solvers and Numerical Methods for Fluid Dynamics, Third Edition*. , Springer-Verlag, 2009.

[39] N. Xiao, J. Alastruey and F. C. Alberto. A systematic comparison between 1-D and 3-D hemodynamics in compliant arterial models. *Int. J. Numer. Meth. Biomed. Engng.*, 2014; 30:204-231.