# Reduced Basis method for finite volume simulations of parabolic PDEs applied to porous media flows

Jana Tarhini, Sébastien Boyaval, Guillaume Enchéry, Quang Huy Tran

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

## To cite this version:

Jana Tarhini, Sébastien Boyaval, Guillaume Enchéry, Quang Huy Tran. Reduced Basis method for finite volume simulations of parabolic PDEs applied to porous media flows. 2024. hal-04608663

HAL Id: hal-04608663

https://enpc.hal.science/hal-04608663

Preprint submitted on 11 Jun 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Reduced Basis method for finite volume simulations of parabolic PDEs
applied to porous media flows

Jana TARHINI*

Sébastien BOYAVAL†

Guillaume ENCHÉRY*

Quang-Huy TRAN*

June 11, 2024

Abstract

Numerical simulations are a highly valuable tool to evaluate the impact of the uncertainties of various model parameters, and to optimize e.g. injection-production scenarios in the context of underground storage (of CO₂ typically). Finite volume approximations of Darcy's parabolic model for flows in porous media are typically run many times, for many values of parameters like permeability and porosity, at costly computational efforts.

We study the relevance of reduced basis methods as a way to lower the overall simulation cost of finite volume approximations to Darcy's parabolic model for flows in porous media for different values of the parameters such as permeability. In the context of underground gas storage (of CO₂ typically) in saline aquifers, our aim is to evaluate quickly, for many parameter values, the flux along some interior boundaries near the well injection area—regarded as a quantity of interest—. To this end, we construct reduced bases by a standard POD-Greedy algorithm. Our POD-Greedy algorithm uses a new goal-oriented error estimator designed from a discrete space-time energy norm independent of the parameter. We provide some numerical experiments that validate the efficiency of the proposed estimator.

Keywords

single-phase flow, porous media, finite volumes, reduced basis, goal-oriented error estimate

Mathematics subject classification

35J50, 65M08, 65N15, 76S05

# 1 Introduction

In the context of geological storage of gases such as CO2, computational models of single phase Darcy flow are useful to optimize the efficiency of injection, and to quantify uncertainties with a view to assessing the enduring stability of the storage site. As concerns CO2, it is usually injected in underground storage sites such as depleted oil and gas reservoirs or saline aquifers in sedimentary basins. In this work, we are mainly interested in the case of saline aquifers.

In saline aquifers, the numerical simulation of *single phase* Darcy flows is very meaningful, in particular in a large domain at basin scale where it is computationally costly. Indeed, brine is moved by the gas (CO2) injected outside the storage area. In risk assessment studies, one needs to evaluate the pressure field in the surrounding aquifer (typically along faults far from the storage domain) *many times, for many values of the uncertain parameters*. Quantifying the impact of the uncertainties of model parameters, on the time evolution of the flux at underground boundaries of the storage area, is also desired to optimize the injection process. For both purposes and given values of the model parameters, the flow simulator computes the solution of a large linear system many times. This multi-query setting induces costly computational efforts especially for large domains. To lower the overall simulation cost of computations at many parameter values, we consider a Reduced Basis (RB) approach.

Many other methods have been proposed to reduce the time calculations in basin modeling and reservoir simulation. For instance, in case of LGR methods, the grid is locally refined only in a region of interest depending on the

<footnote>*IFP Energies nouvelles, 1 et 4 avenue de Bois Préau, 92852 Rueil-Malmaison Cedex, France. jana.tarhini@ifpen.fr, guillaume.enchery@ifpen.fr, quang-huy.tran@ifpen.fr</footnote>

<footnote>†Laboratoire d'hydraulique Saint-Venant, École des Ponts, EDF R&D, 6 quai Watier, 78401 Chatou Cedex, France & Materials, Inria, Paris, France. sebastien.boyaval@enpc.fr</footnote>
Reduced Basis method for finite volume simulations of parabolic PDEs
applied to porous media flows

Jana TARHINI*

Sébastien BOYAVAL†

Guillaume ENCHÉRY*

Quang-Huy TRAN*

June 11, 2024

Abstract

Numerical simulations are a highly valuable tool to evaluate the impact of the uncertainties of various model parameters, and to optimize e.g. injection-production scenarios in the context of underground storage (of CO₂ typically). Finite volume approximations of Darcy’s parabolic model for flows in porous media are typically run many times, for many values of parameters like permeability and porosity, at costly computational efforts.

We study the relevance of reduced basis methods as a way to lower the overall simulation cost of finite volume approximations to Darcy’s parabolic model for flows in porous media for different values of the parameters such as permeability. In the context of underground gas storage (of CO₂ typically) in saline aquifers, our aim is to evaluate quickly, for many parameter values, the flux along some interior boundaries near the well injection area—regarded as a quantity of interest—. To this end, we construct reduced bases by a standard POD-Greedy algorithm. Our POD-Greedy algorithm uses a new goal-oriented error estimator designed from a discrete space-time energy norm independent of the parameter. We provide some numerical experiments that validate the efficiency of the proposed estimator.

Keywords

single-phase flow, porous media, finite volumes, reduced basis, goal-oriented error estimate

Mathematics subject classification

35J50, 65M08, 65N15, 76S05

# 1 Introduction

In the context of geological storage of gases such as CO2, computational models of single phase Darcy flow are useful to optimize the efficiency of injection, and to quantify uncertainties with a view to assessing the enduring stability of the storage site. As concerns CO2, it is usually injected in underground storage sites such as depleted oil and gas reservoirs or saline aquifers in sedimentary basins. In this work, we are mainly interested in the case of saline aquifers.

In saline aquifers, the numerical simulation of *single phase* Darcy flows is very meaningful, in particular in a large domain at basin scale where it is computationally costly. Indeed, brine is moved by the gas (CO2) injected outside the storage area. In risk assessment studies, one needs to evaluate the pressure field in the surrounding aquifer (typically along faults far from the storage domain) *many times, for many values of the uncertain parameters*. Quantifying the impact of the uncertainties of model parameters, on the time evolution of the flux at underground boundaries of the storage area, is also desired to optimize the injection process. For both purposes and given values of the model parameters, the flow simulator computes the solution of a large linear system many times. This multi-query setting induces costly computational efforts especially for large domains. To lower the overall simulation cost of computations at many parameter values, we consider a Reduced Basis (RB) approach.

Many other methods have been proposed to reduce the time calculations in basin modeling and reservoir simulation. For instance, in case of LGR methods, the grid is locally refined only in a region of interest depending on the

<footnote>*IFP Energies nouvelles, 1 et 4 avenue de Bois Préau, 92852 Rueil-Malmaison Cedex, France. jana.tarhini@ifpen.fr, guillaume.enchery@ifpen.fr, quang-huy.tran@ifpen.fr
†Laboratoire d'hydraulique Saint-Venant, École des Ponts, EDF R&D, 6 quai Watier, 78401 Chatou Cedex, France & Materials, Inria, Paris, France. sebastien.boyaval@enpc.fr</footnote>
local solution properties, whereas in other areas, where the solution is relatively smooth or uniform, grid cells can
be larger, leading to a smaller computational cost. Likewise, the Adaptive Mesh Refinement (AMR) method [4, 29]
divides the computational domain into a hierarchy of grids and each grid is refined or coarsened by considering an
error estimate as the simulation progresses.

In this work, our approach consists rather in considering a reduced basis (RB) approach to replace many calls
to a parametrized High-Fidelity (HF) simulator at many parameter values, by calls to a less expensive Low-Fidelity
(LF) surrogate model with a certification of the error.

A RB procedure relies on an existing computational model, a parametrized HF simulator which can provide one
with numerical approximations of the model solutions at fixed parameter values. When the HF model consists in
large linear systems (one at each time step in a nonstationary flow simulation e.g.), a LF (reduced) computational
model is usually constructed by Galerkin projection of the HF model onto a linear subspace.

In the context of porous media flows, the choice of the appropriate numerical scheme to discretize the governing
equations in space is crucial to obtain a consistent approximation of the fluxes. Finite-difference [5, 26], finite-
volume [15] or finite-element [11, 14] methods have been classically used in industrial contexts such as reservoir
engineering. In particular the finite-volume two-point flux approximation is a reference method in this field because
of its simplicity and its stability properties (the discrete operator turns out to be an M-matrix). However, once the
grids are not Λ-orthogonal <fnref n="1" />, this scheme is no more consistent. Over the past years, new discretization methods
have been proposed to satisfy this property: multi-point flux approximations [1], mimetic finite-differences [6],
virtual elements [3], hybrid [16] or vertex-centred finite-volumes [17] to quote just a few of them. We also mention
non-linear schemes [23, 24, 27] that were designed in order to obtain monotone approximations and properties such
as the positivity of the solutions or the maximum principle on these grids. In this work, we consider the average
multi-point flux approximation (MPFA-FV) method which was for instance studied in [27]. That approximation
does not preserve the positivity of the solutions or the maximum principle, but it is consistent on grids that are not
Λ-orthogonal.

In the present work, given a parametrized HF model that discretizes Darcy’s parabolic model by a MPFA-FV
method resulting in a time-series of (large) linear systems, see Section 2, we then standardly construct a parametrized
LF model by projecting (the pressure field solution to) Darcy flows at each time step, whatever the parameter value,
onto one (single) linear subspace by Galerkin method.

To that aim, we adopt a standard two-stage procedure [18, 21]. First, in a costly offline stage, we identify a
linear approximation space spanned by (snapshots of) simulations at relevant parameter values. During this stage,
HF simulations with many degrees of freedom $N$ (several thousands of cells) are run at least $N$ times, $N$ being
the dimension of the linear approximation space. During the offline stage, a LF reduced computational model is
also numerically constructed. Next, in an online stage, the values of solutions and the quantities of interests at
yet-unexplored parameter values are evaluated numerically using the LF (reduced) computational model, if possible
with a computational complexity independent from $N$.

The offline selection of a good linear subspace for Galerkin projection is crucial to the quality of the LF model,
i.e. to control the approximation error of the LF model with respect to the HF model at every parameter values.
Regarding the applications of the RB method to parametrized HF simulator that are based on (finite-volume
discretizations of) parabolic PDEs, a standard selection technique is the POD-Greedy method based on a reliable
a posteriori error estimation, see e.g. [21].

In the present work, we propose new a posteriori error estimators. Like e.g. [21], our a posteriori error estimators
evaluate the approximation by a LF model of a HF model discretizing parabolic PDEs by the finite volume method (a
MPFA-FV discretization in our application case to single phase Darcy flows). However, by contrast with standard
RB literature, we use a discrete space-time energy norm of $L^2([0, T]; H^1(\Omega))$-type that is *independent from the
parameter*. We also provide goal-oriented estimators for a linear QOI, like in [18, 20] (where a parametrized HF
model based on finite element approximations is reduced using another algorithm than POD-Greedy to construct
the LF model).

To assess the accuracy of our error estimate, and the performance of our new certified RB approach (i.e. the
dimension of the reduced LF model in a multi-query scenario for monophasic Darcy flows parametrized by the
permeability), we also perform numerical simulations.

Our HF model based on MPFA-FV discretization is non-affine in the permeability parameter: we use the
Empirical Interpolation Method (EIM) [10, 25] for the construction of a LF model independent from $N$. Moreover,

<footnote n="1">¹A grid is Λ-orthogonal if the product of the permeability tensor Λ with the face normal is orthogonal to the line joining the cell and face centers.</footnote>
local solution properties, whereas in other areas, where the solution is relatively smooth or uniform, grid cells can
be larger, leading to a smaller computational cost. Likewise, the Adaptive Mesh Refinement (AMR) method [4, 29]
divides the computational domain into a hierarchy of grids and each grid is refined or coarsened by considering an
error estimate as the simulation progresses.

In this work, our approach consists rather in considering a reduced basis (RB) approach to replace many calls
to a parametrized High-Fidelity (HF) simulator at many parameter values, by calls to a less expensive Low-Fidelity
(LF) surrogate model with a certification of the error.

A RB procedure relies on an existing computational model, a parametrized HF simulator which can provide one
with numerical approximations of the model solutions at fixed parameter values. When the HF model consists in
large linear systems (one at each time step in a nonstationary flow simulation e.g.), a LF (reduced) computational
model is usually constructed by Galerkin projection of the HF model onto a linear subspace.

In the context of porous media flows, the choice of the appropriate numerical scheme to discretize the governing
equations in space is crucial to obtain a consistent approximation of the fluxes. Finite-difference [5, 26], finite-
volume [15] or finite-element [11, 14] methods have been classically used in industrial contexts such as reservoir
engineering. In particular the finite-volume two-point flux approximation is a reference method in this field because
of its simplicity and its stability properties (the discrete operator turns out to be an M-matrix). However, once the
grids are not Λ-orthogonal <fnref n="1" />, this scheme is no more consistent. Over the past years, new discretization methods
have been proposed to satisfy this property: multi-point flux approximations [1], mimetic finite-differences [6],
virtual elements [3], hybrid [16] or vertex-centred finite-volumes [17] to quote just a few of them. We also mention
non-linear schemes [23, 24, 27] that were designed in order to obtain monotone approximations and properties such
as the positivity of the solutions or the maximum principle on these grids. In this work, we consider the average
multi-point flux approximation (MPFA-FV) method which was for instance studied in [27]. That approximation
does not preserve the positivity of the solutions or the maximum principle, but it is consistent on grids that are not
Λ-orthogonal.

In the present work, given a parametrized HF model that discretizes Darcy’s parabolic model by a MPFA-FV
method resulting in a time-series of (large) linear systems, see Section 2, we then standardly construct a parametrized
LF model by projecting (the pressure field solution to) Darcy flows at each time step, whatever the parameter value,
onto one (single) linear subspace by Galerkin method.

To that aim, we adopt a standard two-stage procedure [18, 21]. First, in a costly offline stage, we identify a
linear approximation space spanned by (snapshots of) simulations at relevant parameter values. During this stage,
HF simulations with many degrees of freedom $N$ (several thousands of cells) are run at least $N$ times, $N$ being
the dimension of the linear approximation space. During the offline stage, a LF reduced computational model is
also numerically constructed. Next, in an online stage, the values of solutions and the quantities of interests at
yet-unexplored parameter values are evaluated numerically using the LF (reduced) computational model, if possible
with a computational complexity independent from $N$.

The offline selection of a good linear subspace for Galerkin projection is crucial to the quality of the LF model,
i.e. to control the approximation error of the LF model with respect to the HF model at every parameter values.
Regarding the applications of the RB method to parametrized HF simulator that are based on (finite-volume
discretizations of) parabolic PDEs, a standard selection technique is the POD-Greedy method based on a reliable
a posteriori error estimation, see e.g. [21].

In the present work, we propose new a posteriori error estimators. Like e.g. [21], our a posteriori error estimators
evaluate the approximation by a LF model of a HF model discretizing parabolic PDEs by the finite volume method (a
MPFA-FV discretization in our application case to single phase Darcy flows). However, by contrast with standard
RB literature, we use a discrete space-time energy norm of $L^2([0, T]; H^1(\Omega))$-type that is *independent from the
parameter*. We also provide goal-oriented estimators for a linear QOI, like in [18, 20] (where a parametrized HF
model based on finite element approximations is reduced using another algorithm than POD-Greedy to construct
the LF model).

To assess the accuracy of our error estimate, and the performance of our new certified RB approach (i.e. the
dimension of the reduced LF model in a multi-query scenario for monophasic Darcy flows parametrized by the
permeability), we also perform numerical simulations.

Our HF model based on MPFA-FV discretization is non-affine in the permeability parameter: we use the
Empirical Interpolation Method (EIM) [10, 25] for the construction of a LF model independent from $N$. Moreover,

<footnote n="1">¹A grid is Λ-orthogonal if the product of the permeability tensor Λ with the face normal is orthogonal to the line joining the cell and face centers.</footnote>
to construct a lower bound of the coercivity constant (required by a rigorous error estimation) we use the Successive Constraint Method (SCM) [12, 22]. It is noteworthy that, for accurate a posteriori error estimation taking care of machine precision, we numerically compute quadratic-in-the-parameter forms (in the dual norm of the residual) as in [7], using a specific orthonormal basis to represent the residual, while a traditional offline/online decomposition [8] leads to numerical precision issues.

Our numerical results show that the proposed a posteriori error estimation guarantees a reliable evaluation of a single-phase Darcy flow model and accurately quantify the solution error. We also proved that the goal-oriented estimation for a given linear QOI offers for a small dimension of the reduced model a very good precision for the output error.

This paper is outlined as follows. In Section 2, we present a discretization of single-phase flow (SPF) equations in porous media based on the multi-point flux approximation that defines our HF numerical model. Section 3 discusses the concept of a goal-oriented method applied to the SPF problem. We derive the a posteriori error estimation for the primal and dual problems as well as for the output model. We also present the POD-Greedy algorithm to construct the reduced basis. In Section 3.5, we additionally elaborate on the computation of the residual dual norm in order to avoid the impact of round-off errors on the error bound. We then numerically study the behavior and efficiency of the proposed estimate in Section 4. Finally, concluding remarks are given in Section 5.

# 2 A parametrized High-Fidelity model

In this section, we introduce the PDE modelling of single phase porous media flows, and its discretization by a finite volume method which defines our parametrized HF simulator in the sequel.

## 2.1 A Darcy model of single phase porous media flows

We consider the flow of a slightly compressible fluid saturating a porous rock within a connected and bounded polygonal domain $\Omega$ of $\mathbb{R}^3$ and a time $T > 0$. The boundary $\partial\Omega = \Gamma_D \cup \Gamma_N$ of $\Omega$ is partitioned into a part where Dirichlet boundary conditions are applied, and a part where homogeneous Neumann boundary conditions are used.

The balance of the water volume combined with Darcy's law and with initial and boundary data leads to

$$ \phi c_t \partial_t p - \nabla \cdot (\Lambda(\nabla p + \rho g \nabla z)) = q, \quad \text{in } (0, T) \times \Omega, \quad (2.1a) $$

$$ \Lambda(\nabla p + \rho g \nabla z) \cdot \mathbf{n} = 0, \quad \text{on } (0, T) \times \Gamma_N, \quad (2.1b) $$

$$ p = p_D, \quad \text{on } (0, T) \times \Gamma_D, \quad (2.1c) $$

$$ p(x, t = 0) = p^0(x), \quad \text{in } \Omega, \quad (2.1d) $$

where $p$ denotes the fluid pressure, $\Lambda = \bar{\kappa}/\mu$ the mobility tensor, $\bar{\kappa}$ the rock permeability tensor, $\mu$ the fluid viscosity, $\phi$ the rock porosity, $c_t$ the total compressibility, $\rho$ the fluid density, $g$ the gravity constant and $q(p)$ a well source term to be precised later (in Section 2.2 after discretization). We designate by $\mathbf{n}$ the unit normal vector outside the domain.

A typical multi-query setting (with parameter variations to be addressed by the RB method), occurs when the permeability tensor $\bar{\kappa}$ is uncertain. We here assume that the domain is made up of two areas, corresponding to two rock types (a reservoir one and a cap rock) each with constant isotropic permeabilities $\bar{\kappa}_1$ or $\bar{\kappa}_2$ (see e.g. Fig. 2.1) so that

$$ \Lambda = \begin{bmatrix} \Lambda & 0 & 0 \\ 0 & \Lambda & 0 \\ 0 & 0 & \Lambda \end{bmatrix} $$

where $\Lambda(x) = \sum_{i=1}^2 \frac{\kappa_i}{\mu} 1_{(i)}(x)$ and $1_i$ denotes the indicator function of the region $i$.
to construct a lower bound of the coercivity constant (required by a rigorous error estimation) we use the Successive Constraint Method (SCM) [12, 22]. It is noteworthy that, for accurate a posteriori error estimation taking care of machine precision, we numerically compute quadratic-in-the-parameter forms (in the dual norm of the residual) as in [7], using a specific orthonormal basis to represent the residual, while a traditional offline/online decomposition [8] leads to numerical precision issues.

Our numerical results show that the proposed a posteriori error estimation guarantees a reliable evaluation of a single-phase Darcy flow model and accurately quantify the solution error. We also proved that the goal-oriented estimation for a given linear QOI offers for a small dimension of the reduced model a very good precision for the output error.

This paper is outlined as follows. In Section 2, we present a discretization of single-phase flow (SPF) equations in porous media based on the multi-point flux approximation that defines our HF numerical model. Section 3 discusses the concept of a goal-oriented method applied to the SPF problem. We derive the a posteriori error estimation for the primal and dual problems as well as for the output model. We also present the POD-Greedy algorithm to construct the reduced basis. In Section 3.5, we additionally elaborate on the computation of the residual dual norm in order to avoid the impact of round-off errors on the error bound. We then numerically study the behavior and efficiency of the proposed estimate in Section 4. Finally, concluding remarks are given in Section 5.

# 2 A parametrized High-Fidelity model

In this section, we introduce the PDE modelling of single phase porous media flows, and its discretization by a finite volume method which defines our parametrized HF simulator in the sequel.

## 2.1 A Darcy model of single phase porous media flows

We consider the flow of a slightly compressible fluid saturating a porous rock within a connected and bounded polygonal domain $\Omega$ of $\mathbb{R}^3$ and a time $T > 0$. The boundary $\partial\Omega = \Gamma_D \cup \Gamma_N$ of $\Omega$ is partitioned into a part where Dirichlet boundary conditions are applied, and a part where homogeneous Neumann boundary conditions are used. The balance of the water volume combined with Darcy's law and with initial and boundary data leads to

$$ \phi c_t \partial_t p - \nabla \cdot (\Lambda(\nabla p + \rho g \nabla z)) = q, \qquad \text{in } (0, T) \times \Omega, \qquad (2.1a) $$

$$ \Lambda(\nabla p + \rho g \nabla z) \cdot \mathbf{n} = 0, \qquad \text{on } (0, T) \times \Gamma_N, \qquad (2.1b) $$

$$ p = p_D, \qquad \text{on } (0, T) \times \Gamma_D, \qquad (2.1c) $$

$$ p(x, t = 0) = p^0(x), \qquad \text{in } \Omega, \qquad (2.1d) $$

where $p$ denotes the fluid pressure, $\Lambda = \bar{\kappa}/\mu$ the mobility tensor, $\bar{\kappa}$ the rock permeability tensor, $\mu$ the fluid viscosity, $\phi$ the rock porosity, $c_t$ the total compressibility, $\rho$ the fluid density, $g$ the gravity constant and $q(p)$ a well source term to be precised later (in Section 2.2 after discretization). We designate by $\mathbf{n}$ the unit normal vector outside the domain.

A typical multi-query setting (with parameter variations to be addressed by the RB method), occurs when the permeability tensor $\bar{\kappa}$ is uncertain. We here assume that the domain is made up of two areas, corresponding to two rock types (a reservoir one and a cap rock) each with constant isotropic permeabilities $\bar{\kappa}_1$ or $\bar{\kappa}_2$ (see e.g. Fig. 2.1) so that

$$ \Lambda = \begin{bmatrix} \Lambda & 0 & 0 \\ 0 & \Lambda & 0 \\ 0 & 0 & \Lambda \end{bmatrix} $$

where $\Lambda(x) = \sum_{i=1}^2 \frac{\kappa_i}{\mu} 1_{(i)}(x)$ and $1_i$ denotes the indicator function of the region $i$.
<figure><img src="image_1.png" /><figcaption>Figure 2.1: Typical domain configuration.</figcaption></figure>

Throughout this work, we consider the storage area (S) with boundaries

$$ \Gamma_{\text{int}} = \{y_1\} \times [0, z_1] \cup \{y_2\} \times [0, z_1] \cup [y_1, y_2] \times \{z_1\}, $$

where we seek to predict the time evolution of the flux s defined by

$$ s = - \int_{\Gamma_{\text{int}}} \mathbf{\Lambda}(\nabla p + \rho g \nabla z) \cdot \mathbf{n} \, dS \qquad (2.2) $$

over $\Gamma_{\text{int}}$ for many values of $\kappa_1$ and $\kappa_2$.

## 2.2 Finite volume discretization

The single-phase flow equations are first discretized in time using the implicit Euler method, next in space using a finite volume method. Choosing a constant time step

$$ \Delta t = \frac{T}{N}, \quad N \in \mathbb{N}^*, $$

we consider at each time iteration $n \in \{0, \dots, N-1\}$ an approximation $p^{n+1} \approx p(t^{n+1})$ at $t^{n+1} = (n+1) \Delta t$ solution to

$$ \phi c_t \frac{p^{n+1} - p^n}{\Delta t} - \nabla \cdot (\mathbf{\Lambda}(\nabla p^{n+1} + \rho g \nabla z)) = q^{n+1}, \qquad (2.3a) $$

$$ \mathbf{\Lambda}(\nabla p^{n+1} + \rho g \nabla z) \cdot \mathbf{n} = 0, \qquad (2.3b) $$

$$ p^{n+1} = p_D. \qquad (2.3c) $$

The space discretization is performed using an admissible mesh of $\Omega$, defined by a triplet $\mathcal{D} = (\mathcal{T}, \mathcal{E}, \mathcal{P})$ where:

* $\mathcal{T}$ is a finite set of non empty compact convex polygonal sub-domains of $\Omega$ (the set of cells), called control volumes such that $\bar{\Omega} = \bigcup_{K \in \mathcal{T}} \bar{K}$. For all $K \in \mathcal{T}$, we denote by $m_K > 0$ its measure and set $\partial K \stackrel{\text{def}}{=} \bar{K} \setminus K$.

* $\mathcal{E}$ is a family of subsets of $\bar{\Omega}$ (the set of faces) such that for any $K \in \mathcal{T}$, there exists a subset $\mathcal{E}_K$ of $\mathcal{E}$ where $\partial K = \bigcup_{\sigma \in \mathcal{E}_K} \sigma$. For any $(K, L) \in \mathcal{T}^2$ with $K \neq L$, either the $(d-1)$ Lebesgue measure of $\bar{K} \cap \bar{L}$ is 0 or $\bar{K} \cap \bar{L} = \bar{\sigma}$ for some $\sigma \in \mathcal{E}$, with $\sigma = K|L$ (an interior face). We denote by $m_\sigma$ the $(d-1)$-dimensional measure of $\sigma$. The sets of inner and boundary faces are denoted by $\mathcal{E}_{\text{int}}$ and $\mathcal{E}_{\text{ext}}$ respectively.

* $\mathcal{P} = \{\mathbf{x}_K\}_{K \in \mathcal{T}}$ is a collection of points within $\Omega$ indexed by $\mathcal{T}$ (called the cell centers, not required to be the barycenters) s.t. $\mathbf{x}_K \in K$ and $K$ is star-shaped with respect to $\mathbf{x}_K$.

For each cell $K \in \mathcal{T}$ and face $\sigma \in \mathcal{E}_K$, $\mathbf{n}_{K,\sigma}$ denotes the unit vector normal to $\sigma$ and pointing outward to $K$. Additionally, for any cell $K \in \mathcal{T}$ and any function $\Phi$ belonging to $L^1(K)$, we define $\langle \Phi \rangle_K \stackrel{\text{def}}{=} m_K^{-1} \int_K \Phi \, dx$.

Given an admissible mesh, numerically computable approximations $p_K^n \approx \langle p^n \rangle_K$ are defined after space discretization by a finite volume method. We first integrate (2.3a) over a cell $K$ to obtain

$$ \int_K \phi c_t \frac{p^{n+1} - p^n}{\Delta t} \, dx - \int_K \nabla \cdot (\mathbf{\Lambda}(\nabla p^{n+1} + \rho g \nabla z)) \, dx = \int_K q^{n+1} \, dx. \qquad (2.4) $$
<figure><img src="image_3.png" /><figcaption>Figure 2.1: Typical domain configuration.</figcaption></figure>

Throughout this work, we consider the storage area (S) with boundaries

$$ \Gamma_{\text{int}} = \{y_1\} \times [0, z_1] \cup \{y_2\} \times [0, z_1] \cup [y_1, y_2] \times \{z_1\}, $$

where we seek to predict the time evolution of the flux s defined by

$$ s = - \int_{\Gamma_{\text{int}}} \mathbf{\Lambda}(\nabla p + \rho g \nabla z) \cdot \mathbf{n} \, dS \qquad (2.2) $$

over $\Gamma_{\text{int}}$ for many values of $\kappa_1$ and $\kappa_2$.

## 2.2 Finite volume discretization

The single-phase flow equations are first discretized in time using the implicit Euler method, next in space using a finite volume method. Choosing a constant time step

$$ \Delta t = \frac{T}{N}, \quad N \in \mathbb{N}^*, $$

we consider at each time iteration $n \in \{0, \dots, N-1\}$ an approximation $p^{n+1} \approx p(t^{n+1})$ at $t^{n+1} = (n+1) \Delta t$ solution to

$$ \phi c_t \frac{p^{n+1} - p^n}{\Delta t} - \nabla \cdot (\mathbf{\Lambda}(\nabla p^{n+1} + \rho g \nabla z)) = q^{n+1}, \qquad (2.3a) $$

$$ \mathbf{\Lambda}(\nabla p^{n+1} + \rho g \nabla z) \cdot \mathbf{n} = 0, \qquad (2.3b) $$

$$ p^{n+1} = p_{\text{D}}. \qquad (2.3c) $$

The space discretization is performed using an admissible mesh of $\Omega$, defined by a triplet $\mathcal{D} = (\mathcal{T}, \mathcal{E}, \mathcal{P})$ where:

* $\mathcal{T}$ is a finite set of non empty compact convex polygonal sub-domains of $\Omega$ (the set of cells), called control volumes such that $\bar{\Omega} = \bigcup_{K \in \mathcal{T}} \bar{K}$. For all $K \in \mathcal{T}$, we denote by $m_K > 0$ its measure and set $\partial K \stackrel{\text{def}}{=} \bar{K} \setminus K$.

* $\mathcal{E}$ is a family of subsets of $\bar{\Omega}$ (the set of faces) such that for any $K \in \mathcal{T}$, there exists a subset $\mathcal{E}_K$ of $\mathcal{E}$ where $\partial K = \bigcup_{\sigma \in \mathcal{E}_K} \sigma$. For any $(K, L) \in \mathcal{T}^2$ with $K \neq L$, either the $(d-1)$ Lebesgue measure of $\bar{K} \cap \bar{L}$ is 0 or $\bar{K} \cap \bar{L} = \bar{\sigma}$ for some $\sigma \in \mathcal{E}$, with $\sigma = K|L$ (an interior face). We denote by $m_\sigma$ the $(d-1)$-dimensional measure of $\sigma$. The sets of inner and boundary faces are denoted by $\mathcal{E}_{\text{int}}$ and $\mathcal{E}_{\text{ext}}$ respectively.

* $\mathcal{P} = \{\mathbf{x}_K\}_{K \in \mathcal{T}}$ is a collection of points within $\Omega$ indexed by $\mathcal{T}$ (called the cell centers, not required to be the barycenters) s.t. $\mathbf{x}_K \in K$ and $K$ is star-shaped with respect to $\mathbf{x}_K$.

For each cell $K \in \mathcal{T}$ and face $\sigma \in \mathcal{E}_K$, $\mathbf{n}_{K,\sigma}$ denotes the unit vector normal to $\sigma$ and pointing outward to $K$. Additionally, for any cell $K \in \mathcal{T}$ and any function $\Phi$ belonging to $L^1(K)$, we define $\langle \Phi \rangle_K \stackrel{\text{def}}{=} m_K^{-1} \int_K \Phi \, dx$.

Given an admissible mesh, numerically computable approximations $p_K^n \approx \langle p^n \rangle_K$ are defined after space discretization by a finite volume method. We first integrate (2.3a) over a cell $K$ to obtain

$$ \int_K \phi c_t \frac{p^{n+1} - p^n}{\Delta t} \, dx - \int_K \nabla \cdot (\mathbf{\Lambda}(\nabla p^{n+1} + \rho g \nabla z)) \, dx = \int_K q^{n+1} \, dx. \qquad (2.4) $$
Applying Green's formula, we can transform the first two integrals and recast (2.4) as

$$ m_K \phi_K c_t \frac{\langle p^{n+1} \rangle_K - \langle p^n \rangle_K}{\Delta t} - \int_{\partial K} \Lambda(\nabla p^{n+1} + \rho g \nabla z) \cdot \mathbf{n}_K \, d\gamma = \int_K q^{n+1} \, dx. \quad (2.5) $$

By decomposing the boundary $\partial K$ into faces, we get

$$ m_K c_t \phi_K \frac{\langle p^{n+1} \rangle_K - \langle p^n \rangle_K}{\Delta t} - \sum_{\sigma \in \mathcal{E}_K} \int_{\sigma} \Lambda(\nabla p^{n+1} + \rho g \nabla z) \cdot \mathbf{n}_{K,\sigma} \, dS = \int_K q^{n+1} \, dx, $$

which leads one to consider the following numerical scheme

$$ m_K \phi_K c_t (p_K^{n+1} - p_K^n) + \Delta t \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma}^{n+1} = \Delta t m_K q_K^{n+1} \quad (2.6) $$

with numerical fluxes $F_{K,\sigma}^{n+1} \approx -\int_{\sigma} \Lambda(\nabla p^{n+1} + \rho g \nabla z) \cdot \mathbf{n}_{K,\sigma} \, dS$ and numerical source terms $q_K^{n+1}$ that allow for the numerical computation of $(p_K^{n+1})_{K \in \mathcal{T}}$ given $(p_K^n)_{K \in \mathcal{T}}$.

The Peaceman model [13] is used here for the well source term $q_K^{n+1} = m_K^{-1} \int_K q^{n+1} \, dx$. We suppose that the well is vertical and the perforations are oriented in the z-direction. The well model is then given by

$$ q_K^{n+1} = W I_K (p_{bh} - p_K^{n+1} - \rho g (z_{bh} - z_K)), \quad (2.7) $$

where $p_{bh}$ is the bottom hole pressure and $W I_K$ is the Peaceman well index<fnref n="2" /> in a perforated cell $K$.

The flux $-\int_{\sigma} \Lambda(\nabla p^{n+1} + \rho g \nabla z) \cdot \mathbf{n}_{K,\sigma} \, dS$ is numerically approximated using the average multi-point flux scheme studied in [27]. For each interior edge $\sigma \in \mathcal{E}_{\text{int}}$, with $\mathcal{T}_{\sigma} = \{K, L\}$ the approximated flux $F_{K,\sigma}^{n+1}$ is defined as a convex combination of two linear fluxes $\tilde{F}_{K,\sigma}^{n+1}$ and $\tilde{F}_{L,\sigma}^{n+1}$ such that

$$ F_{K,\sigma}^{n+1} = \mu_{K,\sigma} \tilde{F}_{K,\sigma}^{n+1} - \mu_{L,\sigma} \tilde{F}_{L,\sigma}^{n+1}, \quad \text{with } \mu_{K,\sigma} \ge 0, \mu_{L,\sigma} \ge 0, \mu_{K,\sigma} + \mu_{L,\sigma} = 1. \quad (2.8) $$

A numerical flux formula as given by (2.8) is clearly conservative, i.e,

$$ F_{K,\sigma}^{n+1} + F_{L,\sigma}^{n+1} = 0. \quad (2.9) $$

To build the linear fluxes $\tilde{F}_{K,\sigma}^{n+1}$ in (2.8), we approximate the pressure gradient $\nabla p$ in the direction of the conormal vector $\langle \Lambda \rangle_K \mathbf{n}_{K,\sigma}$ after expressing the conormal as a linear combination of the vectors $(\mathbf{x}_{\sigma'} - \mathbf{x}_K)_{\{\sigma' \in S_{K,\sigma}\}}$

$$ \langle \Lambda \rangle_K \mathbf{n}_{K,\sigma} \approx \sum_{\sigma' \in S_{K,\sigma}} \alpha_{K,\sigma\sigma'} (\mathbf{x}_{\sigma'} - \mathbf{x}_K). \quad (2.10) $$

The decomposition (2.10) is achieved numerically by means of an optimization procedure which aims at reducing the sum of the coefficients $\alpha_{K,\sigma\sigma'}$ and the size of the stencil $S_{K,\sigma} \stackrel{\text{def}}{=} \{\sigma' \in \mathcal{E}_K | \alpha_{K,\sigma\sigma'} \neq 0\}$ within $\mathcal{E}_K$ [28]. In (2.10), $\mathbf{x}_{\sigma}$ is not the face center but a harmonic averaging interpolation point

$$ \mathbf{x}_{\sigma} = \omega_{K,\sigma} \mathbf{y}_K + \omega_{L,\sigma} \mathbf{y}_L + \frac{d_{K,\sigma} d_{L,\sigma}}{d_{L,\sigma} \tau_{K,\sigma} + d_{K,\sigma} \tau_{L,\sigma}} (\tau_K^{\sigma} - \tau_L^{\sigma}) \quad (2.11) $$

$$ W I = \frac{2\pi h_3 \sqrt{\lambda_1 \lambda_2}}{\ln(r_e/r_w) + s_d} $$

when

$$ \mathbf{A} = \begin{bmatrix} \lambda_1 & 0 & 0 \\ 0 & \lambda_2 & 0 \\ 0 & 0 & \lambda_3 \end{bmatrix} $$

$h_3$ is the perforation height, $r_w$ is the well radius, $s_d$ is the skin factor i.e. a dimensionless number modeling the formation damage caused by drilling, and $r_e$ is the Peaceman radius defined as

$$ r_e = \frac{0.14[(\lambda_2/\lambda_1)^{1/2}h_1^2 + (\lambda_1/\lambda_2)^{1/2}h_2^2]^{1/2}}{0.5[(\lambda_2/\lambda_1)^{1/4} + (\lambda_1/\lambda_2)^{1/4}]}, $$

where $h_1$ and $h_2$ are the grid sizes in x and y directions.

<footnote n="2">²We adopt here a usual definition of the well index [13]</footnote>
Applying Green's formula, we can transform the first two integrals and recast (2.4) as

$$ m_K \phi_K c_t \frac{\langle p^{n+1} \rangle_K - \langle p^n \rangle_K}{\Delta t} - \int_{\partial K} \mathbf{\Lambda}(\nabla p^{n+1} + \rho g \nabla z) \cdot \mathbf{n}_K \, d\gamma = \int_K q^{n+1} \, dx. \quad (2.5) $$

By decomposing the boundary $\partial K$ into faces, we get

$$ m_K c_t \phi_K \frac{\langle p^{n+1} \rangle_K - \langle p^n \rangle_K}{\Delta t} - \sum_{\sigma \in \mathcal{E}_K} \int_{\sigma} \mathbf{\Lambda}(\nabla p^{n+1} + \rho g \nabla z) \cdot \mathbf{n}_{K,\sigma} \, dS = \int_K q^{n+1} \, dx, $$

which leads one to consider the following numerical scheme

$$ m_K \phi_K c_t (p_K^{n+1} - p_K^n) + \Delta t \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma}^{n+1} = \Delta t m_K q_K^{n+1} \quad (2.6) $$

with numerical fluxes $F_{K,\sigma}^{n+1} \approx -\int_{\sigma} \mathbf{\Lambda}(\nabla p^{n+1} + \rho g \nabla z) \cdot \mathbf{n}_{K,\sigma} \, dS$ and numerical source terms $q_K^{n+1}$ that allow for the numerical computation of $(p_K^{n+1})_{K \in \mathcal{T}}$ given $(p_K^n)_{K \in \mathcal{T}}$.

The Peaceman model [13] is used here for the well source term $q_K^{n+1} = m_K^{-1} \int_K q^{n+1} \, dx$. We suppose that the well is vertical and the perforations are oriented in the z-direction. The well model is then given by

$$ q_K^{n+1} = \mathbf{W}\mathbf{I}_K(p_{bh} - p_K^{n+1} - \rho g(z_{bh} - z_K)), \quad (2.7) $$

where $p_{bh}$ is the bottom hole pressure and $\mathbf{W}\mathbf{I}_K$ is the Peaceman well index<fnref n="2" /> in a perforated cell $K$.

The flux $-\int_{\sigma} \mathbf{\Lambda}(\nabla p^{n+1} + \rho g \nabla z) \cdot \mathbf{n}_{K,\sigma} \, dS$ is numerically approximated using the average multi-point flux scheme studied in [27]. For each interior edge $\sigma \in \mathcal{E}_{\text{int}}$, with $\mathcal{T}_{\sigma} = \{K, L\}$ the approximated flux $F_{K,\sigma}^{n+1}$ is defined as a convex combination of two linear fluxes $\tilde{F}_{K,\sigma}^{n+1}$ and $\tilde{F}_{L,\sigma}^{n+1}$ such that

$$ F_{K,\sigma}^{n+1} = \mu_{K,\sigma} \tilde{F}_{K,\sigma}^{n+1} - \mu_{L,\sigma} \tilde{F}_{L,\sigma}^{n+1}, \quad \text{with } \mu_{K,\sigma} \ge 0, \mu_{L,\sigma} \ge 0, \mu_{K,\sigma} + \mu_{L,\sigma} = 1. \quad (2.8) $$

A numerical flux formula as given by (2.8) is clearly conservative, i.e,

$$ F_{K,\sigma}^{n+1} + F_{L,\sigma}^{n+1} = 0. \quad (2.9) $$

To build the linear fluxes $\tilde{F}_{K,\sigma}^{n+1}$ in (2.8), we approximate the pressure gradient $\nabla p$ in the direction of the conormal vector $\langle \mathbf{\Lambda} \rangle_K \mathbf{n}_{K,\sigma}$ after expressing the conormal as a linear combination of the vectors $(\mathbf{x}_{\sigma'} - \mathbf{x}_K)_{\{\sigma' \in S_{K,\sigma}\}}$

$$ \langle \mathbf{\Lambda} \rangle_K \mathbf{n}_{K,\sigma} \approx \sum_{\sigma' \in S_{K,\sigma}} \alpha_{K,\sigma\sigma'} (\mathbf{x}_{\sigma'} - \mathbf{x}_K). \quad (2.10) $$

The decomposition (2.10) is achieved numerically by means of an optimization procedure which aims at reducing the sum of the coefficients $\alpha_{K,\sigma\sigma'}$ and the size of the stencil $S_{K,\sigma} \stackrel{\text{def}}{=} \{\sigma' \in \mathcal{E}_K \alpha_{K,\sigma\sigma'} \neq 0\}$ within $\mathcal{E}_K$ [28]. In (2.10), $\mathbf{x}_{\sigma}$ is not the face center but a harmonic averaging interpolation point

$$ \mathbf{x}_{\sigma} = \omega_{K,\sigma} \mathbf{y}_K + \omega_{L,\sigma} \mathbf{y}_L + \frac{d_{K,\sigma} d_{L,\sigma}}{d_{L,\sigma} \tau_{K,\sigma} + d_{K,\sigma} \tau_{L,\sigma}} (\tau_K^{\sigma} - \tau_L^{\sigma}) \quad (2.11) $$

$$ \mathbf{W}\mathbf{I} = \frac{2\pi h_3 \sqrt{\lambda_1 \lambda_2}}{\ln(r_e/r_w) + s_d} $$

when

$$ \mathbf{\Lambda} = \begin{bmatrix} \lambda_1 & 0 & 0 \\ 0 & \lambda_2 & 0 \\ 0 & 0 & \lambda_3 \end{bmatrix} $$

$h_3$ is the perforation height, $r_w$ is the well radius, $s_d$ is the skin factor i.e. a dimensionless number modeling the formation damage caused by drilling, and $r_e$ is the Peaceman radius defined as

$$ r_e = \frac{0.14[(\lambda_2/\lambda_1)^{1/2}h_1^2 + (\lambda_1/\lambda_2)^{1/2}h_2^2]^{1/2}}{0.5[(\lambda_2/\lambda_1)^{1/4} + (\lambda_1/\lambda_2)^{1/4}]}, $$

where $h_1$ and $h_2$ are the grid sizes in x and y directions.

<footnote n="2">²We adopt here a usual definition of the well index [13]</footnote>
where

$$ \omega_{K,\sigma} = \frac{d_{L,\sigma}\tau_{K,\sigma}}{d_{L,\sigma}\tau_{K,\sigma} + d_{K,\sigma}\tau_{L,\sigma}}, \quad \omega_{L,\sigma} = \frac{d_{K,\sigma}\tau_{L,\sigma}}{d_{L,\sigma}\tau_{K,\sigma} + d_{K,\sigma}\tau_{L,\sigma}}, \qquad (2.12) $$

$$ \tau_{K,\sigma} = \mathbf{n}_{K,\sigma}\langle\mathbf{\Lambda}\rangle_K \mathbf{n}_{K,\sigma}, \quad \tau_{L,\sigma} = \mathbf{n}_{L,\sigma}\langle\mathbf{\Lambda}\rangle_L \mathbf{n}_{L,\sigma}, \qquad (2.13) $$

$$ \boldsymbol{\tau}_{K}^{\sigma} = (\mathbf{\Lambda}_{K} - \tau_{K,\sigma}\mathrm{Id})\mathbf{n}_{K,\sigma}, \quad \boldsymbol{\tau}_{L}^{\sigma} = (\mathbf{\Lambda}_{L} - \tau_{L,\sigma}\mathrm{Id})\mathbf{n}_{L,\sigma}, \qquad (2.14) $$

$d_{K,\sigma}$, $d_{L,\sigma}$ are the distances of the cell centers to $\sigma$, $\mathbf{y}_K$, $\mathbf{y}_L$ their projection on $\sigma$ defined by (see Figure 2.2)

$$ \mathbf{y}_K = \mathbf{x}_K + d_{K,\sigma}\mathbf{n}_{K,\sigma}, \quad \mathbf{y}_L = \mathbf{x}_L + d_{L,\sigma}\mathbf{n}_{L,\sigma}. $$

<figure><img src="image_2.png" /><figcaption>Figure 2.2: Harmonic averaging point.</figcaption></figure>

The pressure trace at $\sigma \in \mathcal{E}_K$ is consistently reconstructed as

$$ I_{\sigma}p = \sum_{M \in \{K,L\}} \omega_{M,\sigma}p_{M}, \quad \text{where} \quad \sum_{M \in \{K,L\}} \omega_{M,\sigma} = 1, \quad \omega_{M,\sigma} \ge 0, \qquad (2.15) $$

using the same weights $\omega_{M,\sigma}$ as in (2.12) (for more details see [2]). With the previous approximations, the linear fluxes read

$$ \tilde{F}_{K,\sigma}^{n+1} = m_{\sigma} \sum_{\sigma' \in S_{K,\sigma}} \alpha_{K,\sigma\sigma'} [(p^{n+1} + \rho g z)_{K,\sigma'} - (p_K^{n+1} + \rho g z_K)], \qquad (2.16) $$

where we used the notation

$$ (u)_{K,\sigma} = \begin{cases} I_{\sigma}u & \text{if } \sigma = K|L, \\ u_{\sigma} & \text{otherwise.} \end{cases} \qquad (2.17) $$

Finally, to define the numerical flux $F_{K,\sigma}^{n+1}$, the weights $\mu_{K,\sigma}$ and $\mu_{L,\sigma}$ can be chosen in different ways. In this work, we set $\mu_{K,\sigma} = \mu_{L,\sigma} = \frac{1}{2}$ if $\sigma = K|L$ and $\mu_{K,\sigma} = 1$ if $\sigma \subset \partial\Omega$. Let us remark that in (2.16), the average value of the boundary condition (2.1c) is used on the face $\sigma'$ if $\sigma' \subset \Gamma_D$. An additional unknown $p_{\sigma'}$ is added if $\sigma' \subset \Gamma_N$ which can be computed by adding the homogeneous Neumann condition to the discrete system, namely:

$$ F_{K,\sigma}^{n+1} = 0. $$

Our numerical output quantity of interest (QOI) $s$ is then defined as

$$ s^{n+1} = \sum_{\substack{\sigma=K|L \\ \sigma \subset \Gamma_{\text{int}}}} F_{K,\sigma}^{n+1}. \qquad (2.18) $$

Our HF model consists in the solutions of the $N$ discrete linear systems obtained after assembling equations (2.6) for all $K \in \mathcal{T}$, $n \in \{0, \dots, N-1\}$

$$ (M + \Delta t \mathbf{A}) \mathbf{p}_{M}^{n+1} = M \mathbf{p}_{M}^{n} + \Delta t \mathbf{b}, \qquad (2.19) $$

see also (3.1) below, where

* $\mathbf{A} \in \mathbb{R}^{N \times N}$ is a matrix containing the terms $\alpha_{K,\sigma\sigma'}$, with $N = N_c + N_b$ ($N_c$ is the total number of cells and $N_b$ the number of boundaries where we have imposed Neumann boundary condition),
where

$$ \omega_{K,\sigma} = \frac{d_{L,\sigma}\tau_{K,\sigma}}{d_{L,\sigma}\tau_{K,\sigma} + d_{K,\sigma}\tau_{L,\sigma}}, \quad \omega_{L,\sigma} = \frac{d_{K,\sigma}\tau_{L,\sigma}}{d_{L,\sigma}\tau_{K,\sigma} + d_{K,\sigma}\tau_{L,\sigma}}, \qquad (2.12) $$

$$ \tau_{K,\sigma} = \mathbf{n}_{K,\sigma} \langle \mathbf{\Lambda} \rangle_K \mathbf{n}_{K,\sigma}, \quad \tau_{L,\sigma} = \mathbf{n}_{L,\sigma} \langle \mathbf{\Lambda} \rangle_L \mathbf{n}_{L,\sigma}, \qquad (2.13) $$

$$ \boldsymbol{\tau}_{\mathrm{K}}^{\sigma} = (\boldsymbol{\Lambda}_{\mathrm{K}} - \tau_{\mathrm{K},\sigma} \mathrm{Id})\boldsymbol{n}_{\mathrm{K},\sigma}, \quad \boldsymbol{\tau}_{\mathrm{L}}^{\sigma} = (\boldsymbol{\Lambda}_{\mathrm{L}} - \tau_{\mathrm{L},\sigma} \mathrm{Id})\boldsymbol{n}_{\mathrm{L},\sigma}, \qquad (2.14) $$

$d_{K,\sigma}$, $d_{L,\sigma}$ are the distances of the cell centers to $\sigma$, $\mathbf{y}_K$, $\mathbf{y}_L$ their projection on $\sigma$ defined by (see Figure 2.2)

$$ \mathbf{y}_K = \mathbf{x}_K + d_{K,\sigma}\mathbf{n}_{K,\sigma}, \quad \mathbf{y}_L = \mathbf{x}_L + d_{L,\sigma}\mathbf{n}_{L,\sigma}. $$

<figure><img src="image_4.png" /><figcaption>Figure 2.2: Harmonic averaging point.</figcaption></figure>

The pressure trace at $\sigma \in \mathcal{E}_K$ is consistently reconstructed as

$$ I_{\sigma}p = \sum_{M \in \{K,L\}} \omega_{M,\sigma}p_{M}, \quad \text{where} \quad \sum_{M \in \{K,L\}} \omega_{M,\sigma} = 1, \quad \omega_{M,\sigma} \ge 0, \qquad (2.15) $$

using the same weights $\omega_{M,\sigma}$ as in (2.12) (for more details see [2]). With the previous approximations, the linear fluxes read

$$ \tilde{F}_{K,\sigma}^{n+1} = m_{\sigma} \sum_{\sigma' \in S_{K,\sigma}} \alpha_{K,\sigma\sigma'} [(p^{n+1} + \rho g z)_{K,\sigma'} - (p_K^{n+1} + \rho g z_K)], \qquad (2.16) $$

where we used the notation

$$ (u)_{K,\sigma} = \begin{cases} I_{\sigma}u & \text{if } \sigma = K|L, \\ u_{\sigma} & \text{otherwise.} \end{cases} \qquad (2.17) $$

Finally, to define the numerical flux $F_{K,\sigma}^{n+1}$, the weights $\mu_{K,\sigma}$ and $\mu_{L,\sigma}$ can be chosen in different ways. In this work, we set $\mu_{K,\sigma} = \mu_{L,\sigma} = \frac{1}{2}$ if $\sigma = K|L$ and $\mu_{K,\sigma} = 1$ if $\sigma \subset \partial\Omega$. Let us remark that in (2.16), the average value of the boundary condition (2.1c) is used on the face $\sigma'$ if $\sigma' \subset \Gamma_D$. An additional unknown $p_{\sigma'}$ is added if $\sigma' \subset \Gamma_N$ which can be computed by adding the homogeneous Neumann condition to the discrete system, namely:

$$ F_{K,\sigma}^{n+1} = 0. $$

Our numerical output quantity of interest (QOI) $s$ is then defined as

$$ s^{n+1} = \sum_{\substack{\sigma=K|L \\ \sigma \subset \Gamma_{\text{int}}}} F_{K,\sigma}^{n+1}. \qquad (2.18) $$

Our HF model consists in the solutions of the $N$ discrete linear systems obtained after assembling equations (2.6) for all $K \in \mathcal{T}$, $n \in \{0, \dots, N-1\}$

$$ (M + \Delta t \mathbf{A})\mathbf{p}_{\mathrm{M}}^{n+1} = \mathbf{M}\mathbf{p}_{\mathrm{M}}^{n} + \Delta t \mathbf{b}, \qquad (2.19) $$

see also (3.1) below, where

* $\mathbf{A} \in \mathbb{R}^{N \times N}$ is a matrix containing the terms $\alpha_{K,\sigma\sigma'}$, with $\mathcal{N} = N_c + N_b$ ($N_c$ is the total number of cells and $N_b$ the number of boundaries where we have imposed Neumann boundary condition),
* $M \in \mathbb{R}^{N \times N}$ is such that

$$M = \begin{bmatrix} M_c & 0 \\ 0 & 0 \end{bmatrix},$$

where $M_c \in \mathbb{R}^{N_c \times N_c}$ is a diagonal matrix made of the quantities $m_K \phi_K c_t$,

* $p_M^{n+1} \in \mathbb{R}^N$ is a vector composed of the values of the pressure $p^{n+1}$ in each element $K \in \mathcal{T}$ and on the edges $\sigma \subset \Gamma_N \cap \partial K$,

* $b \in \mathbb{R}^N$ contains the Dirichlet condition values as well as the source terms $Q_K^{n+1}$.

The QOI (2.18) can be rewritten in

$$s^{n+1} = l^T p_M^{n+1} + c, \qquad (2.20)$$

with $l \in \mathbb{R}^N$. We want to reduce the computational cost associated with the multiple resolutions of equations (2.6) and (2.20) that occur when changing the permeability values $\kappa_1$ and $\kappa_2$ and assembling the corresponding values of $A$, $b$ and $l$. We equip $\mathbb{R}^N$ with the inner product $\langle \cdot, \cdot \rangle_{G^*}$ and the corresponding norm $\| \cdot \|_{G^*}$, where $G^*$ is a symmetric positive definite matrix in $\mathbb{R}^{N \times N}$ that will be defined in the sequel. The Euclidean norm on $\mathbb{R}^N$ is denoted by $\| \cdot \|$.

# 3 Reduced Basis technique in time-dependent setting

We now develop a goal-oriented RB procedure for HF model (2.19) and the QOI (2.20).

In Section 3.1, we define a LF model based on Galerkin projection. Then, in Section 3.2, we explain how to compute a Galerkin projection subspace by a POD-Greedy approach. The crux of our POD-Greedy approach is the minimization of an a posteriori error estimator of the QOI that is described in Section 3.4.

## 3.1 Low-fidelity model

We consider the discrete space-time problem

$$\begin{bmatrix} M + \Delta t A & 0 & 0 & \dots & 0 \\ -M & M + \Delta t A & \dots & \dots & \dots \\ 0 & -M & \dots & \dots & \dots \\ \vdots & \vdots & \ddots & \ddots & \ddots \\ 0 & \dots & \dots & 0 & -M \quad M + \Delta t A \end{bmatrix} \begin{bmatrix} p_M^1 \\ p_M^2 \\ \vdots \\ p_M^{N-1} \\ p_M^N \end{bmatrix} = \Delta t \begin{bmatrix} b + M p_M^0 \\ b \\ \vdots \\ b \\ b \end{bmatrix}, \quad (3.1)$$

in a multi-query setting where (3.1) has to be solved for many values of a parameter contained in $A$ and $b$. In the sequel, we generically denote $\xi$ that parameter – which is simply $(\kappa_1, \kappa_2) \in \mathbb{R}^2$ in our application.

Numerical reduction (i.e. reduction of the computational cost without loss of accuracy) can be achieved by replacing (3.1), for many parameter values, by its Galerkin projection onto a linear subspace spanned by a reduced basis $Z_{pr} \in \mathbb{R}^{N \times N_{pr}}$, $N_{pr} \ll N$. The reduced solution $p^{N_{pr},n+1}$ at time $n+1$ is defined as

$$p^{N_{pr},n+1} = Z_{pr} \tilde{p}^{n+1}, \qquad (3.2)$$

where $\tilde{p}^{n+1}$ is the solution to

$$(Z_{pr}^T M Z_{pr} + \Delta t Z_{pr}^T A Z_{pr}) \tilde{p}^{n+1} = Z_{pr}^T M Z_{pr} \tilde{p}^n + \Delta t Z_{pr}^T b \quad (3.3)$$

in the so-called **online phase**, once $Z_{pr}$ has been computed. To that aim, a reduced basis $Z_{pr}$ can first be computed in a so-called **offline phase**, e.g. using a POD-Greedy method.
* $M \in \mathbb{R}^{N \times N}$ is such that

$$M = \begin{bmatrix} M_c & 0 \\ 0 & 0 \end{bmatrix},$$

where $M_c \in \mathbb{R}^{N_c \times N_c}$ is a diagonal matrix made of the quantities $m_K \phi_K c_t$,

* $p_M^{n+1} \in \mathbb{R}^N$ is a vector composed of the values of the pressure $p^{n+1}$ in each element $K \in \mathcal{T}$ and on the edges $\sigma \subset \Gamma_N \cap \partial K$,

* $b \in \mathbb{R}^N$ contains the Dirichlet condition values as well as the source terms $Q_K^{n+1}$.

The QOI (2.18) can be rewritten in

$$s^{n+1} = l^T p_M^{n+1} + c, \qquad (2.20)$$

with $l \in \mathbb{R}^N$. We want to reduce the computational cost associated with the multiple resolutions of equations (2.6) and (2.20) that occur when changing the permeability values $\kappa_1$ and $\kappa_2$ and assembling the corresponding values of $A$, $b$ and $l$. We equip $\mathbb{R}^N$ with the inner product $\langle \cdot, \cdot \rangle_{G^*}$ and the corresponding norm $\| \cdot \|_{G^*}$, where $G^*$ is a symmetric positive definite matrix in $\mathbb{R}^{N \times N}$ that will be defined in the sequel. The Euclidean norm on $\mathbb{R}^N$ is denoted by $\| \cdot \|$.

# 3 Reduced Basis technique in time-dependent setting

We now develop a goal-oriented RB procedure for HF model (2.19) and the QOI (2.20).

In Section 3.1, we define a LF model based on Galerkin projection. Then, in Section 3.2, we explain how to compute a Galerkin projection subspace by a POD-Greedy approach. The crux of our POD-Greedy approach is the minimization of an a posteriori error estimator of the QOI that is described in Section 3.4.

## 3.1 Low-fidelity model

We consider the discrete space-time problem

$$ \begin{bmatrix} M + \Delta t A & 0 & 0 & \dots & 0 \\ -M & M + \Delta t A & \dots & \dots & \dots \\ 0 & -M & \dots & \dots & \dots \\ \vdots & \vdots & \ddots & \ddots & \ddots \\ 0 & \dots & \dots & 0 & -M \end{bmatrix} \begin{bmatrix} p_M^1 \\ p_M^2 \\ \vdots \\ p_M^{N-1} \\ p_M^N \end{bmatrix} = \Delta t \begin{bmatrix} b + M p_M^0 \\ b \\ \vdots \\ b \\ b \end{bmatrix}, \quad (3.1) $$

in a multi-query setting where (3.1) has to be solved for many values of a parameter contained in $A$ and $b$. In the sequel, we generically denote $\xi$ that parameter – which is simply $(\kappa_1, \kappa_2) \in \mathbb{R}^2$ in our application.

Numerical reduction (i.e. reduction of the computational cost without loss of accuracy) can be achieved by replacing (3.1), for many parameter values, by its Galerkin projection onto a linear subspace spanned by a reduced basis $Z_{pr} \in \mathbb{R}^{N \times N_{pr}}$, $N_{pr} \ll N$. The reduced solution $p^{N_{pr},n+1}$ at time $n+1$ is defined as

$$p^{N_{pr},n+1} = Z_{pr} \tilde{p}^{n+1}, \qquad (3.2)$$

where $\tilde{p}^{n+1}$ is the solution to

$$(Z_{pr}^T M Z_{pr} + \Delta t Z_{pr}^T A Z_{pr}) \tilde{p}^{n+1} = Z_{pr}^T M Z_{pr} \tilde{p}^n + \Delta t Z_{pr}^T b \qquad (3.3)$$

in the so-called **online phase**, once $Z_{pr}$ has been computed. To that aim, a reduced basis $Z_{pr}$ can first be computed in a so-called **offline phase**, e.g. using a POD-Greedy method.
## 3.2 POD-Greedy method

The aim of a POD-Greedy algorithm is to iteratively construct matrices $Z_{pr}^{N_1} \in \mathbb{R}^{N \times N_1}$, $Z_{pr}^{N_2} \in \mathbb{R}^{N \times N_2}$, ..., of rank $N_1 < N_2 < ...$ such that Span $Z_{pr}^{N_i}$ is a subspace of the vector space Span $Z_{pr}^{N_j}$ spanned by the column vectors of $Z_{pr}^{N_j}$ as soon as $N_i \le N_j$. At the end of the algorithm, $Z_{pr}^{N_{pr}} \equiv Z_{pr}$ is the so-called reduced basis useful for Galerkin projection in (3.3), of dimension $N_{pr} \ll N$.

The final iteration is reached when some projection error is under a fixed tolerance $\epsilon_{tol} > 0$, for instance

$$ \|e^n\|_{pr}(\xi) \le \epsilon_{tol} \quad \forall \xi \in \Xi, $$

where $\|\cdot\|_{pr}$ is a norm on the HF space $\mathbb{R}^N$ (see proposition 3.1),

$$ e^n(\xi) = p_M^n(\xi) - p^{N_{pr},n}(\xi) \quad (3.4) $$

is the "primal" reduction error at parameter value $\xi$, and $\Xi := \{\xi_1, \dots, \xi_L\}$ is a training set of parameter values. Between two iterations, a scalar $ric \in [0, 1]$ controls the increase $N_{i+1} - N_i$: one only adds to $Z_{pr}^{N_i} \in \mathbb{R}^{N \times N_i}$ the largest $N_{i+1} - N_i$ POD modes of a "snapshot" matrix (collecting the time evolution of the projection error at a new selected parameter $\xi_\ell$) (see Algorithm 1). This type of basis-increase between two iterations has been introduced in [21] and has remained a standard since, when one iteratively constructs a reduced basis with a greedy-type algorithm that iteratively selects parameter values $\xi_1, \xi_2, \dots \in \Xi$ so as to increase $Z_{pr}^{N_i} \in \mathbb{R}^{N \times N_i}$ using the time trajectory $\{p_M^{n+1}(\xi_\ell)\}_{n=0}^{N-1}$ at iteration $i$.

The quality of the POD-Greedy selection (i.e. the accuracy reached by the approximation $p^{N_{pr},n}(\xi) \approx p_M^n(\xi)$ for all $n = 1, \dots, N$ and $\xi \in \Xi$) strongly depends on the parameter $\xi_\ell$ selected at iteration $i$. To that aim, we propose to choose $\xi_\ell$ as a maximizer of a new a posteriori estimator $\Delta_{pr}^N(\xi)$ of the reduction error

$$ \|\mathbf{p}_{\mathcal{M}}^{N}(\xi) - \mathbf{p}^{N_{\text{pr}},N}(\xi)\|_{\text{pr}} \le \Delta_{\text{pr}}^{N}(\xi) $$

among all training parameter values $\xi \in \Xi$.

**Algorithm 1** POD-greedy algorithm using $\Delta_{pr}^N$

1: **Procedure** $Z_{pr} = \text{POD-Greedy}(N_{\text{max}}, \epsilon_{\text{tol}}, \Xi, ric)$.

2: $N_{pr} = 1$, $\delta^{N_{pr}} = \epsilon_{\text{tol}} + 1$.

3: Take $\xi_1 \in \Xi$, $\ell = 1$ and set $\Xi^\ell = \{\xi_1\}$.

4: Define $Z_{pr} = \emptyset$.

5: **while** $\delta^{N_{pr}} > \epsilon_{\text{tol}}$ and $N_{pr} < N_{\text{max}}$ **do**.

6: Compute $p_M^n(\xi_\ell)$ for $1 \le n \le N$.

7: Set $S_{pr} := [p_M^1(\xi_\ell) - \text{Proj}_{Z_{pr}}(p_M^1(\xi_\ell))|\dots|p_M^N(\xi_\ell) - \text{Proj}_{Z_{pr}}(p_M^N(\xi_\ell))]$.

8: Compute $[z_1|\dots|z_{\delta N_{pr}}] = \text{POD}(S_{pr}, ric)$ using Algorithm 2.

9: Define $Z_{pr}^{N_{pr}+\delta N_{pr}} := \text{orthonormalize}(Z_{pr}^{N_{pr}} \cup [z_1|\dots|z_{\delta N_{pr}}])$ using Algorithm 3.

10: Compute $\delta^{N_{pr}} = \max_{\xi \in \Xi} \Delta_{pr}^N$.

11: Set $\xi_{\ell+1} = \arg \max_{\xi \in \Xi} \Delta_{pr}^N$.

12: $\Xi^{\ell+1} \leftarrow \Xi^\ell \cup \{\xi_{\ell+1}\}$.

13: $N_{pr} \leftarrow N_{pr} + \delta N_{pr}$.

14: $\ell \leftarrow \ell + 1$.

15: **end while**
## 3.2 POD-Greedy method

The aim of a POD-Greedy algorithm is to iteratively construct matrices $Z_{pr}^{N_1} \in \mathbb{R}^{N_1 \times N_1}$, $Z_{pr}^{N_2} \in \mathbb{R}^{N_2 \times N_2}$, ..., of rank $N_1 < N_2 < ...$ such that $\text{Span} Z_{pr}^{N_i}$ is a subspace of the vector space $\text{Span} Z_{pr}^{N_j}$ spanned by the column vectors of $Z_{pr}^{N_j}$ as soon as $N_i \le N_j$. At the end of the algorithm, $Z_{pr}^{N_{pr}} \equiv Z_{pr}$ is the so-called reduced basis useful for Galerkin projection in (3.3), of dimension $N_{pr} \ll N$.

The final iteration is reached when some projection error is under a fixed tolerance $\epsilon_{tol} > 0$, for instance

$$ \|e^n\|_{pr}(\xi) \le \epsilon_{tol} \quad \forall \xi \in \Xi, $$

where $\|\cdot\|_{pr}$ is a norm on the HF space $\mathbb{R}^N$ (see proposition 3.1),

$$ e^n(\xi) = p_M^n(\xi) - p^{N_{pr},n}(\xi) \quad (3.4) $$

is the "primal" reduction error at parameter value $\xi$, and $\Xi := \{\xi_1, \dots, \xi_L\}$ is a training set of parameter values. Between two iterations, a scalar $ric \in [0, 1]$ controls the increase $N_{i+1} - N_i$: one only adds to $Z_{pr}^{N_i} \in \mathbb{R}^{N_i \times N_i}$ the largest $N_{i+1} - N_i$ POD modes of a "snapshot" matrix (collecting the time evolution of the projection error at a new selected parameter $\xi_\ell$) (see Algorithm 1). This type of basis-increase between two iterations has been introduced in [21] and has remained a standard since, when one iteratively constructs a reduced basis with a greedy-type algorithm that iteratively selects parameter values $\xi_\ell, \xi_{\ell+1}, \dots \in \Xi$ so as to increase $Z_{pr}^{N_i} \in \mathbb{R}^{N_i \times N_i}$ using the time trajectory $\{p_M^{n+1}(\xi_\ell)\}_{n=0}^{N_i-1}$ at iteration $i$.

The quality of the POD-Greedy selection (i.e. the accuracy reached by the approximation $p^{N_{pr},n}(\xi) \approx p_M^n(\xi)$ for all $n = 1, \dots, N$ and $\xi \in \Xi$) strongly depends on the parameter $\xi_\ell$ selected at iteration $i$. To that aim, we propose to choose $\xi_\ell$ as a maximizer of a new a posteriori estimator $\Delta_{pr}^N(\xi)$ of the reduction error

$$ \|\mathbf{p}_M^N(\xi) - \mathbf{p}^{N_{pr},N}(\xi)\|_{pr} \le \Delta_{pr}^N(\xi) $$

among all training parameter values $\xi \in \Xi$.

**Algorithm 1** POD-greedy algorithm using $\Delta_{pr}^N$

1: **Procedure** $Z_{pr} = \text{POD-Greedy}(N_{max}, \epsilon_{tol}, \Xi, ric)$.

2: $N_{pr} = 1$, $\delta^{N_{pr}} = \epsilon_{tol} + 1$.

3: Take $\xi_1 \in \Xi$, $\ell = 1$ and set $\Xi^\ell = \{\xi_1\}$.

4: Define $Z_{pr} = \emptyset$.

5: **while** $\delta^{N_{pr}} > \epsilon_{tol}$ and $N_{pr} < N_{max}$ **do**.

6: Compute $p_M^n(\xi_\ell)$ for $1 \le n \le N$.

7: Set $\mathbf{S}_{pr} := [p_M^1(\xi_\ell) - \text{Proj}_{Z_{pr}}(p_M^1(\xi_\ell))|\dots|p_M^N(\xi_\ell) - \text{Proj}_{Z_{pr}}(p_M^N(\xi_\ell))]$.

8: Compute $[\mathbf{z}_1|\dots|\mathbf{z}_{\delta N_{pr}}] = \text{POD}(\mathbf{S}_{pr}, ric)$ using Algorithm 2.

9: Define $Z_{pr}^{N_{pr}+\delta N_{pr}} := \text{orthonormalize}(Z_{pr}^{N_{pr}} \cup [\mathbf{z}_1|\dots|\mathbf{z}_{\delta N_{pr}}])$ using Algorithm 3.

10: Compute $\delta^{N_{pr}} = \max_{\xi \in \Xi} \Delta_{pr}^N$.

11: Set $\xi_{\ell+1} = \arg \max_{\xi \in \Xi} \Delta_{pr}^N$.

12: $\Xi^{\ell+1} \leftarrow \Xi^\ell \cup \{\xi_{\ell+1}\}$.

13: $N_{pr} \leftarrow N_{pr} + \delta N_{pr}$.

14: $\ell \leftarrow \ell + 1$.

15: **end while**
**Algorithm 2 POD with ric**

Input: $S_{pr} \in \mathbb{R}^{N \times N}$, ric $\in (0, 1)$

Output: $G^*$-orthonormal $[z_1|...|z_{\delta N_{pr}}]$.

1: **for** i $\leftarrow 1, ..., N$ **do**

2: $V_i$, $\sigma_i$ i largest vector and eigenvalue pair of $C = S_{pr}^T G^* S_{pr}$

3: **end for**

4: $\delta N_{pr} \leftarrow 1$

$$ \sum^{\delta N_{pr}} \lambda_n $$

5: **while** $E_{\delta N_{pr}} = \frac{\sum_{n=1}^{n=1} \lambda_n}{N} < ric$ **do**

$$ \sum^{\delta N_{pr}} \lambda_n $$

6: $\delta N_{pr} \leftarrow 1 + \delta N_{pr}$

7: $z_{\delta N_{pr}} = \frac{1}{\sqrt{\sigma_{\delta N_{pr}}}} S_{pr} \tilde{V}_{\delta N_{pr}}$

8: **end while**

**Algorithm 3** Gram-Schmidt with re-iteration

Input: vectors $v_i$, $i \in 1, ..., N_{pr}$.

Output: orthonormal vectors $v_i$.

1: **for** $\ell = 1, 2$ **do**

2: **for** $i \leftarrow 1, ..., N_{pr}$ **do**

3: **for** $j \leftarrow 1, ..., (i-1)$ **do**

4: $v_i \leftarrow v_i - \langle v_i, v_j \rangle_{G^*} v_j$.

5: **end for**

6: **end for**

7: $v_i = v_i / \|v_i\|_{G^*}$.

8: **end for**

### 3.3 A posteriori estimation of the primal error

We define the residue of $p^{N_{pr},n+1}$ as the following linear form

$$ \langle r(p^{N_{pr},n+1}), v \rangle = \frac{1}{\Delta t} \langle (M + \Delta t A)p^{N_{pr},n+1} - Mp^{N_{pr},n} - \Delta t b, v \rangle, \forall v \in \mathbb{R}^N \quad (3.5) $$

which we also denote $r(p^{N_{pr},n+1}) = r^{n+1}$ for the sake of simplicity. The residual dual norm is next defined as

$$ \|r^{n+1}\|_{-1} = \sup_{v \in \mathbb{R}^N} \frac{\langle r^{n+1}, v \rangle}{\|v\|_{G^*}} \quad (3.6) $$

and we now propose to evaluate the primal reduction error $e^N$ a posteriori (with a computable estimator $\Delta_{pr}^N$) using a new space-time energy norm $\|\cdot\|_{pr}$ independent from the parameter $\xi$.

**Proposition 3.1** (Energy a posteriori error estimate for the primal problem). Denote $A = \frac{1}{2}(A + A^T) + \frac{1}{2}(A - A^T) := A_{sym} + A_{skew}$ the symmetric and skew-symmetric of matrices $A$. For any $\xi$, given lower bounds

$$ \alpha_{A_{sym}, LB}(\xi) \le \inf_{v \in \mathbb{R}^N} \frac{v^T A_{sym} v}{\|v\|_{G^*}^2} := \alpha_{A_{sym}}(\xi) \quad (3.7) $$

$$ \alpha_{G, LB}(\xi) \le \inf_{v \in \mathbb{R}^N} \frac{v^T (M + \Delta t A_{sym}) v}{\|v\|_{G^*}^2} := \alpha_G(\xi) \quad (3.8) $$

there holds

$$ \|e^N\|_{pr} := \left( \sum_{m=1}^{N} \langle e^m, Me^m \rangle + \Delta t \sum_{m=1}^{N} \langle e^m, A_{sym}^* e^m \rangle \right)^{1/2} \le \Delta_{pr}^N, \quad (3.9) $$
**Algorithm 2 POD with ric**

Input: $S_{pr} \in \mathbb{R}^{N \times N}$, ric $\in (0, 1)$

Output: $G^*$-orthonormal $[z_1|...|z_{\delta N_{pr}}]$.

1: **for** i $\leftarrow 1, ..., N$ **do**

2: $V_i$, $\sigma_i$ i largest vector and eigenvalue pair of $C = S_{pr}^T G^* S_{pr}$

3: **end for**

4: $\delta N_{pr} \leftarrow 1$

$$ \sum^{\delta N_{pr}} \lambda_n $$

5: **while** $E_{\delta N_{pr}} = \frac{\sum_{n=1}^{n=1} \lambda_n}{N} < ric$ **do**

$$ \sum^{\delta N_{pr}} \lambda_n $$

6: $\delta N_{pr} \leftarrow 1 + \delta N_{pr}$

7: $z_{\delta N_{pr}} = \frac{1}{\sqrt{\sigma_{\delta N_{pr}}}} S_{pr} \tilde{V}_{\delta N_{pr}}$

8: **end while**

**Algorithm 3** Gram-Schmidt with re-iteration

Input: vectors $v_i$, $i \in 1, ..., N_{pr}$.

Output: orthonormal vectors $v_i$.

1: **for** $\ell = 1, 2$ **do**

2: **for** $i \leftarrow 1, ..., N_{pr}$ **do**

3: **for** $j \leftarrow 1, ..., (i-1)$ **do**

4: $v_i \leftarrow v_i - \langle v_i, v_j \rangle_{G^*} v_j$.

5: **end for**

6: **end for**

7: $v_i = v_i / \|v_i\|_{G^*}$.

8: **end for**

### 3.3 A posteriori estimation of the primal error

We define the residue of $p^{N_{pr},n+1}$ as the following linear form

$$ \langle r(p^{N_{pr},n+1}), v \rangle = \frac{1}{\Delta t} \langle (M + \Delta t A)p^{N_{pr},n+1} - Mp^{N_{pr},n} - \Delta t b, v \rangle, \forall v \in \mathbb{R}^N \quad (3.5) $$

which we also denote $r(p^{N_{pr},n+1}) = r^{n+1}$ for the sake of simplicity. The residual dual norm is next defined as

$$ \|r^{n+1}\|_{-1} = \sup_{v \in \mathbb{R}^N} \frac{\langle r^{n+1}, v \rangle}{\|v\|_{G^*}} \quad (3.6) $$

and we now propose to evaluate the primal reduction error $e^N$ a posteriori (with a computable estimator $\Delta_{pr}^N$) using a new space-time energy norm $\|\cdot\|_{pr}$ independent from the parameter $\xi$.

**Proposition 3.1** (Energy a posteriori error estimate for the primal problem). Denote $A = \frac{1}{2}(A + A^T) + \frac{1}{2}(A - A^T) := A_{sym} + A_{skew}$ the symmetric and skew-symmetric of matrices $A$. For any $\xi$, given lower bounds

$$ \alpha_{A_{sym}, LB}(\xi) \le \inf_{v \in \mathbb{R}^N} \frac{v^T A_{sym} v}{\|v\|_{G^*}^2} := \alpha_{A_{sym}}(\xi) \quad (3.7) $$

$$ \alpha_{G, LB}(\xi) \le \inf_{v \in \mathbb{R}^N} \frac{v^T (M + \Delta t A_{sym}) v}{\|v\|_{G^*}^2} := \alpha_G(\xi) \quad (3.8) $$

there holds

$$ \|e^N\|_{pr} := \left( \sum_{m=1}^{N} \langle e^m, Me^m \rangle + \Delta t \sum_{m=1}^{N} \langle e^m, A_{sym}^* e^m \rangle \right)^{1/2} \le \Delta_{pr}^N, \quad (3.9) $$
where the upper bound is taken as

$$ \Delta_{\text{pr}}^{N} := \left( \frac{T + \Delta t}{\alpha_{G, \text{LB}} \alpha_{A_{\text{sym}}, \text{LB}}} \sum_{m=1}^{N} \|\boldsymbol{r}^{m}\|_{-1}^{2} \right)^{1/2}, \qquad (3.10) $$

and where $A_{\text{sym}}^*$ is the symmetric part of $A^*$ for a specific parameter $\xi^*$.

*Proof.* See Appendix A.1. □

Note that with Prop.3.1, the primal reduction error can be evaluated numerically using the same norm of $L^2([0, T]; H^1(\Omega))$-type for all values of the parameter $\xi$, as opposed to [21, Prop. 4.3] e.g.

A typical choice for the symmetric and positive definite matrix $G^*$ (in our numerical results of Section 4 e.g.) is

$$ G^* = M + \Delta t A_{\text{sym}}^* $$

which allows for the simplifications $|||e^N|||_{\text{pr}}^2 = \sum_{m=1}^{N} \|e^m\|_{G^*}^2$ and

$$ \alpha_G(\xi) = \Delta t \alpha_{A_{\text{sym}}}(\xi) + \alpha_M \qquad (3.11) $$

where $\alpha_M := \inf_{v \in \mathbb{R}^N} \frac{v^T M v}{\|v\|_{G^*}^2}$ can be computed once for all, independently of the parameter $\xi$.

## 3.4 A posteriori estimation of the QOI error

In our goal-oriented setting, one is mostly interested by the QOI (2.20) for many values of $\xi$

$$ \langle l, p_{M}^{n} \rangle = s^{n}. \qquad (3.12) $$

Then, a reduced basis can in fact be constructed after modifying the POD-Greedy Algorithm 1 based on $\Delta_{\text{pr}}^N$ into a POD-Greedy based on an a posteriori estimator $\Delta_s^N$ of the QOI reduction error, see Algorithm 4.

To define an a posteriori estimator $\Delta_s^N$ for a QOI linear in the primal problem, let us now introduce for all $n = 1, \dots, N$ a dual problem which evolves backward in time

$$ M \psi_{M,n}^n = -l, \qquad (3.13a) $$

$$ (M + \Delta t A^T) \psi_{M,n}^m = M \psi_{M,n}^{m+1} \quad m = 0, \dots, n-1. \qquad (3.13b) $$

Since $M$, $A^T$ and $l$ do not depend on time, we only solve once the following problem:

$$ M \Psi_{M}^{N} = -l, \qquad (3.14a) $$

$$ (M + \Delta t A^T) \Psi_{M}^{n} = M \Psi_{M}^{n+1} \quad n = 0, \dots, N-1. \qquad (3.14b) $$

Then we appropriately shift the results by defining

$$ \psi_{M,n}^m = \Psi_{M}^{N-n+m} \quad m = 0, \dots, n. \qquad (3.15) $$

Let us also introduce a reduced basis $Z_{\text{du}}$ for the dual problem, such that at time $t = n$ one computes an approximation of the dual solution in (3.14) as

$$ \Psi^{N, n} = Z_{\text{du}} \tilde{\Psi}^{n}, \qquad (3.16) $$

by a Galerkin projection, with $\tilde{\Psi}^n$ solution to

$$ (Z_{\text{du}}^T M Z_{\text{du}} + \Delta t Z_{\text{du}}^T A^T Z_{\text{du}}) \tilde{\Psi}^n = Z_{\text{du}}^T M Z_{\text{du}} \tilde{\Psi}^{n+1}. \qquad (3.17) $$

We denote by $\boldsymbol{\varepsilon}^n = \Psi_{M}^n - \Psi^{N, n}$ the dual reduction error at $t = n$ and $\boldsymbol{\varrho}^n$ the residue associated with the dual problem

$$ \langle \boldsymbol{\varrho}^n, \boldsymbol{v} \rangle = \frac{1}{\Delta t} \langle (M + \Delta t A^T) \Psi^{N, n} - M \Psi^{N, n+1}, \boldsymbol{v} \rangle, \qquad (3.18) $$

with a residual dual norm taken as

$$ \| \boldsymbol{\varrho}^n \|_{-1} = \sup_{\boldsymbol{v} \in \mathbb{R}^N} \frac{\langle \boldsymbol{\varrho}^n, \boldsymbol{v} \rangle}{\| \boldsymbol{v} \|_{G^*}}. \qquad (3.19) $$
where the upper bound is taken as

$$ \Delta_{\text{pr}}^{N} := \left( \frac{T + \Delta t}{\alpha_{G, \text{LB}} \alpha_{A_{\text{sym}}, \text{LB}}} \sum_{m=1}^{N} \|\boldsymbol{r}^{m}\|_{-1}^{2} \right)^{1/2}, \qquad (3.10) $$

and where $A_{\text{sym}}^*$ is the symmetric part of $A^*$ for a specific parameter $\xi^*$.

*Proof.* See Appendix A.1. □

Note that with Prop.3.1, the primal reduction error can be evaluated numerically using the same norm of $L^2([0, T]; H^1(\Omega))$-type for all values of the parameter $\xi$, as opposed to [21, Prop. 4.3] e.g.

A typical choice for the symmetric and positive definite matrix $G^*$ (in our numerical results of Section 4 e.g.) is

$$ G^* = M + \Delta t A_{\text{sym}}^* $$

which allows for the simplifications $\|e^N\|_{\text{pr}}^2 = \sum_{m=1}^{N} \|e^m\|_{G^*}^2$ and

$$ \alpha_G(\xi) = \Delta t \alpha_{A_{\text{sym}}}(\xi) + \alpha_M \qquad (3.11) $$

where $\alpha_M := \inf_{v \in \mathbb{R}^N} \frac{v^T M v}{\|v\|_{G^*}^2}$ can be computed once for all, independently of the parameter $\xi$.

## 3.4 A posteriori estimation of the QOI error

In our goal-oriented setting, one is mostly interested by the QOI (2.20) for many values of $\xi$

$$ \langle l, p_{M}^{n} \rangle = s^{n}. \qquad (3.12) $$

Then, a reduced basis can in fact be constructed after modifying the POD-Greedy Algorithm 1 based on $\Delta_{\text{pr}}^N$ into a POD-Greedy based on a posteriori estimator $\Delta_s^N$ of the QOI reduction error, see Algorithm 4.

To define a posteriori estimator $\Delta_s^N$ for a QOI linear in the primal problem, let us now introduce for all $n = 1, \dots, N$ a dual problem which evolves backward in time

$$ M \psi_{M,n}^n = -l, \qquad (3.13a) $$

$$ (M + \Delta t A^T) \psi_{M,n}^m = M \psi_{M,n}^{m+1} \quad m = 0, \dots, n-1. \qquad (3.13b) $$

Since $M$, $A^T$ and $l$ do not depend on time, we only solve once the following problem:

$$ M \Psi_{M}^{N} = -l, \qquad (3.14a) $$

$$ (M + \Delta t A^T) \Psi_{M}^{n} = M \Psi_{M}^{n+1} \quad n = 0, \dots, N-1. \qquad (3.14b) $$

Then we appropriately shift the results by defining

$$ \psi_{M,n}^{m} = \Psi_{M}^{N-n+m} \quad m = 0, \dots, n. \qquad (3.15) $$

Let us also introduce a reduced basis $Z_{\text{du}}$ for the dual problem, such that at time $t = n$ one computes an approximation of the dual solution in (3.14) as

$$ \Psi^{N_{\text{du}}, n} = Z_{\text{du}} \tilde{\Psi}^{n}, \qquad (3.16) $$

by a Galerkin projection, with $\tilde{\Psi}^n$ solution to

$$ (Z_{\text{du}}^T M Z_{\text{du}} + \Delta t Z_{\text{du}}^T A^T Z_{\text{du}}) \tilde{\Psi}^n = Z_{\text{du}}^T M Z_{\text{du}} \tilde{\Psi}^{n+1}. \qquad (3.17) $$

We denote by $\boldsymbol{\varepsilon}^n = \Psi_{M}^n - \Psi^{N_{\text{du}}, n}$ the dual reduction error at $t = n$ and $\boldsymbol{\varrho}^n$ the residue associated with the dual problem

$$ \langle \boldsymbol{\varrho}^n, \boldsymbol{v} \rangle = \frac{1}{\Delta t} \langle (M + \Delta t A^T) \Psi^{N_{\text{du}}, n} - M \Psi^{N_{\text{du}}, n+1}, \boldsymbol{v} \rangle, \qquad (3.18) $$

with a residual dual norm taken as

$$ \| \boldsymbol{\varrho}^n \|_{-1} = \sup_{\boldsymbol{v} \in \mathbb{R}^N} \frac{\langle \boldsymbol{\varrho}^n, \boldsymbol{v} \rangle}{\| \boldsymbol{v} \|_{G^*}}. \qquad (3.19) $$
**Proposition 3.2** (Energy a posteriori error estimate for the dual problem). Given the same data as in Prop. 3.1, there holds for all $\xi$

$$ \|\boldsymbol{\epsilon}^N\|_{\text{du}} := \left( \sum_{m=0}^{N-1} \langle \boldsymbol{\epsilon}^m, M\boldsymbol{\epsilon}^m \rangle + \Delta t \sum_{m=0}^{N-1} \langle \boldsymbol{\epsilon}^m, \boldsymbol{A}_{\text{sym}}^* \boldsymbol{\epsilon}^m \rangle \right)^{1/2} \le \Delta_{\text{du}}^N \quad (3.20a) $$

in the space-time energy norm $\|\cdot\|_{\text{du}}$ independent of the parameter $\xi$, with

$$ \Delta_{\text{du}}^N := \left( \frac{T + \Delta t}{\alpha_{\mathbf{G},\text{LB}} \alpha_{\mathbf{A}_{\text{sym}},\text{LB}}} \sum_{m=0}^{N-1} \|\boldsymbol{\varrho}^m\|_{-1}^2 \right)^{1/2} . \quad (3.20b) $$

*Proof.* See Appendix A.2. □

**Proposition 3.3** (Output error evaluation). Given the same data as in Prop. 3.1, one can define two reduced outputs:

$$ s^{\mathbf{N},n} = \langle \mathbf{l}, \mathbf{p}^{\mathbf{N}_{\text{pr}},n} \rangle + \Delta t \sum_{n'=0}^{n-1} \langle \mathbf{r}^{n'+1}, \Psi^{\mathbf{N}_{\text{du}},N-n+n'} \rangle \quad (3.21) $$

with approximation error bounded as

$$ |s^N - s^{\mathbf{N},N}| \le \Delta t \left( \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2} \Delta_{\text{du}}^N := \Delta_s^N, \quad (3.22) $$

or

$$ \tilde{s}^{\mathbf{N},n} = \langle \mathbf{l}, \mathbf{p}^{\mathbf{N}_{\text{pr}},n} \rangle \quad (3.23) $$

with approximation error bounded as

$$ |s^N - \tilde{s}^{\mathbf{N},N}| \le \Delta t \left( \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2} \Delta_{\text{du}}^N + \Delta t \sum_{n=0}^{N-1} \left| \langle \mathbf{r}^{n+1}, \Psi^{\mathbf{N}_{\text{du}},n} \rangle \right| := \tilde{\Delta}_s^N. \quad (3.24) $$

*Proof.* See Appendix A.3. □

The optimal selection between these two definitions, based on their accuracy and efficiency, will be elucidated in Section 4. There, we present a comparative analysis of the numerical results obtained by POD-Greedy algorithms with $\Delta_s^N$ and $\tilde{\Delta}_s^N$, where the construction of $Z_{\text{du}}$ is simultaneous to that of $Z_{\text{pr}}$, see e.g. Algorithm 4.

## 3.5 Computational aspects

Having addressed the offline computation of the reduced basis, we notice that for online computations, the reduced systems (3.3) and (3.17), required to evaluate (3.22), (3.24), depend on $\mathcal{N}$. An affine decomposition strategy is a classical way to ensure the rapid assembly of the reduced system in the online stage. In the following, we also discuss the practical computation of the residual dual norm and the coercivity constant.

**Affine decomposition.** To construct the reduced matrix $\mathbf{A}^{\mathbf{N}_{\text{pr}}} \in \mathbb{R}^{\mathbf{N}_{\text{pr}} \times \mathbf{N}_{\text{pr}}}$ and reduced vector $\mathbf{b}^{\mathbf{N}_{\text{pr}}} \in \mathbb{R}^{\mathbf{N}_{\text{pr}}}$ defined as $\mathbf{A}^{\mathbf{N}_{\text{pr}}} = \mathbf{Z}_{\text{pr}}^T A \mathbf{Z}_{\text{pr}}$ and $\mathbf{b}^{\mathbf{N}_{\text{pr}}} = \mathbf{Z}_{\text{pr}}^T \mathbf{b}$ in equation (3.3), we still need to compute $\mathbf{A}$ and $\mathbf{b}$ depending on $\xi$. This evaluation, which depends on $\mathcal{N}$, is detrimental to the rapid online evaluation of the reduced basis solution when varying the parameter values. To accelerate the construction, we rewrite $\mathbf{A}$ and $\mathbf{b}$ as

$$ \mathbf{A} = \sum_{d=1}^{\mathbf{D}_\mathbf{a}} \theta_d^a(\xi) \mathbf{A}_d, \quad \mathbf{b} = \sum_{d=1}^{\mathbf{D}_\mathbf{b}} \theta_d^b(\xi) \mathbf{b}_d. \quad (3.25) $$

In (3.25), $\mathbf{A}_d \in \mathbb{R}^{\mathcal{N} \times \mathcal{N}}$, $1 \le d \le \mathbf{D}_\mathbf{a}$ and $\mathbf{b}_d \in \mathbb{R}^{\mathcal{N}}$, $1 \le d \le \mathbf{D}_\mathbf{b}$ do not depend on $\xi$ and they are computed and stored once during the whole offline stage. Thus, for each new parameter $\xi$, we only have to compute the two sets of scalars $\{\theta_d^a(\xi)\}_{d=1}^{\mathbf{D}_\mathbf{a}}$ and $\{\theta_d^b(\xi)\}_{d=1}^{\mathbf{D}_\mathbf{b}}$ and assemble $\mathbf{A}^{\mathbf{N}_{\text{pr}}}$ and $\mathbf{b}^{\mathbf{N}_{\text{pr}}}$. This operation only depends on the dimension $\mathbf{N}_{\text{pr}}$ of the reduced basis. Note that, in our case this affine decomposition does not exist. We therefore use the Empirical Interpolation Method (EIM) [19, 25] (see also Appendix C) to build such an approximation for $\mathbf{A}$ and $\mathbf{b}$. Taking into account the definition of the numerical flux given by (2.8), (2.15) and (2.17), we consider the
**Proposition 3.2** (Energy a posteriori error estimate for the dual problem). Given the same data as in Prop. 3.1, there holds for all $\xi$

$$ \|\boldsymbol{\epsilon}^N\|_{\text{du}} := \left( \sum_{m=0}^{N-1} \langle \boldsymbol{\epsilon}^m, M\boldsymbol{\epsilon}^m \rangle + \Delta t \sum_{m=0}^{N-1} \langle \boldsymbol{\epsilon}^m, \boldsymbol{A}_{\text{sym}}^* \boldsymbol{\epsilon}^m \rangle \right)^{1/2} \le \Delta_{\text{du}}^N \quad (3.20a) $$

in the space-time energy norm $\|\cdot\|_{\text{du}}$ independent of the parameter $\xi$, with

$$ \Delta_{\text{du}}^N := \left( \frac{T + \Delta t}{\alpha_{\mathbf{G},\text{LB}} \alpha_{\mathbf{A}_{\text{sym}},\text{LB}}} \sum_{m=0}^{N-1} \|\boldsymbol{\varrho}^m\|_{-1}^2 \right)^{1/2}. \quad (3.20b) $$

*Proof.* See Appendix A.2. □

**Proposition 3.3** (Output error evaluation). Given the same data as in Prop. 3.1, one can define two reduced outputs:

$$ s^{\mathbf{N},n} = \langle \mathbf{l}, \mathbf{p}^{\mathbf{N}_{\text{pr}},n} \rangle + \Delta t \sum_{n'=0}^{n-1} \langle \mathbf{r}^{n'+1}, \Psi^{\mathbf{N}_{\text{du}},N-n+n'} \rangle \quad (3.21) $$

with approximation error bounded as

$$ |s^N - s^{\mathbf{N},N}| \le \Delta t \left( \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2} \Delta_{\text{du}}^N := \Delta_s^N, \quad (3.22) $$

or

$$ \tilde{s}^{\mathbf{N},n} = \langle \mathbf{l}, \mathbf{p}^{\mathbf{N}_{\text{pr}},n} \rangle \quad (3.23) $$

with approximation error bounded as

$$ |s^N - \tilde{s}^{\mathbf{N},N}| \le \Delta t \left( \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2} \Delta_{\text{du}}^N + \Delta t \sum_{n=0}^{N-1} \left| \langle \mathbf{r}^{n+1}, \Psi^{\mathbf{N}_{\text{du}},n} \rangle \right| := \tilde{\Delta}_s^N. \quad (3.24) $$

*Proof.* See Appendix A.3. □

The optimal selection between these two definitions, based on their accuracy and efficiency, will be elucidated in Section 4. There, we present a comparative analysis of the numerical results obtained by POD-Greedy algorithms with $\Delta_s^N$ and $\tilde{\Delta}_s^N$, where the construction of $Z_{\text{du}}$ is simultaneous to that of $Z_{\text{pr}}$, see e.g. Algorithm 4.

## 3.5 Computational aspects

Having addressed the offline computation of the reduced basis, we notice that for online computations, the reduced systems (3.3) and (3.17), required to evaluate (3.22), (3.24), depend on $\mathcal{N}$. An affine decomposition strategy is a classical way to ensure the rapid assembly of the reduced system in the online stage. In the following, we also discuss the practical computation of the residual dual norm and the coercivity constant.

**Affine decomposition.** To construct the reduced matrix $\mathbf{A}^{\mathbf{N}_{\text{pr}}} \in \mathbb{R}^{\mathbf{N}_{\text{pr}} \times \mathbf{N}_{\text{pr}}}$ and reduced vector $\mathbf{b}^{\mathbf{N}_{\text{pr}}} \in \mathbb{R}^{\mathbf{N}_{\text{pr}}}$ defined as $\mathbf{A}^{\mathbf{N}_{\text{pr}}} = \mathbf{Z}_{\text{pr}}^T A \mathbf{Z}_{\text{pr}}$ and $\mathbf{b}^{\mathbf{N}_{\text{pr}}} = \mathbf{Z}_{\text{pr}}^T \mathbf{b}$ in equation (3.3), we still need to compute $\mathbf{A}$ and $\mathbf{b}$ depending on $\xi$. This evaluation, which depends on $\mathcal{N}$, is detrimental to the rapid online evaluation of the reduced basis solution when varying the parameter values. To accelerate the construction, we rewrite $\mathbf{A}$ and $\mathbf{b}$ as

$$ \mathbf{A} = \sum_{d=1}^{\mathbf{D}_\mathbf{a}} \theta_d^a(\xi) \mathbf{A}_d, \quad \mathbf{b} = \sum_{d=1}^{\mathbf{D}_\mathbf{b}} \theta_d^b(\xi) \mathbf{b}_d. \quad (3.25) $$

In (3.25), $\mathbf{A}_d \in \mathbb{R}^{\mathcal{N} \times \mathcal{N}}$, $1 \le d \le \mathbf{D}_\mathbf{a}$ and $\mathbf{b}_d \in \mathbb{R}^{\mathcal{N}}$, $1 \le d \le \mathbf{D}_\mathbf{b}$ do not depend on $\xi$ and they are computed and stored once during the whole offline stage. Thus, for each new parameter $\xi$, we only have to compute the two sets of scalars $\{\theta_d^a(\xi)\}_{d=1}^{\mathbf{D}_\mathbf{a}}$ and $\{\theta_d^b(\xi)\}_{d=1}^{\mathbf{D}_\mathbf{b}}$ and assemble $\mathbf{A}^{\mathbf{N}_{\text{pr}}}$ and $\mathbf{b}^{\mathbf{N}_{\text{pr}}}$. This operation only depends on the dimension $\mathbf{N}_{\text{pr}}$ of the reduced basis. Note that, in our case this affine decomposition does not exist. We therefore use the Empirical Interpolation Method (EIM) [19, 25] (see also Appendix C) to build such an approximation for $\mathbf{A}$ and $\mathbf{b}$. Taking into account the definition of the numerical flux given by (2.8), (2.15) and (2.17), we consider the
**Algorithm 4 POD-Greedy Algorithm with $\Delta_s^N$**

1: Procedure $[Z_{pr}, Z_{du}] = \text{POD-Greedy}(\mathbf{N}_{\text{max}}, \epsilon_{\text{tol}}, \Xi, \text{ric})$.

2: $\mathbf{N}_{pr} = 0, \mathbf{N}_{du} = 0$.

3: $\delta^{\mathbf{N}_s} = \epsilon_{\text{tol}} + 1$.

4: Take $\xi_1 \in \Xi$, $\ell = 1$ and set $\Xi^\ell = \{\xi_1\}$.

5: Define $Z_{pr} = \emptyset$ and $Z_{du} = \emptyset$.

6: **while** $\delta^{\mathbf{N}_s} > \epsilon_{\text{tol}}$ and $\mathbf{N}_{pr} < \mathbf{N}_{\text{max}}$ **do**.

7: Compute $\mathbf{p}_{\mathcal{M}}^n(\xi_\ell)$ for $1 \le n \le N$.

8: Compute $\mathbf{\Psi}_{\mathcal{M}}^n(\xi_\ell)$ for $0 \le n \le N - 1$.

9: Set $\mathbf{S}_{pr} := [\mathbf{p}_{\mathcal{M}}^1(\xi_\ell) - \text{Proj}_{Z_{pr}}(\mathbf{p}_{\mathcal{M}}^1(\xi_\ell)) | \dots | \mathbf{p}_{\mathcal{M}}^N(\xi_\ell) - \text{Proj}_{Z_{pr}}(\mathbf{p}_{\mathcal{M}}^N(\xi_\ell))]$.

10: Set $\mathbf{S}_{du} := [\mathbf{\Psi}_{\mathcal{M}}^0(\xi_\ell) - \text{Proj}_{Z_{du}}(\mathbf{\Psi}_{\mathcal{M}}^0(\xi_\ell)) | \dots | \mathbf{\Psi}_{\mathcal{M}}^{N-1}(\xi_\ell) - \text{Proj}_{Z_{du}}(\mathbf{\Psi}_{\mathcal{M}}^{N-1}(\xi_\ell))]$.

11: Compute $\tilde{Z}_{pr}^{\delta \mathbf{N}_{pr}} = \text{POD}(\mathbf{S}_{pr}, \text{ric})$.

12: Compute $\tilde{Z}_{du}^{\delta \mathbf{N}_{du}} = \text{POD}(\mathbf{S}_{du}, \text{ric})$.

13: Define $\mathbf{Z}_{pr}^{\mathbf{N}_{pr}+\delta \mathbf{N}_{pr}} := \text{orthonormalize}(\mathbf{Z}_{pr}^{\mathbf{N}_{pr}} \cup \{\tilde{Z}_{pr}^{\delta \mathbf{N}_{pr}}\})$ using Algorithm 3.

14: Define $\mathbf{Z}_{du}^{\mathbf{N}_{du}+\delta \mathbf{N}_{du}} := \text{orthonormalize}(\mathbf{Z}_{du}^{\mathbf{N}_{du}} \cup \{\tilde{Z}_{du}^{\delta \mathbf{N}_{du}}\})$ using Algorithm 3.

15: Compute $\delta^{\mathbf{N}_s} = \max_{\xi \in \Xi} \Delta_s^N$.

16: Set $\xi_{\ell+1} = \arg\max_{\xi \in \Xi} \Delta_s^N$.

17: $\Xi^{\ell+1} \leftarrow \Xi^\ell \cup \{\xi_{\ell+1}\}$.

18: $\mathbf{N}_{pr} \leftarrow \mathbf{N}_{pr} + \delta \mathbf{N}_{pr}$.

19: $\mathbf{N}_{du} \leftarrow \mathbf{N}_{du} + \delta \mathbf{N}_{du}$.

20: $\ell \leftarrow \ell + 1$.

21: **end while**

vector $\hat{\mathbf{v}} = ((\alpha_{K,\sigma\sigma'})_{K \in \mathcal{T}, \sigma \in \mathcal{E}_K, \sigma' \in \mathcal{S}_{K,\sigma}}, (\alpha_{K,\sigma\sigma'}\omega_{M,\sigma'})_{K \in \mathcal{T}, M \in \mathcal{T}_{\sigma'}, \sigma \in \mathcal{E}_K, \sigma' \in \mathcal{S}_{K,\sigma}, \sigma' \in \mathcal{E}_{\text{int}}})$ and seek for a linearization of it depending on the parameter $\xi \in \Xi$ through the operator $\mathcal{I}_{\mathbf{M}_{\text{EIM}}}$ such that

$$ \hat{\mathbf{v}}(\xi) \approx \sum_{d=1}^{\mathbf{M}_{\text{EIM}}} \theta_d(\xi) \tilde{\mathbf{v}}^d := \mathcal{I}_{\mathbf{M}_{\text{EIM}}} [\hat{\mathbf{v}}(\xi)], $$

where $\theta_d(\xi) \in \mathbb{R}$. If such an approximation exists, we can then replace the terms of each vector $\tilde{\mathbf{v}}^d$, $1 \le d \le \mathbf{M}_{\text{EIM}}$ in the flux formula and obtain the matrices $\mathbf{A}_d$ and the vectors $\mathbf{b}_d$, $1 \le d \le \mathbf{M}_{\text{EIM}}$ independently from the parameter $\xi$. In terms of online cost, we need $\mathcal{O}((\mathbf{D}_a + 1)\mathbf{N}_{pr}^2)$ and $\mathcal{O}(\mathbf{D}_b\mathbf{N}_{pr})$ to assemble the left-hand side and right-hand side respectively in (3.3). The reduced system is then solved with $\mathcal{O}(\mathbf{N}\mathbf{N}_{pr}^3)$.

**Residual norm evaluation.** Using the affine decomposition and the fact that $\|\mathbf{r}^{n+1}\|_{-1}^2 = (\mathbf{r}^{n+1})^T(\mathbf{G}^*)^{-1}\mathbf{r}^{n+1}$ (which results from Cauchy-Schwartz inequality), we can now rewrite the dual norm of the residual as

$$ \begin{aligned} \|\mathbf{r}^{n+1}\|_{-1}^2 = & \sum_{d=1}^{\mathbf{D}_b} \sum_{d'=1}^{\mathbf{D}_b} \theta_d^b(\xi) \theta_{d'}^b(\xi) \mathbf{b}_d^T(\mathbf{G}^*)^{-1} \mathbf{b}_{d'} - 2 \sum_{d=1}^{\mathbf{D}_a} \sum_{d'=1}^{\mathbf{D}_b} \theta_d^a(\xi) \theta_{d'}^b(\xi) (\tilde{\mathbf{p}}^{n+1})^T \mathbf{Z}_{pr}^T \mathbf{A}_d^T(\mathbf{G}^*)^{-1} \mathbf{b}_{d'} \\ & + \sum_{d=1}^{\mathbf{D}_a} \sum_{d'=1}^{\mathbf{D}_a} \theta_d^a(\xi) \theta_{d'}^a(\xi) (\tilde{\mathbf{p}}^{n+1})^T \mathbf{Z}_{pr}^T \mathbf{A}_d^T(\mathbf{G}^*)^{-1} \mathbf{A}_{d'} \mathbf{Z}_{pr} \tilde{\mathbf{p}}^{n+1} \\ & - \frac{2}{\Delta t} \sum_{d=1}^{\mathbf{D}_b} \theta_d^b(\xi) (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n)^T \mathbf{Z}_{pr}^T \mathbf{M}(\mathbf{G}^*)^{-1} \mathbf{b}_d + \frac{2}{\Delta t} \sum_{d=1}^{\mathbf{D}_a} \theta_d^a(\xi) (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n)^T \mathbf{Z}_{pr}^T \mathbf{M}(\mathbf{G}^*)^{-1} \mathbf{A}_d \mathbf{Z}_{pr} \tilde{\mathbf{p}}^{n+1} \\ & + \frac{1}{\Delta t^2} (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n)^T \mathbf{Z}_{pr}^T \mathbf{M}(\mathbf{G}^*)^{-1} \mathbf{M} \mathbf{Z}_{pr} (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n). \end{aligned} \quad (3.26) $$

Its evaluation is very sensitive to round-off errors as stressed in [9]. Hence, a naive implementation of (3.26) may suffer from accuracy issues. These can be circumvented by following the method of [7]. We first introduce the Riesz's representative of the residual $\mathcal{R}$ such that

$$ \langle \mathcal{R}(\mathbf{r}^{n+1}), \mathbf{v} \rangle_{\mathbf{G}^*} = \mathbf{r}^{n+1}(\mathbf{v}), \quad \forall \mathbf{v} \in \mathbb{R}^N. $$
**Algorithm 4 POD-Greedy Algorithm with $\Delta_s^N$**

1: **Procedure** $[Z_{pr}, Z_{du}] = \text{POD-Greedy}(\mathbf{N}_{\text{max}}, \epsilon_{\text{tol}}, \Xi, \text{ric})$.

2: $\mathbf{N}_{pr} = 0, \mathbf{N}_{du} = 0$.

3: $\delta^{\mathbf{N}_s} = \epsilon_{\text{tol}} + 1$.

4: Take $\xi_1 \in \Xi$, $\ell = 1$ and set $\Xi^\ell = \{\xi_1\}$.

5: Define $Z_{pr} = \emptyset$ and $Z_{du} = \emptyset$.

6: **while** $\delta^{\mathbf{N}_s} > \epsilon_{\text{tol}}$ and $\mathbf{N}_{pr} < \mathbf{N}_{\text{max}}$ **do**.

7: Compute $\mathbf{p}_{\mathcal{M}}^n(\xi_\ell)$ for $1 \le n \le N$.

8: Compute $\mathbf{\Psi}_{\mathcal{M}}^n(\xi_\ell)$ for $0 \le n \le N - 1$.

9: Set $\mathbf{S}_{pr} := [\mathbf{p}_{\mathcal{M}}^1(\xi_\ell) - \text{Proj}_{Z_{pr}}(\mathbf{p}_{\mathcal{M}}^1(\xi_\ell)) | \dots | \mathbf{p}_{\mathcal{M}}^N(\xi_\ell) - \text{Proj}_{Z_{pr}}(\mathbf{p}_{\mathcal{M}}^N(\xi_\ell))]$.

10: Set $\mathbf{S}_{du} := [\mathbf{\Psi}_{\mathcal{M}}^0(\xi_\ell) - \text{Proj}_{Z_{du}}(\mathbf{\Psi}_{\mathcal{M}}^0(\xi_\ell)) | \dots | \mathbf{\Psi}_{\mathcal{M}}^{N-1}(\xi_\ell) - \text{Proj}_{Z_{du}}(\mathbf{\Psi}_{\mathcal{M}}^{N-1}(\xi_\ell))]$.

11: Compute $\tilde{Z}_{pr}^{\delta\mathbf{N}_{pr}} = \text{POD}(\mathbf{S}_{pr}, ric)$.

12: Compute $\tilde{Z}_{du}^{\delta\mathbf{N}_{du}} = \text{POD}(\mathbf{S}_{du}, ric)$.

13: Define $\mathbf{Z}_{pr}^{\mathbf{N}_{pr}+\delta\mathbf{N}_{pr}} := \text{orthonormalize}(\mathbf{Z}_{pr}^{\mathbf{N}_{pr}} \cup \{\tilde{Z}_{pr}^{\delta\mathbf{N}_{pr}}\})$ using Algorithm 3.

14: Define $\mathbf{Z}_{du}^{\mathbf{N}_{du}+\delta\mathbf{N}_{du}} := \text{orthonormalize}(\mathbf{Z}_{du}^{\mathbf{N}_{du}} \cup \{\tilde{Z}_{du}^{\delta\mathbf{N}_{du}}\})$ using Algorithm 3.

15: Compute $\delta^{\mathbf{N}_s} = \max_{\xi \in \Xi} \Delta_s^N$.

16: Set $\xi_{\ell+1} = \arg\max_{\xi \in \Xi} \Delta_s^N$.

17: $\Xi^{\ell+1} \leftarrow \Xi^\ell \cup \{\xi_{\ell+1}\}$.

18: $\mathbf{N}_{pr} \leftarrow \mathbf{N}_{pr} + \delta\mathbf{N}_{pr}$.

19: $\mathbf{N}_{du} \leftarrow \mathbf{N}_{du} + \delta\mathbf{N}_{du}$.

20: $\ell \leftarrow \ell + 1$.

21: **end while**

vector $\hat{\mathbf{v}} = ((\alpha_{K,\sigma\sigma'})_{K \in \mathcal{T}, \sigma \in \mathcal{E}_K, \sigma' \in \mathcal{S}_{K,\sigma}}, (\alpha_{K,\sigma\sigma'}\omega_{M,\sigma'})_{K \in \mathcal{T}, M \in \mathcal{T}_{\sigma'}, \sigma \in \mathcal{E}_K, \sigma' \in \mathcal{S}_{K,\sigma}, \sigma' \in \mathcal{E}_{\text{int}}})$ and seek for a linearization of it depending on the parameter $\xi \in \Xi$ through the operator $\mathcal{I}_{\mathbf{M}_{\text{EIM}}}$ such that

$$ \hat{\mathbf{v}}(\xi) \approx \sum_{d=1}^{\mathbf{M}_{\text{EIM}}} \theta_d(\xi) \tilde{\mathbf{v}}^d := \mathcal{I}_{\mathbf{M}_{\text{EIM}}} [\hat{\mathbf{v}}(\xi)], $$

where $\theta_d(\xi) \in \mathbb{R}$. If such an approximation exists, we can then replace the terms of each vector $\tilde{\mathbf{v}}^d$, $1 \le d \le \mathbf{M}_{\text{EIM}}$ in the flux formula and obtain the matrices $\mathbf{A}_d$ and the vectors $\mathbf{b}_d$, $1 \le d \le \mathbf{M}_{\text{EIM}}$ independently from the parameter $\xi$. In terms of online cost, we need $\mathcal{O}((\mathbf{D}_a + 1)\mathbf{N}_{pr}^2)$ and $\mathcal{O}(\mathbf{D}_b\mathbf{N}_{pr})$ to assemble the left-hand side and right-hand side respectively in (3.3). The reduced system is then solved with $\mathcal{O}(\mathbf{N}\mathbf{N}_{pr}^3)$.

**Residual norm evaluation.** Using the affine decomposition and the fact that $\|\mathbf{r}^{n+1}\|_{-1}^2 = (\mathbf{r}^{n+1})^T(\mathbf{G}^*)^{-1}\mathbf{r}^{n+1}$ (which results from Cauchy-Schwartz inequality), we can now rewrite the dual norm of the residual as

$$ \begin{aligned} \|\mathbf{r}^{n+1}\|_{-1}^2 = & \sum_{d=1}^{\mathbf{D}_b} \sum_{d'=1}^{\mathbf{D}_b} \theta_d^b(\xi) \theta_{d'}^b(\xi) \mathbf{b}_d^T(\mathbf{G}^*)^{-1} \mathbf{b}_{d'} - 2 \sum_{d=1}^{\mathbf{D}_a} \sum_{d'=1}^{\mathbf{D}_b} \theta_d^a(\xi) \theta_{d'}^b(\xi) (\tilde{\mathbf{p}}^{n+1})^T \mathbf{Z}_{pr}^T \mathbf{A}_d^T(\mathbf{G}^*)^{-1} \mathbf{b}_{d'} \\ & + \sum_{d=1}^{\mathbf{D}_a} \sum_{d'=1}^{\mathbf{D}_a} \theta_d^a(\xi) \theta_{d'}^a(\xi) (\tilde{\mathbf{p}}^{n+1})^T \mathbf{Z}_{pr}^T \mathbf{A}_d^T(\mathbf{G}^*)^{-1} \mathbf{A}_{d'} \mathbf{Z}_{pr} \tilde{\mathbf{p}}^{n+1} \\ & - \frac{2}{\Delta t} \sum_{d=1}^{\mathbf{D}_b} \theta_d^b(\xi) (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n)^T \mathbf{Z}_{pr}^T \mathbf{M}(\mathbf{G}^*)^{-1} \mathbf{b}_d + \frac{2}{\Delta t} \sum_{d=1}^{\mathbf{D}_a} \theta_d^a(\xi) (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n)^T \mathbf{Z}_{pr}^T \mathbf{M}(\mathbf{G}^*)^{-1} \mathbf{A}_d \mathbf{Z}_{pr} \tilde{\mathbf{p}}^{n+1} \\ & + \frac{1}{\Delta t^2} (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n)^T \mathbf{Z}_{pr}^T \mathbf{M}(\mathbf{G}^*)^{-1} \mathbf{M} \mathbf{Z}_{pr} (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n). \end{aligned} \quad (3.26) $$

Its evaluation is very sensitive to round-off errors as stressed in [9]. Hence, a naive implementation of (3.26) may suffer from accuracy issues. These can be circumvented by following the method of [7]. We first introduce the Riesz's representative of the residual $\mathcal{R}$ such that

$$ \langle \mathcal{R}(\mathbf{r}^{n+1}), \mathbf{v} \rangle_{\mathbf{G}^*} = \mathbf{r}^{n+1}(\mathbf{v}), \quad \forall \mathbf{v} \in \mathbb{R}^{\mathcal{N}}. $$
Using the affine decomposition (3.25), the Riesz representation of the primal residual is given by

$$ \mathcal{R}(\mathbf{r}^{n+1}) = \frac{1}{\Delta t} (\mathbf{G}^*)^{-1} \mathbf{M} \mathbf{Z}_{\text{pr}} (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n) + \sum_{d=1}^{\mathcal{D}_a} \theta_d^a(\xi) (\mathbf{G}^*)^{-1} \mathbf{A}_d \mathbf{Z}_{\text{pr}} \tilde{\mathbf{p}}^{n+1} - \sum_{d=1}^{\mathcal{D}_b} \theta_d^b(\xi) (\mathbf{G}^*)^{-1} \mathbf{b}_d. \quad (3.27) $$

Setting $\mathcal{D}_r = \mathcal{D}_b + \mathcal{D}_a \mathbf{N}_{\text{pr}} + \mathbf{N}_{\text{pr}}$, we define the coefficient vector $\hat{\mathbf{r}}^{n+1} \in \mathbb{R}^{\mathcal{D}_r}$ as

$$ \hat{\mathbf{r}}^{n+1} = \left( \frac{1}{\Delta t} (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n)^T, \theta_1^a(\xi)(\tilde{\mathbf{p}}^{n+1})^T, \dots, \theta_{\mathcal{D}_a}^a(\xi)(\tilde{\mathbf{p}}^{n+1})^T, -\theta_1^b(\xi), \dots, -\theta_{\mathcal{D}_b}^b(\xi) \right)^T, $$

and the vector $\hat{\boldsymbol{\eta}} \in \mathbb{R}^{\mathcal{D}_r}$ as

$$ \hat{\boldsymbol{\eta}} = ((\mathbf{G}^*)^{-1} \mathbf{M} \mathbf{Z}_{\text{pr}}, (\mathbf{G}^*)^{-1} \mathbf{A}_1 \mathbf{Z}_{\text{pr}}, \dots, (\mathbf{G}^*)^{-1} \mathbf{A}_{\mathcal{D}_a} \mathbf{Z}_{\text{pr}}, (\mathbf{G}^*)^{-1} \mathbf{b}_1, \dots, (\mathbf{G}^*)^{-1} \mathbf{b}_{\mathcal{D}_b})^T. $$

The Riesz representative is then written as

$$ \mathcal{R}(\mathbf{r}^{n+1}) = \sum_{d=1}^{\mathcal{D}_r} \hat{\mathbf{r}}_d^{n+1} \hat{\boldsymbol{\eta}}_d, \quad (3.28) $$

and the norm is given by

$$ \| \mathcal{R}(\mathbf{r}^{n+1}) \|_{\mathbf{G}^*}^2 = \left\langle \sum_{d=1}^{\mathcal{D}_r} \hat{\mathbf{r}}_d^{n+1} \hat{\boldsymbol{\eta}}_d, \sum_{d=1}^{\mathcal{D}_r} \hat{\mathbf{r}}_d^{n+1} \hat{\boldsymbol{\eta}}_d \right\rangle_{\mathbf{G}^*}. \quad (3.29) $$

The evaluation of (3.29) is divided into three steps:

1. We construct an orthonormal basis $\boldsymbol{\zeta}$ of $\hat{\boldsymbol{\eta}}$ by applying a modified Gram-Schmidt algorithm with reorthogonalization (see Algorithm 3).

2. We evaluate each term $\hat{\boldsymbol{\eta}}_d = \sum_{i=1}^{\mathcal{D}_r} \bar{\boldsymbol{\eta}}_{d,i} \boldsymbol{\zeta}_i$, with $\bar{\boldsymbol{\eta}}_{d,i} = \langle \hat{\boldsymbol{\eta}}_d, \boldsymbol{\zeta}_i \rangle_{\mathbf{G}^*}$.

3. We compute (3.29) using

$$ \begin{align} \| \mathcal{R}(\mathbf{r}^{n+1}) \|_{\mathbf{G}^*}^2 &= \left\langle \sum_{d=1}^{\mathcal{D}_r} \hat{\mathbf{r}}_d^{n+1} \left( \sum_{i=1}^{\mathcal{D}_r} \bar{\boldsymbol{\eta}}_{d,i} \boldsymbol{\zeta}_i \right), \sum_{d=1}^{\mathcal{D}_r} \hat{\mathbf{r}}_d^{n+1} \left( \sum_{i=1}^{\mathcal{D}_r} \bar{\boldsymbol{\eta}}_{d,i} \boldsymbol{\zeta}_i \right) \right\rangle_{\mathbf{G}^*} \nonumber \\ &= \sum_{i=1}^{\mathcal{D}_r} \sum_{j=1}^{\mathcal{D}_r} \left( \sum_{d=1}^{\mathcal{D}_r} \hat{\mathbf{r}}_d^{n+1} \bar{\boldsymbol{\eta}}_{d,i} \right) \left( \sum_{d=1}^{\mathcal{D}_r} \hat{\mathbf{r}}_d^{n+1} \bar{\boldsymbol{\eta}}_{d,j} \right) \langle \boldsymbol{\zeta}_i, \boldsymbol{\zeta}_j \rangle_{\mathbf{G}^*} \tag{3.30} \\ &= \sum_{i=1}^{\mathcal{D}_r} \left( \sum_{d=1}^{\mathcal{D}_r} \hat{\mathbf{r}}_d^{n+1} \bar{\boldsymbol{\eta}}_{d,i} \right)^2. \nonumber \end{align} $$

Steps 1 and 2 are performed in the offline stage, while steps 3 is completed during the online stage. We follow the same strategy to define to residual dual norm of the dual problem.

**Coercivity constant computation.** The coercivity constant defined by (3.7) is the minimum of the generalized Rayleigh quotient and we have that $\alpha_{\mathbf{A}_{\text{sym}}}$ is the smallest eigenvalue of the following generalized eigenvalue problem

$$ \mathbf{A}_{\text{sym}} \mathbf{v} = \lambda \mathbf{G}^* \mathbf{v}. \quad (3.31) $$

To avoid the resolution of the generalized eigenvalue problem (3.31) which requires for instance $\mathcal{O}(\mathcal{N}^3)$ using a QR algorithm, we consider the successive constraint method (SCM) [12, 22] (see also Appendix B) which, using (3.25), provides an upper bound $\alpha_{\mathbf{A}_{\text{sym}},\text{UB}}(\xi) \in \mathbb{R}$ and a lower bound $\alpha_{\mathbf{A}_{\text{sym}},\text{LB}}(\xi) \in \mathbb{R}$ for the coercivity constant such that

$$ \alpha_{\mathbf{A}_{\text{sym}},\text{LB}}(\xi) \le \alpha_{\mathbf{A}_{\text{sym}}}(\xi) \le \alpha_{\mathbf{A}_{\text{sym}},\text{UB}}(\xi). $$

The evaluation of these bounds do not depend on $\mathcal{N}$. The coercivity constant $\alpha_{\mathbf{A}_{\text{sym}}}$ in the a posteriori estimation formula is then replaced by its corresponding lower bounds. Noting that, once $\alpha_{\mathbf{A}_{\text{sym}},\text{LB}}(\xi)$ is computed, we can replace it in (3.11) to obtain a lower bound for $\alpha_{\mathbf{G}}$.
Using the affine decomposition (3.25), the Riesz representation of the primal residual is given by

$$ \mathcal{R}(\mathbf{r}^{n+1}) = \frac{1}{\Delta t} (\mathbf{G}^*)^{-1} \mathbf{M} \mathbf{Z}_{\text{pr}} (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n) + \sum_{d=1}^{\text{D}_a} \theta_d^a(\xi) (\mathbf{G}^*)^{-1} \mathbf{A}_d \mathbf{Z}_{\text{pr}} \tilde{\mathbf{p}}^{n+1} - \sum_{d=1}^{\text{D}_b} \theta_d^b(\xi) (\mathbf{G}^*)^{-1} \mathbf{b}_d. \quad (3.27) $$

Setting $D_r = D_b + D_a N_{\text{pr}} + N_{\text{pr}}$, we define the coefficient vector $\hat{\mathbf{r}}^{n+1} \in \mathbb{R}^{\text{D}_r}$ as

$$ \hat{\mathbf{r}}^{n+1} = \left( \frac{1}{\Delta t} (\tilde{\mathbf{p}}^{n+1} - \tilde{\mathbf{p}}^n)^T, \theta_1^a(\xi)(\tilde{\mathbf{p}}^{n+1})^T, \dots, \theta_{\text{D}_a}^a(\xi)(\tilde{\mathbf{p}}^{n+1})^T, -\theta_1^b(\xi), \dots, -\theta_{\text{D}_b}^b(\xi) \right)^T, $$

and the vector $\hat{\boldsymbol{\eta}} \in \mathbb{R}^{\text{D}_r}$ as

$$ \hat{\boldsymbol{\eta}} = ((\mathbf{G}^*)^{-1} \mathbf{M} \mathbf{Z}_{\text{pr}}, (\mathbf{G}^*)^{-1} \mathbf{A}_1 \mathbf{Z}_{\text{pr}}, \dots, (\mathbf{G}^*)^{-1} \mathbf{A}_{\text{D}_a} \mathbf{Z}_{\text{pr}}, (\mathbf{G}^*)^{-1} \mathbf{b}_1, \dots, (\mathbf{G}^*)^{-1} \mathbf{b}_{\text{D}_b})^T. $$

The Riesz representative is then written as

$$ \mathcal{R}(\mathbf{r}^{n+1}) = \sum_{d=1}^{\text{D}_r} \hat{\mathbf{r}}_d^{n+1} \hat{\boldsymbol{\eta}}_d, \quad (3.28) $$

and the norm is given by

$$ \| \mathcal{R}(\mathbf{r}^{n+1}) \|_{\mathbf{G}^*}^2 = \left\langle \sum_{d=1}^{\text{D}_r} \hat{\mathbf{r}}_d^{n+1} \hat{\boldsymbol{\eta}}_d, \sum_{d=1}^{\text{D}_r} \hat{\mathbf{r}}_d^{n+1} \hat{\boldsymbol{\eta}}_d \right\rangle_{\mathbf{G}^*} \quad (3.29) $$

The evaluation of (3.29) is divided into three steps:

1. We construct an orthonormal basis $\boldsymbol{\zeta}$ of $\hat{\boldsymbol{\eta}}$ by applying a modified Gram-Schmidt algorithm with reorthogonalization (see Algorithm 3).

2. We evaluate each term $\hat{\boldsymbol{\eta}}_d = \sum_{i=1}^{\text{D}_r} \bar{\boldsymbol{\eta}}_{d,i} \boldsymbol{\zeta}_i$, with $\bar{\boldsymbol{\eta}}_{d,i} = \langle \hat{\boldsymbol{\eta}}_d, \boldsymbol{\zeta}_i \rangle_{\mathbf{G}^*}$.

3. We compute (3.29) using

$$ \begin{aligned} \| \mathcal{R}(\mathbf{r}^{n+1}) \|_{\mathbf{G}^*}^2 &= \left\langle \sum_{d=1}^{\text{D}_r} \hat{\mathbf{r}}_d^{n+1} \left( \sum_{i=1}^{\text{D}_r} \bar{\boldsymbol{\eta}}_{d,i} \boldsymbol{\zeta}_i \right), \sum_{d=1}^{\text{D}_r} \hat{\mathbf{r}}_d^{n+1} \left( \sum_{i=1}^{\text{D}_r} \bar{\boldsymbol{\eta}}_{d,i} \boldsymbol{\zeta}_i \right) \right\rangle_{\mathbf{G}^*} \\ &= \sum_{i=1}^{\text{D}_r} \sum_{j=1}^{\text{D}_r} \left( \sum_{d=1}^{\text{D}_r} \hat{\mathbf{r}}_d^{n+1} \bar{\boldsymbol{\eta}}_{d,i} \right) \left( \sum_{d=1}^{\text{D}_r} \hat{\mathbf{r}}_d^{n+1} \bar{\boldsymbol{\eta}}_{d,j} \right) \langle \boldsymbol{\zeta}_i, \boldsymbol{\zeta}_j \rangle_{\mathbf{G}^*} \quad (3.30) \\ &= \sum_{i=1}^{\text{D}_r} \left( \sum_{d=1}^{\text{D}_r} \hat{\mathbf{r}}_d^{n+1} \bar{\boldsymbol{\eta}}_{d,i} \right)^2. \end{aligned} $$

Steps 1 and 2 are performed in the offline stage, while steps 3 is completed during the online stage. We follow the same strategy to define to residual dual norm of the dual problem.

**Coercivity constant computation.** The coercivity constant defined by (3.7) is the minimum of the generalized Rayleigh quotient and we have that $\alpha_{\mathbf{A}_{\text{sym}}}$ is the smallest eigenvalue of the following generalized eigenvalue problem

$$ \mathbf{A}_{\text{sym}} \mathbf{v} = \lambda \mathbf{G}^* \mathbf{v}. \quad (3.31) $$

To avoid the resolution of the generalized eigenvalue problem (3.31) which requires for instance $\mathcal{O}(N^3)$ using a QR algorithm, we consider the successive constraint method (SCM) [12, 22] (see also Appendix B) which, using (3.25), provides an upper bound $\alpha_{\mathbf{A}_{\text{sym}},\text{UB}}(\xi) \in \mathbb{R}$ and a lower bound $\alpha_{\mathbf{A}_{\text{sym}},\text{LB}}(\xi) \in \mathbb{R}$ for the coercivity constant such that

$$ \alpha_{\mathbf{A}_{\text{sym}},\text{LB}}(\xi) \le \alpha_{\mathbf{A}_{\text{sym}}}(\xi) \le \alpha_{\mathbf{A}_{\text{sym}},\text{UB}}(\xi). $$

The evaluation of these bounds do not depend on $\mathcal{N}$. The coercivity constant $\alpha_{\mathbf{A}_{\text{sym}}}$ in the a posteriori estimation formula is then replaced by its corresponding lower bounds. Noting that, once $\alpha_{\mathbf{A}_{\text{sym}},\text{LB}}(\xi)$ is computed, we can replace it in (3.11) to obtain a lower bound for $\alpha_{\mathbf{G}}$.
# 4 Numerical results

In this section, we numerically validate the theoretical results obtained for the reduction of problem (2.1). Our main goals are to study both efficiency and computation cost of the proposed estimators.

The following parameters are considered:

$$ 
\begin{array}{l}
\mu = 1.5 \times 10^{-5} \text{ Pa.s,} \quad c_t = 1.4 \times 10^{-7} \text{ Pa}^{-1}, \quad g = 9.81 \text{ m/s}^2, \quad \rho = 700 \text{ kg/m}^3, \\
p_{bh} = 4.13 \times 10^7 \text{ Pa,} \quad \phi = 0.2, \quad r_w = 0.1, \quad z_{bh} = 0 \text{ m.}
\end{array}
$$

We use $p_D = 10^5$ Pa as Dirichlet boundary condition. The total duration of the simulation is $T = 200$ days and the time step is $\Delta t = 10$ days. The initial pressure is defined by

$$ p_K^0 = p_D - \rho g (z_K - z_D), $$

where $z_D = 80$ m.

We consider a three-dimensional domain

$$ \Omega = [-2.686 \cdot 10^{-3} \text{ m}, 1996 \text{ m}] \times [6.1 \cdot 10^{-5} \text{ m}, 1996 \text{ m}] \times [-1000.13 \text{ m}, 2.686 \cdot 10^{-3} \text{ m}] $$

(see Figure 4.2), where an anticline is located in the middle. In depth, a high permeability zone whose values $\kappa_1$ belong to $[10^{-13}, 10^{-12}]$ is surrounded by two impermeable over- and under- burdens where the permeability $\kappa_2$ is in the range $[10^{-17}, 10^{-15}]$. To represent this geometry, a Corner Point Grid (CPG) with hexahedra and non-planar faces, is used. The number of cells $N$ is equal to 15210. A well is located in the center of $\Omega$ and perforated along 27 cells whose centers lie within the bounding box $[945.819, 1049.83] \times [946.32, 1049.62] \times [-715.73, -537.624]$. The boundary $\Gamma_{int}$ is given by

$$ 
\begin{aligned}
\Gamma_{int} = & \{818.87\} \times [818.87, 1125.95] \times [-816.148, -490.622] \\
& \cup \{1125.95\} \times [818.87, 1125.95] \times [-816.148, -490.622] \\
& \cup [818.87, 1125.95] \times \{818.87\} \times [-816.148, -490.622] \\
& \cup [818.87, 1125.95] \times \{1125.95\} \times [-816.148, -490.622] \\
& \cup [818.87, 1125.95] \times [818.87, 1125.95] \times \{-816.148\} \\
& \cup [818.87, 1125.95] \times [818.87, 1125.95] \times \{-490.622\}.
\end{aligned}
$$

To apply the EIM, SCM and Greedy processes, a sample of parameter values $\Xi_{training} = \{\xi_1, \dots, \xi_L\}$ is generated by randomly choosing $\xi = \{\kappa_1, \kappa_2\}$ from their ranges and by taking $L = 100$. The distribution of the values is shown in Figure 4.1. $\Xi_{training}$ is used in the offline stage and a new sampling set $\Xi_{test}$ is introduced in the online stage to validate the previous processes and, in particular, control the quality of the EIM and SCM. $\Xi_{test}$ is constructed using unexplored parameters $\kappa_1$ and $\kappa_2$ from the same range of values given above.

<figure><img src="image_3.png" /><figcaption>Figure 4.1: Permeabilities’ distribution used for the offline stage</figcaption></figure>
# 4 Numerical results

In this section, we numerically validate the theoretical results obtained for the reduction of problem (2.1). Our main goals are to study both efficiency and computation cost of the proposed estimators.

The following parameters are considered:

$$ 
\begin{array}{l}
\mu = 1.5 \times 10^{-5} \text{ Pa.s,} \quad c_t = 1.4 \times 10^{-7} \text{ Pa}^{-1}, \quad g = 9.81 \text{ m/s}^2, \quad \rho = 700 \text{ kg/m}^3, \\
p_{bh} = 4.13 \times 10^7 \text{ Pa,} \quad \phi = 0.2, \quad r_w = 0.1, \quad z_{bh} = 0 \text{ m.}
\end{array}
$$

We use $p_D = 10^5$ Pa as Dirichlet boundary condition. The total duration of the simulation is $T = 200$ days and the time step is $\Delta t = 10$ days. The initial pressure is defined by

$$ p_K^0 = p_D - \rho g (z_K - z_D), $$

where $z_D = 80$ m.

We consider a three-dimensional domain

$$ \Omega = [-2.686 \cdot 10^{-3} \text{ m}, 1996 \text{ m}] \times [6.1 \cdot 10^{-5} \text{ m}, 1996 \text{ m}] \times [-1000.13 \text{ m}, 2.686 \cdot 10^{-3} \text{ m}] $$

(see Figure 4.2), where an anticline is located in the middle. In depth, a high permeability zone whose values $\kappa_1$ belong to $[10^{-13}, 10^{-12}]$ is surrounded by two impermeable over- and under- burdens where the permeability $\kappa_2$ is in the range $[10^{-17}, 10^{-15}]$. To represent this geometry, a Corner Point Grid (CPG) with hexahedra and non-planar faces, is used. The number of cells $\mathcal{N}$ is equal to 15210. A well is located in the center of $\Omega$ and perforated along 27 cells whose centers lie within the bounding box $[945.819, 1049.83] \times [946.32, 1049.62] \times [-715.73, -537.624]$. The boundary $\Gamma_{\text{int}}$ is given by

$$ 
\begin{aligned}
\Gamma_{\text{int}} = & \{818.87\} \times [818.87, 1125.95] \times [-816.148, -490.622] \\
& \cup \{1125.95\} \times [818.87, 1125.95] \times [-816.148, -490.622] \\
& \cup [818.87, 1125.95] \times \{818.87\} \times [-816.148, -490.622] \\
& \cup [818.87, 1125.95] \times \{1125.95\} \times [-816.148, -490.622] \\
& \cup [818.87, 1125.95] \times [818.87, 1125.95] \times \{-816.148\} \\
& \cup [818.87, 1125.95] \times [818.87, 1125.95] \times \{-490.622\}.
\end{aligned}
$$

To apply the EIM, SCM and Greedy processes, a sample of parameter values $\Xi_{\text{training}} = \{\xi_1, \dots, \xi_\mathcal{L}\}$ is generated by randomly choosing $\xi = \{\kappa_1, \kappa_2\}$ from their ranges and by taking $\mathcal{L} = 100$. The distribution of the values is shown in Figure 4.1. $\Xi_{\text{training}}$ is used in the offline stage and a new sampling set $\Xi_{\text{test}}$ is introduced in the online stage to validate the previous processes and, in particular, control the quality of the EIM and SCM. $\Xi_{\text{test}}$ is constructed using unexplored parameters $\kappa_1$ and $\kappa_2$ from the same range of values given above.

<figure><img src="image_5.png" /><figcaption>Figure 4.1: Permeabilities’ distribution used for the offline stage</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 4.2: Spatial repartition of the permeabilities within Ω: the yellow zone includes cells having the high permeability value κ₁ and the low permeability value κ₂ is used in the blue zone.</figcaption></figure>

## 4.1 Affine decomposition of the scheme coefficients

To construct the reduced model, we start by applying the EIM as discussed in Section 3.5. We plot in Figure 4.3 the evolution of the interpolation error defined as

$$e_{M,\max}^{\infty} = \max_{\xi \in \Xi_{\text{training}}} \frac{\|\hat{v}(\xi) - \mathcal{I}_{M}[\hat{v}(\xi)]\|_{L^{\infty}}}{\|\hat{v}(\xi)\|_{L^{\infty}}}, \quad (4.1)$$

with respect to the number of parameters M. The final number of selected parameters is M<sub>EIM</sub> = 10.

<figure><img src="image_5.png" /><figcaption>Figure 4.3: Evolution of the EIM interpolation error (4.1) with respect to the number of selected parameters M.</figcaption></figure>

To evaluate the accuracy of the linearization with M<sub>EIM</sub> = 10, we compute the maximum relative error for new values of ξ ∈ Ξ<sub>test</sub>. On this sampling, we obtained

$$\max_{\xi \in \Xi_{\text{test}}} \frac{\|\hat{v}(\xi) - \mathcal{I}_{M_{\text{EIM}}}[\hat{v}(\xi)]\|_{L^{\infty}}}{\|\hat{v}(\xi)\|_{L^{\infty}}} = 4.25 \cdot 10^{-16}.$$

## 4.2 Bounds on the coercivity constants

We consider the successive constraint method to compute a lower bound for the coercivity constant α<sub>A<sub>sym</sub></sub> defined in (3.7). We use the same training set Ξ<sub>training</sub> and the affine decomposition obtained with the EIM. For this test, we have used M₁ = M₂ = 5 and tol = 1e-04 (see Appendix B). Since M<sub>EIM</sub> = 10, we have to solve 10 eigenvalue problems to define B. The offline greedy algorithm 5 generates a parameter set Ξ<sub>M</sub> of dimension M = 20. Again to check the quality of the SCM result in the online stage, we compute lower and upper bounds for α<sub>A<sub>sym</sub></sub> for both samplings. More precisely, we compute the values of the ratio

$$r_{A_{\text{sym}}}(\xi) = \frac{\alpha_{A_{\text{sym}}}(\xi) - \alpha_{A_{\text{sym}}, \text{LB}}(\xi)}{\alpha_{A_{\text{sym}}, \text{UB}}(\xi) - \alpha_{A_{\text{sym}}, \text{LB}}(\xi)},$$

and observe that r<sub>A<sub>sym</sub></sub> is in the range [0.999, 1.00479] for all ξ belonging to Ξ<sub>training</sub> and Ξ<sub>test</sub>.
<figure><img src="image_6.png" /><figcaption>Figure 4.2: Spatial repartition of the permeabilities within Ω: the yellow zone includes cells having the high permeability value κ₁ and the low permeability value κ₂ is used in the blue zone.</figcaption></figure>

## 4.1 Affine decomposition of the scheme coefficients

To construct the reduced model, we start by applying the EIM as discussed in Section 3.5. We plot in Figure 4.3 the evolution of the interpolation error defined as

$$e_{M,\max}^{\infty} = \max_{\xi \in \Xi_{\text{training}}} \frac{\|\hat{v}(\xi) - \mathcal{I}_{M}[\hat{v}(\xi)]\|_{L^{\infty}}}{\|\hat{v}(\xi)\|_{L^{\infty}}}, \quad (4.1)$$

with respect to the number of parameters M. The final number of selected parameters is M<sub>EIM</sub> = 10.

<figure><img src="image_7.png" /><figcaption>Figure 4.3: Evolution of the EIM interpolation error (4.1) with respect to the number of selected parameters M.</figcaption></figure>

To evaluate the accuracy of the linearization with M<sub>EIM</sub> = 10, we compute the maximum relative error for new values of ξ ∈ Ξ<sub>test</sub>. On this sampling, we obtained

$$\max_{\xi \in \Xi_{\text{test}}} \frac{\|\hat{v}(\xi) - \mathcal{I}_{M_{\text{EIM}}}[\hat{v}(\xi)]\|_{L^{\infty}}}{\|\hat{v}(\xi)\|_{L^{\infty}}} = 4.25 \cdot 10^{-16}.$$

## 4.2 Bounds on the coercivity constants

We consider the successive constraint method to compute a lower bound for the coercivity constant α<sub>A<sub>sym</sub></sub> defined in (3.7). We use the same training set Ξ<sub>training</sub> and the affine decomposition obtained with the EIM. For this test, we have used M₁ = M₂ = 5 and tol = 1e-04 (see Appendix B). Since M<sub>EIM</sub> = 10, we have to solve 10 eigenvalue problems to define B. The offline greedy algorithm 5 generates a parameter set Ξ<sub>M</sub> of dimension M = 20. Again to check the quality of the SCM result in the online stage, we compute lower and upper bounds for α<sub>A<sub>sym</sub></sub> for both samplings. More precisely, we compute the values of the ratio

$$r_{A_{\text{sym}}}(\xi) = \frac{\alpha_{A_{\text{sym}}}(\xi) - \alpha_{A_{\text{sym}}, \text{LB}}(\xi)}{\alpha_{A_{\text{sym}}, \text{UB}}(\xi) - \alpha_{A_{\text{sym}}, \text{LB}}(\xi)},$$

and observe that r<sub>A<sub>sym</sub></sub> is in the range [0.999, 1.00479] for all ξ belonging to Ξ<sub>training</sub> and Ξ<sub>test</sub>.
## 4.3 POD-Greedy algorithm

As a first test, we construct the reduced model obtained with a POD-Greedy algorithm and the a posteriori estimator $\Delta_{pr}^N$. We define the following errors

$$e_{pr,max}^N = \max_{\xi \in \Xi} \|e^N\|_{pr}, \quad \Delta_{pr,max}^N = \max_{\xi \in \Xi} \Delta_{pr}^N, \quad \eta_{pr,max}^N = \max_{\xi \in \Xi} \frac{\Delta_{pr}^N}{\|e^N\|_{pr}},$$

where $e^N$ and $\Delta_{pr,max}^N$ are given by (3.4) and (3.10) respectively. Figure 4.4 shows the evolution of the a posteriori error estimator $\Delta_{pr,max}^N$ along with the true error $e_{pr,max}^N$ with respect to the basis dimension $N_{pr}$. In Figure 4.4a, the training parameter set $\Xi_{training}$ is used to evaluate these error indicators, while $\Xi_{test}$ is used in Figure 4.4b following the same sequence of introduction of basis vectors as in POD-Greedy process. The results confirm that the proposed estimator is reliable as it forms an upper bound of the true error in both cases. In the offline phase, the POD-Greedy algorithm generates a basis of dimension $N_{pr} = 92$, for which the maximum relative error defined as

$$E_{pr,max}^N = \max_{\xi \in \Xi} \frac{\|e^N\|_{pr}}{\|p_M^N\|_{pr}},$$

reaches $4 \cdot 10^{-10}$. To analyse the efficiency of the estimator, we detail in Table 1 the value of the effectivities $\eta_{pr,max}^N$ in the offline stage for different basis dimensions. We can see that the effectivities are quite good $O(3)$ and the estimator can be safely used to replace the true error.

<figure><img src="image_6.png" /><figcaption>Figure 4.4: Maximum true and estimated errors as functions of the basis dimension for the primal problem using the parameter samplings $\Xi_{training}$ on the left and $\Xi_{test}$ on the right.</figcaption></figure>

<table><caption>Table 1: Effectivities of the primal a posteriori error estimate with respect to the basis dimension $N_{pr}$ in the offline stage.</caption><thead><tr><th>N<sub>pr</sub></th><th>E<sub>pr,max</sub><sup>N</sup></th><th>η<sub>pr,max</sub><sup>N</sup></th></tr></thead><tbody><tr><td>5</td><td>1.04e-03</td><td>658.1</td></tr><tr><td>13</td><td>7.94e-05</td><td>1061.3</td></tr><tr><td>23</td><td>5.05e-05</td><td>2327.8</td></tr><tr><td>33</td><td>9.12e-07</td><td>1549.5</td></tr><tr><td>43</td><td>2.81e-07</td><td>1400.9</td></tr><tr><td>54</td><td>7.39e-08</td><td>1376.7</td></tr><tr><td>68</td><td>1.08e-08</td><td>677.4</td></tr><tr><td>80</td><td>7.46e-09</td><td>447.1</td></tr><tr><td>92</td><td>4e-10</td><td>392</td></tr></tbody></table>
## 4.3 POD-Greedy algorithm

As a first test, we construct the reduced model obtained with a POD-Greedy algorithm and the a posteriori estimator $\Delta_{pr}^N$. We define the following errors

$$e_{pr,max}^N = \max_{\xi \in \Xi} \|e^N\|_{pr}, \quad \Delta_{pr,max}^N = \max_{\xi \in \Xi} \Delta_{pr}^N, \quad \eta_{pr,max}^N = \max_{\xi \in \Xi} \frac{\Delta_{pr}^N}{\|e^N\|_{pr}},$$

where $e^N$ and $\Delta_{pr,max}^N$ are given by (3.4) and (3.10) respectively. Figure 4.4 shows the evolution of the a posteriori error estimator $\Delta_{pr,max}^N$ along with the true error $e_{pr,max}^N$ with respect to the basis dimension $N_{pr}$. In Figure 4.4a, the training parameter set $\Xi_{training}$ is used to evaluate these error indicators, while $\Xi_{test}$ is used in Figure 4.4b following the same sequence of introduction of basis vectors as in POD-Greedy process. The results confirm that the proposed estimator is reliable as it forms an upper bound of the true error in both cases. In the offline phase, the POD-Greedy algorithm generates a basis of dimension $N_{pr} = 92$, for which the maximum relative error defined as

$$E_{pr,max}^N = \max_{\xi \in \Xi} \frac{\|e^N\|_{pr}}{\|p_M^N\|_{pr}},$$

reaches $4 \cdot 10^{-10}$. To analyse the efficiency of the estimator, we detail in Table 1 the value of the effectivities $\eta_{pr,max}^N$ in the offline stage for different basis dimensions. We can see that the effectivities are quite good $\mathcal{O}(3)$ and the estimator can be safely used to replace the true error.

<figure><img src="image_8.png" /><figcaption>Figure 4.4: Maximum true and estimated errors as functions of the basis dimension for the primal problem using the parameter samplings $\Xi_{training}$ on the left and $\Xi_{test}$ on the right.</figcaption></figure>

<table><caption>Table 1: Effectivities of the primal a posteriori error estimate with respect to the basis dimension $N_{pr}$ in the offline stage.</caption><thead><tr><th>N<sub>pr</sub></th><th>E<sup>N</sup><sub>pr,max</sub></th><th>η<sup>N</sup><sub>pr,max</sub></th></tr></thead><tbody><tr><td>5</td><td>1.04e-03</td><td>658.1</td></tr><tr><td>13</td><td>7.94e-05</td><td>1061.3</td></tr><tr><td>23</td><td>5.05e-05</td><td>2327.8</td></tr><tr><td>33</td><td>9.12e-07</td><td>1549.5</td></tr><tr><td>43</td><td>2.81e-07</td><td>1400.9</td></tr><tr><td>54</td><td>7.39e-08</td><td>1376.7</td></tr><tr><td>68</td><td>1.08e-08</td><td>677.4</td></tr><tr><td>80</td><td>7.46e-09</td><td>447.1</td></tr><tr><td>92</td><td>4e-10</td><td>392</td></tr></tbody></table>
In addition, we compare in Figure 4.5, the evolution of $e_{pr,max}^N$ with respect to the basis dimension $N_{pr}$ using a POD-Greedy algorithm driven by different choices of the a posteriori error estimator, which are the ones presented in [21] and defined as

$$ \bar{\Delta}_{pr,1,max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{A_{sym}}} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2}, \quad \text{for } \mathbf{G}^* = \mathbf{A}^*, \qquad (4.2) $$

and

$$ \bar{\Delta}_{pr,2,max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{A_{sym}}} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2}, \quad \text{for } \mathbf{G}^* = \mathbf{M} + \Delta t \mathbf{A}^*. \qquad (4.3) $$

For $N_{pr} < 60$, $e_{pr,max}^N$ behaves in the same way as for the three choices. For $N_{pr} > 60$, a slight difference appears between the three curves. This trend occurs for both $\Xi_{training}$ and $\Xi_{test}$.

<figure><img src="image_7.png" /><figcaption>Figure 4.5: Maximum true error as a function of the primal basis dimension using $\Xi_{training}$ on the left and $\Xi_{test}$ on the right for a POD-Greedy algorithm driven by (3.10), (4.2) and (4.3).</figcaption></figure>

Next, we build a reduced output by first considering the choice (3.21). We use a POD-Greedy algorithm detailed in Algorithm 4 along with the a posteriori error estimator (3.22). We compare, in that case, the evolution of $e_{s,max}^N$ and $\Delta_{s,max}^N$ defined as

$$ e_s^N = |s^N - s^{N_s,N}|, \quad e_{s,max}^N = \max_{\xi \in \Xi} |s^N - s^{N_s,N}|, \quad \Delta_{s,max}^N = \max_{\xi \in \Xi} \Delta_s^N, $$

with respect to the primal basis dimension $N_{pr}$ using $\Xi_{training}$ and $\Xi_{test}$. The results are given in Figure 4.6. We notice that, although the reliability of the estimator is verified, $\Delta_s^N$ is not efficient and the effectivities

$$ \eta_{s,max}^N = \max_{\xi \in \Xi} \frac{\Delta_s^N}{e_s^N} $$

are quite large as shown in Table 2.
In addition, we compare in Figure 4.5, the evolution of $e_{pr,max}^N$ with respect to the basis dimension $N_{pr}$ using a POD-Greedy algorithm driven by different choices of the a posteriori error estimator, which are the ones presented in [21] and defined as

$$ \bar{\Delta}_{pr,1,max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{A_{sym}}} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2}, \quad \text{for } \mathbf{G}^* = \mathbf{A}^*, \qquad (4.2) $$

and

$$ \bar{\Delta}_{pr,2,max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{A_{sym}}} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2}, \quad \text{for } \mathbf{G}^* = \mathbf{M} + \Delta t \mathbf{A}^*. \qquad (4.3) $$

For $N_{pr} < 60$, $e_{pr,max}^N$ behaves in the same way as for the three choices. For $N_{pr} > 60$, a slight difference appears between the three curves. This trend occurs for both $\Xi_{training}$ and $\Xi_{test}$.

<figure><img src="image_9.png" /><figcaption>Figure 4.5: Maximum true error as a function of the primal basis dimension using $\Xi_{training}$ on the left and $\Xi_{test}$ on the right for a POD-Greedy algorithm driven by (3.10), (4.2) and (4.3).</figcaption></figure>

Next, we build a reduced output by first considering the choice (3.21). We use a POD-Greedy algorithm detailed in Algorithm 4 along with the a posteriori error estimator (3.22). We compare, in that case, the evolution of $e_{s,max}^N$ and $\Delta_{s,max}^N$ defined as

$$ e_s^N = |s^N - s^{N_s,N}|, \quad e_{s,max}^N = \max_{\xi \in \Xi} |s^N - s^{N_s,N}|, \quad \Delta_{s,max}^N = \max_{\xi \in \Xi} \Delta_s^N, $$

with respect to the primal basis dimension $N_{pr}$ using $\Xi_{training}$ and $\Xi_{test}$. The results are given in Figure 4.6. We notice that, although the reliability of the estimator is verified, $\Delta_s^N$ is not efficient and the effectivities

$$ \eta_{s,max}^N = \max_{\xi \in \Xi} \frac{\Delta_s^N}{e_s^N} $$

are quite large as shown in Table 2.
<figure><img src="image_8.png" /><figcaption>Figure 4.6: Maximum true and estimated errors for the first choice of the reduced output as functions of the primal basis dimension using the parameter samplings $\Xi_{\text{training}}$ on the left and $\Xi_{\text{test}}$ on the right.</figcaption></figure>

<table><caption>Table 2: Effectivities obtained with the posteriori error estimator (3.22) with respect to the primal and dual basis dimensions $N_{\text{pr}}$ and $N_{\text{du}}$ in the offline stage.</caption><thead><tr><th>N<sub>pr</sub></th><th>N<sub>du</sub></th><th>η<sup>N</sup><sub>s,max</sub></th></tr></thead><tbody><tr><td>5</td><td>10</td><td>2.11e+06</td></tr><tr><td>13</td><td>19</td><td>8.08e+06</td></tr><tr><td>24</td><td>32</td><td>9.83e+06</td></tr><tr><td>35</td><td>44</td><td>1.27e+08</td></tr><tr><td>47</td><td>56</td><td>7.21e+06</td></tr><tr><td>59</td><td>69</td><td>4.95e+06</td></tr><tr><td>73</td><td>84</td><td>1.47e+08</td></tr><tr><td>86</td><td>98</td><td>2.72e+08</td></tr><tr><td>103</td><td>114</td><td>1.08e+08</td></tr></tbody></table>

We also observe that $\eta_{s,max}^N$ becomes less efficient as the final simulation time is increased from $T = 10$ days to $T = 100$ days (see Figure 4.7).

<figure><img src="image_9.png" /><figcaption>Figure 4.7: Maximum true and estimated errors for the first choice of the reduced output as functions of the primal basis dimension with $T = 10$ days (on the left) and $T = 100$ days (on the right).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 4.6: Maximum true and estimated errors for the first choice of the reduced output as functions of the primal basis dimension using the parameter samplings $\Xi_{\text{training}}$ on the left and $\Xi_{\text{test}}$ on the right.</figcaption></figure>

<table><caption>Table 2: Effectivities obtained with the posteriori error estimator (3.22) with respect to the primal and dual basis dimensions $N_{\text{pr}}$ and $N_{\text{du}}$ in the offline stage.</caption><thead><tr><th>N<sub>pr</sub></th><th>N<sub>du</sub></th><th>η<sup>N</sup><sub>s,max</sub></th></tr></thead><tbody><tr><td>5</td><td>10</td><td>2.11e+06</td></tr><tr><td>13</td><td>19</td><td>8.08e+06</td></tr><tr><td>24</td><td>32</td><td>9.83e+06</td></tr><tr><td>35</td><td>44</td><td>1.27e+08</td></tr><tr><td>47</td><td>56</td><td>7.21e+06</td></tr><tr><td>59</td><td>69</td><td>4.95e+06</td></tr><tr><td>73</td><td>84</td><td>1.47e+08</td></tr><tr><td>86</td><td>98</td><td>2.72e+08</td></tr><tr><td>103</td><td>114</td><td>1.08e+08</td></tr></tbody></table>

We also observe that $\eta_{s,max}^N$ becomes less efficient as the final simulation time is increased from $T = 10$ days to $T = 100$ days (see Figure 4.7).

<figure><img src="image_11.png" /><figcaption>Figure 4.7: Maximum true and estimated errors for the first choice of the reduced output as functions of the primal basis dimension with $T = 10$ days (on the left) and $T = 100$ days (on the right).</figcaption></figure>
We also consider the output estimator presented in [18] as

$$ \bar{\Delta}_{s,1,\max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\boldsymbol{\varrho}^n\|_{-1}^2 \right)^{1/2}, \quad \text{for } \mathbf{G}^* = \mathbf{A}^*, \qquad (4.4) $$

and

$$ \bar{\Delta}_{s,2,\max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\boldsymbol{\varrho}^n\|_{-1}^2 \right)^{1/2}, \quad \text{for } \mathbf{G}^* = \mathbf{M} + \Delta t \mathbf{A}^*, \qquad (4.5) $$

and compare in Figure 4.8 the evolution of $\tilde{e}_{s,\max}^N$ as a function of primal basis dimension $N_{pr}$ using a POD-Greedy algorithm controlled by (3.22), (4.4) and (4.5). We observe that the true error is exactly the same using (3.22) and (4.5).

<figure><img src="image_10.png" /><figcaption>Figure 4.8: Estimated errors for the first choice of the reduced output as a function of the primal basis dimension for a POD-Greedy algorithm driven by (3.22), (4.4) and (4.5) using $\Xi_{\text{training}}$ on the left and $\Xi_{\text{test}}$ on the right.</figcaption></figure>

We now employ the second definition of the output error (3.23) and estimator (3.24) and analyze the evolution of the POD-Greedy algorithm in the offline and online stages. The results are shown in Figure 4.9. We define the following errors

$$ \tilde{e}_{s}^{N} = |s^{N} - \tilde{s}^{N_{s},N}|, \quad \tilde{e}_{s,\max}^{N} = \max_{\xi \in \Xi} |s^{N} - \tilde{s}^{N_{s},N}|, \quad \tilde{\Delta}_{s,\max}^{N} = \max_{\xi \in \Xi} \tilde{\Delta}_{s}^{N}, \quad \tilde{\eta}_{s,\max}^{N} = \max_{\xi \in \Xi} \frac{\tilde{\Delta}_{s}^{N}}{\tilde{e}_{s}^{N}}. $$

Table 3 shows that the proposed estimator is very close to the true error and behaves better in terms of effectivity compared to $\eta_{s,\max}^N$: indeed, from $N_{pr} = 59$, $\tilde{\eta}_{s,\max}^N$ behaves as $O(1)$. However, we observe that the first choice of the reduced output (3.21) and the corresponding a posteriori error estimate (3.22) provide better results in terms of accuracy and size of the basis dimensions: to obtain a precision of $1e-10$, we need a primal basis of dimension $N_{pr} = 24$ and a dual basis of dimension $N_{du} = 32$ with the first choice while the second choice requires bases whose sizes are $N_{pr} = 72$ and $N_{du} = 84$. Finally, we introduce the output estimator presented in [18] as

$$ \bar{\bar{\Delta}}_{s,1,\max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\boldsymbol{\varrho}^n\|_{-1}^2 \right)^{1/2} + \Delta t \sum_{n=0}^{N-1} |\langle \mathbf{r}^{n+1}, \Psi^{N_{\text{du}},n} \rangle|, \quad \text{for } \mathbf{G}^* = \mathbf{A}^*, \qquad (4.6) $$

and

$$ \bar{\bar{\Delta}}_{s,2,\max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\boldsymbol{\varrho}^n\|_{-1}^2 \right)^{1/2} + \Delta t \sum_{n=0}^{N-1} |\langle \mathbf{r}^{n+1}, \Psi^{N_{\text{du}},n} \rangle|, \quad \text{for } \mathbf{G}^* = \mathbf{M} + \Delta t \mathbf{A}^*, \qquad (4.7) $$

and compare, in Figure 4.10, the evolution of $\tilde{e}_{s,\max}^N$ using a POD-Greedy algorithm controlled by (3.24), (4.6) and (4.7). We observe that the accuracy is the same when employing the estimators (3.24) and (4.7) and that the curve related to (4.6) lies above the other curves for $N_{pr} < 70$ and below them for $N_{pr} > 70$.
We also consider the output estimator presented in [18] as

$$ \bar{\Delta}_{s,1,\max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\boldsymbol{\varrho}^n\|_{-1}^2 \right)^{1/2}, \quad \text{for } \mathbf{G}^* = \mathbf{A}^*, \qquad (4.4) $$

and

$$ \bar{\Delta}_{s,2,\max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\boldsymbol{\varrho}^n\|_{-1}^2 \right)^{1/2}, \quad \text{for } \mathbf{G}^* = \mathbf{M} + \Delta t \mathbf{A}^*, \qquad (4.5) $$

and compare in Figure 4.8 the evolution of $\tilde{e}_{s,\max}^N$ as a function of primal basis dimension $N_{pr}$ using a POD-Greedy algorithm controlled by (3.22), (4.4) and (4.5). We observe that the true error is exactly the same using (3.22) and (4.5).

<figure><img src="image_12.png" /><figcaption>Figure 4.8: Estimated errors for the first choice of the reduced output as a function of the primal basis dimension for a POD-Greedy algorithm driven by (3.22), (4.4) and (4.5) using $\Xi_{\text{training}}$ on the left and $\Xi_{\text{test}}$ on the right.</figcaption></figure>

We now employ the second definition of the output error (3.23) and estimator (3.24) and analyze the evolution of the POD-Greedy algorithm in the offline and online stages. The results are shown in Figure 4.9. We define the following errors

$$ \tilde{e}_{s}^{N} = |s^{N} - \tilde{s}^{N_{s},N}|, \quad \tilde{e}_{s,\max}^{N} = \max_{\xi \in \Xi} |s^{N} - \tilde{s}^{N_{s},N}|, \quad \tilde{\Delta}_{s,\max}^{N} = \max_{\xi \in \Xi} \tilde{\Delta}_{s}^{N}, \quad \tilde{\eta}_{s,\max}^{N} = \max_{\xi \in \Xi} \frac{\tilde{\Delta}_{s}^{N}}{\tilde{e}_{s}^{N}}. $$

Table 3 shows that the proposed estimator is very close to the true error and behaves better in terms of effectivity compared to $\eta_{s,\max}^N$: indeed, from $N_{pr} = 59$, $\tilde{\eta}_{s,\max}^N$ behaves as $O(1)$. However, we observe that the first choice of the reduced output (3.21) and the corresponding a posteriori error estimate (3.22) provide better results in terms of accuracy and size of the basis dimensions: to obtain a precision of $1e-10$, we need a primal basis of dimension $N_{pr} = 24$ and a dual basis of dimension $N_{du} = 32$ with the first choice while the second choice requires bases whose sizes are $N_{pr} = 72$ and $N_{du} = 84$. Finally, we introduce the output estimator presented in [18] as

$$ \bar{\bar{\Delta}}_{s,1,\max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\boldsymbol{\varrho}^n\|_{-1}^2 \right)^{1/2} + \Delta t \sum_{n=0}^{N-1} |\langle \mathbf{r}^{n+1}, \Psi^{N_{\text{du}},n} \rangle|, \quad \text{for } \mathbf{G}^* = \mathbf{A}^*, \qquad (4.6) $$

and

$$ \bar{\bar{\Delta}}_{s,2,\max} = \left( \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 \sum_{n=1}^{N} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\boldsymbol{\varrho}^n\|_{-1}^2 \right)^{1/2} + \Delta t \sum_{n=0}^{N-1} |\langle \mathbf{r}^{n+1}, \Psi^{N_{\text{du}},n} \rangle|, \quad \text{for } \mathbf{G}^* = \mathbf{M} + \Delta t \mathbf{A}^*, \qquad (4.7) $$

and compare, in Figure 4.10, the evolution of $\tilde{e}_{s,\max}^N$ using a POD-Greedy algorithm controlled by (3.24), (4.6) and (4.7). We observe that the accuracy is the same when employing the estimators (3.24) and (4.7) and that the curve related to (4.6) lies above the other curves for $N_{pr} < 70$ and below them for $N_{pr} > 70$.
<figure><img src="image_11.png" /><figcaption>Figure 4.9: Maximum true and estimated errors for the second choice of the reduced output as functions of the primal basis dimension using $\Xi_{\text{training}}$ on the left and $\Xi_{\text{test}}$ on the right.</figcaption></figure>

<table><caption>Table 3: Effectivities for the a posteriori error estimator (3.24) with respect to the primal and dual basis dimensions $N_{\text{pr}}$ and $N_{\text{du}}$ in the offline stage.</caption><thead><tr><th>N<sub>pr</sub></th><th>N<sub>du</sub></th><th>η̃<sup>N</sup><sub>s,max</sub></th></tr></thead><tbody><tr><td>5</td><td>10</td><td>103103</td></tr><tr><td>13</td><td>19</td><td>222673</td></tr><tr><td>24</td><td>32</td><td>2066.3</td></tr><tr><td>35</td><td>44</td><td>23.6</td></tr><tr><td>47</td><td>56</td><td>1266.5</td></tr><tr><td>59</td><td>69</td><td>1.26</td></tr><tr><td>72</td><td>84</td><td>1.06</td></tr><tr><td>85</td><td>99</td><td>1.34</td></tr><tr><td>100</td><td>115</td><td>1.18</td></tr><tr><td>113</td><td>129</td><td>1.0</td></tr><tr><td>125</td><td>143</td><td>1.01</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Figure 4.9: Maximum true and estimated errors for the second choice of the reduced output as functions of the primal basis dimension using $\Xi_{\text{training}}$ on the left and $\Xi_{\text{test}}$ on the right.</figcaption></figure>

<table><caption>Table 3: Effectivities for the a posteriori error estimator (3.24) with respect to the primal and dual basis dimensions $N_{\text{pr}}$ and $N_{\text{du}}$ in the offline stage.</caption><thead><tr><th>N<sub>pr</sub></th><th>N<sub>du</sub></th><th>η̃<sup>N</sup><sub>s,max</sub></th></tr></thead><tbody><tr><td>5</td><td>10</td><td>103103</td></tr><tr><td>13</td><td>19</td><td>222673</td></tr><tr><td>24</td><td>32</td><td>2066.3</td></tr><tr><td>35</td><td>44</td><td>23.6</td></tr><tr><td>47</td><td>56</td><td>1266.5</td></tr><tr><td>59</td><td>69</td><td>1.26</td></tr><tr><td>72</td><td>84</td><td>1.06</td></tr><tr><td>85</td><td>99</td><td>1.34</td></tr><tr><td>100</td><td>115</td><td>1.18</td></tr><tr><td>113</td><td>129</td><td>1.0</td></tr><tr><td>125</td><td>143</td><td>1.01</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>Figure 4.10: Estimated errors for the second choice of the reduced output as a function of the primal basis dimension for a POD-Greedy algorithm driven by (3.24), (4.6) and (4.7) using $\Xi_{\text{training}}$ on the left and $\Xi_{\text{test}}$ on the right.</figcaption></figure>

## 4.4 Computational time effort

Here we give an indication of the time spent at each stage of the construction and evaluation of the reduced model.

**Offline stage.** We plot in Figure 4.11 the evolution of the computation times related to the different stages of the reduced-basis construction as functions of the primal basis dimension $N_{pr}$. The time is calculated as a factor of the time $u$ required to run one single high-fidelity simulation. Concerning the EIM and SCM algorithms, these two procedures are applied once at the beginning of the offline stage before starting the POD-Greedy process. This explains why their cumulated times appear as constant in the plot. These times amount to $100 \times u$ seconds and $39 \times u$ seconds respectively. The cumulated time spent in the Greedy process is represented in blue up to $N_{pr} = 103$. For each greedy iteration, the given values include the times required to assemble and compute the reduced solutions (3.2) and (3.16) and the residuals (3.27)–(3.30) for all parameters of the sampling. These operations depend on the sizes $N_{pr}$ and $N_{du}$ and therefore increase as the Greedy process evolves. The "Greedy" time therefore includes the time spent in assembling the reduced primal and dual systems (3.3) and (3.17) for all parameters of the sampling. This time is represented under the label "LF assembly" too. It is quite significant in the offline stage since all products involving the terms of the affine decompositions with the new bases matrices should be updated. This cost is of course substantially reduced in the online stage once the bases are fixed. The cumulated calculation time of the POD method is represented in Figure 4.11. It is linear with respect to $N_{pr}$. It includes the times required to run the high-fidelity simulations and the extractions of the POD modes. In both cases, for each selected parameter, these times are roughly constant.
<figure><img src="image_14.png" /><figcaption>Figure 4.10: Estimated errors for the second choice of the reduced output as a function of the primal basis dimension for a POD-Greedy algorithm driven by (3.24), (4.6) and (4.7) using $\Xi_{\text{training}}$ on the left and $\Xi_{\text{test}}$ on the right.</figcaption></figure>

## 4.4 Computational time effort

Here we give an indication of the time spent at each stage of the construction and evaluation of the reduced model.

**Offline stage.** We plot in Figure 4.11 the evolution of the computation times related to the different stages of the reduced-basis construction as functions of the primal basis dimension $N_{\text{pr}}$. The time is calculated as a factor of the time $u$ required to run one single high-fidelity simulation. Concerning the EIM and SCM algorithms, these two procedures are applied once at the beginning of the offline stage before starting the POD-Greedy process. This explains why their cumulated times appear as constant in the plot. These times amount to $100 \times u$ seconds and $39 \times u$ seconds respectively. The cumulated time spent in the Greedy process is represented in blue up to $N_{\text{pr}} = 103$. For each greedy iteration, the given values include the times required to assemble and compute the reduced solutions (3.2) and (3.16) and the residuals (3.27)–(3.30) for all parameters of the sampling. These operations depend on the sizes $N_{\text{pr}}$ and $N_{\text{du}}$ and therefore increase as the Greedy process evolves. The "Greedy" time therefore includes the time spent in assembling the reduced primal and dual systems (3.3) and (3.17) for all parameters of the sampling. This time is represented under the label "LF assembly" too. It is quite significant in the offline stage since all products involving the terms of the affine decompositions with the new bases matrices should be updated. This cost is of course substantially reduced in the online stage once the bases are fixed. The cumulated calculation time of the POD method is represented in Figure 4.11. It is linear with respect to $N_{\text{pr}}$. It includes the times required to run the high-fidelity simulations and the extractions of the POD modes. In both cases, for each selected parameter, these times are roughly constant.
<figure><img src="image_13.png" /><figcaption>Figure 4.11: Offline time computation effort.</figcaption></figure>

**Online stage.** Given a new parameter value $\xi \in \Xi_{test}$, with our implementation, the time used in the online stage to compute $p^{N_{pr},N}$ and its corresponding reduced output at $T = 200$ days and with $N_{pr} = 92$ is divided by a factor of 10 compared to one HF run needed to obtain $p_M^N$. The high-fidelity model (2.1) is solved using a stabilized bi-conjugate gradient method with an incomplete LU preconditioner, where the tolerance is set to $10^{-14}$. The reduced model (3.3) is directly solved using an LU-Decomposition.

# 5 Conclusion

In this work, we have discussed a reduced basis method for (finite volume approximations of) parabolic PDEs. We have introduced a new rigorous a posteriori estimator to evaluate the reduction error in a new discrete space-time energy norm independently of the parameter. We have performed numerical simulations in the context of porous media flows (single-phase flows of slightly compressible fluid parametrized by the permeability) to assess the reliability of the a posteriori error bound and its efficiency at selecting a reduced basis within a POD-Greedy algorithm.

Our numerical results show that our new approach can efficiently reduce the computational cost of engineering studies with many parameter values in the context of porous media flows, especially on choosing well the reduced output in goal-oriented cases with linear QOIs. The discussed methodology can be also considered to estimate different types of linear quantity of interests such as the pressure variation along faults far from the well injection area. Indeed, understanding how injection activities affect pressure distribution in fault networks helps in mitigating risks associated with CO2 migration, fault activation, and potential leakage into overlying aquifers.

# A Proofs of various propositions

## A.1 Proof of Proposition 3.1

*Proof.* For each $v \in \mathbb{R}^N$, we have

$$ \langle (M + \Delta t A)e^n, v \rangle = \langle Me^{n-1}, v \rangle - \Delta t \langle r^n, v \rangle. \quad (A.1) $$

We apply $e^n$ to (A.1). We apply Cauchy-Schwarz inequality and use (3.6). This leads to

$$ \langle (M + \Delta t A)e^n, e^n \rangle \le \|M^{1/2}e^{n-1}\| \|M^{1/2}e^n\| + \Delta t \|r^n\|_{-1} \|e^n\|_{G^*}. \quad (A.2) $$

Now, recalling Young's inequality (for $c \in \mathbb{R}$, $d \in \mathbb{R}$, $\rho \in \mathbb{R}_+$):

$$ 2|c| |d| \le \frac{1}{\rho^2} c^2 + \rho^2 d^2, \quad (A.3) $$
<figure><img src="image_15.png" /><figcaption>Figure 4.11: Offline time computation effort.</figcaption></figure>

**Online stage.** Given a new parameter value $\xi \in \Xi_{\text{test}}$, with our implementation, the time used in the online stage to compute $p^{N_{\text{pr}},N}$ and its corresponding reduced output at $T = 200$ days and with $N_{\text{pr}} = 92$ is divided by a factor of 10 compared to one HF run needed to obtain $p_M^N$. The high-fidelity model (2.1) is solved using a stabilized bi-conjugate gradient method with an incomplete LU preconditioner, where the tolerance is set to $10^{-14}$. The reduced model (3.3) is directly solved using an LU-Decomposition.

# 5 Conclusion

In this work, we have discussed a reduced basis method for (finite volume approximations of) parabolic PDEs. We have introduced a new rigorous a posteriori estimator to evaluate the reduction error in a new discrete space-time energy norm independently of the parameter. We have performed numerical simulations in the context of porous media flows (single-phase flows of slightly compressible fluid parametrized by the permeability) to assess the reliability of the a posteriori error bound and its efficiency at selecting a reduced basis within a POD-Greedy algorithm.

Our numerical results show that our new approach can efficiently reduce the computational cost of engineering studies with many parameter values in the context of porous media flows, especially on choosing well the reduced output in goal-oriented cases with linear QOIs. The discussed methodology can be also considered to estimate different types of linear quantity of interests such as the pressure variation along faults far from the well injection area. Indeed, understanding how injection activities affect pressure distribution in fault networks helps in mitigating risks associated with CO2 migration, fault activation, and potential leakage into overlying aquifers.

# A Proofs of various propositions

## A.1 Proof of Proposition 3.1

*Proof.* For each $v \in \mathbb{R}^N$, we have

$$ \langle (M + \Delta t A)e^n, v \rangle = \langle Me^{n-1}, v \rangle - \Delta t \langle r^n, v \rangle. \quad (A.1) $$

We apply $e^n$ to (A.1). We apply Cauchy-Schwarz inequality and use (3.6). This leads to

$$ \langle (M + \Delta t A)e^n, e^n \rangle \le \|M^{1/2}e^{n-1}\| \|M^{1/2}e^n\| + \Delta t \|r^n\|_{-1} \|e^n\|_{G^*}. \quad (A.2) $$

Now, recalling Young's inequality (for $c \in \mathbb{R}$, $d \in \mathbb{R}$, $\rho \in \mathbb{R}_+$):

$$ 2|c| |d| \le \frac{1}{\rho^2} c^2 + \rho^2 d^2, \quad (A.3) $$
and apply it twice: once for $c = \|\mathbf{M}^{1/2}\mathbf{e}^{n-1}\|$, $d = \|\mathbf{M}^{1/2}\mathbf{e}^n\|$ and $\rho = 1$ to get

$$2\|\mathbf{M}^{1/2}\mathbf{e}^{n-1}\| \|\mathbf{M}^{1/2}\mathbf{e}^n\| \le \langle \mathbf{M}\mathbf{e}^{n-1}, \mathbf{e}^{n-1} \rangle + \langle \mathbf{M}\mathbf{e}^n, \mathbf{e}^n \rangle \quad (A.4)$$

and another time for $c = \|\mathbf{r}^n\|_{-1}$, $d = \|\mathbf{e}^n\|_{G^*}$ and $\rho = \sqrt{\alpha_{\mathbf{A}_{\text{sym}}}}$ to obtain

$$2\|\mathbf{r}^n\|_{-1}\|\mathbf{e}^n\|_{G^*} \le \frac{1}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 + \alpha_{\mathbf{A}_{\text{sym}}} \|\mathbf{e}^n\|_{G^*}^2. \quad (A.5)$$

Now the definition of the coercivity constant (3.7) leads to

$$2\|\mathbf{r}^n\|_{-1} \|\mathbf{e}^n\|_{G^*} \le \frac{1}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 + \langle \mathbf{A}_{\text{sym}} \mathbf{e}^n, \mathbf{e}^n \rangle. \quad (A.6)$$

Combining (A.2), (A.4) and (A.6) yields

$$\langle \mathbf{M}\mathbf{e}^n, \mathbf{e}^n \rangle + \Delta t \langle \mathbf{A}\mathbf{e}^n, \mathbf{e}^n \rangle \le \frac{1}{2} \langle \mathbf{M}\mathbf{e}^{n-1}, \mathbf{e}^{n-1} \rangle + \frac{1}{2} \langle \mathbf{M}\mathbf{e}^n, \mathbf{e}^n \rangle \\ + \frac{\Delta t}{2} \langle \mathbf{A}_{\text{sym}} \mathbf{e}^n, \mathbf{e}^n \rangle + \frac{\Delta t}{2\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2. \quad (A.7)$$

Since $\langle \mathbf{A}_{\text{skew}} \mathbf{e}^n, \mathbf{e}^n \rangle = 0$ and $\langle \mathbf{A}_{\text{sym}} \mathbf{e}^n, \mathbf{e}^n \rangle = \langle \mathbf{A}\mathbf{e}^n, \mathbf{e}^n \rangle$, we obtain

$$\langle \mathbf{M}\mathbf{e}^n, \mathbf{e}^n \rangle - \langle \mathbf{M}\mathbf{e}^{n-1}, \mathbf{e}^{n-1} \rangle + \Delta t \langle \mathbf{A}_{\text{sym}} \mathbf{e}^n, \mathbf{e}^n \rangle \le \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2. \quad (A.8)$$

Finally, we sum (A.8) over $\{1, \dots, n\}$ and consider that $\mathbf{e}^0 = 0$ to get

$$\langle \mathbf{M}\mathbf{e}^n, \mathbf{e}^n \rangle + \sum_{m=1}^{n} \Delta t \langle \mathbf{A}_{\text{sym}} \mathbf{e}^m, \mathbf{e}^m \rangle \le \sum_{m=1}^{n} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^m\|_{-1}^2. \quad (A.9)$$

From (A.9), we have

$$\sum_{n=1}^{N} \langle \mathbf{M} \mathbf{e}^n, \mathbf{e}^n \rangle \le \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{n=1}^{N} (N+1-n) \|\mathbf{r}^n\|_{-1}^2 \le \frac{T}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2. \quad (A.10)$$

As a consequence,

$$\sum_{n=1}^{N} \langle \mathbf{M} \mathbf{e}^n, \mathbf{e}^n \rangle \le \frac{T}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2. \quad (A.11)$$

On the other hand, using (A.9), we can write

$$\sum_{m=1}^{n} \langle \mathbf{A}_{\text{sym}} \mathbf{e}^m, \mathbf{e}^m \rangle \le \frac{1}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=1}^{n} \|\mathbf{r}^m\|_{-1}^2, \quad \forall n \in \{1, \dots, N\}. \quad (A.12)$$

Now (A.11) with (A.12) for $n = N$ enable the following inequality

$$\sum_{m=1}^{N} [\langle \mathbf{M} \mathbf{e}^m, \mathbf{e}^m \rangle + \Delta t \langle \mathbf{A}_{\text{sym}} \mathbf{e}^m, \mathbf{e}^m \rangle] \le \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=1}^{N} \|\mathbf{r}^m\|_{-1}^2 \quad (A.13)$$

where it is possible to use (3.8) and write

$$\alpha_{\mathbf{G},\mathbf{LB}} \sum_{m=1}^{N} [\langle \mathbf{M} \mathbf{e}^m, \mathbf{e}^m \rangle + \Delta t \langle \mathbf{A}_{\text{sym}}^* \mathbf{e}^m, \mathbf{e}^m \rangle] \le \sum_{m=1}^{N} [\langle \mathbf{M} \mathbf{e}^m, \mathbf{e}^m \rangle + \Delta t \langle \mathbf{A}_{\text{sym}} \mathbf{e}^m, \mathbf{e}^m \rangle] \\ \le \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=1}^{N} \|\mathbf{r}^m\|_{-1}^2 \quad (A.14) \\ \le \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym},\mathbf{LB}}}} \sum_{m=1}^{N} \|\mathbf{r}^m\|_{-1}^2$$

i.e. an upper bound of the error that is independent of the parameter $\xi$ as opposed to [21, Prop. 4.3]. $\square$
and apply it twice: once for $c = \|\mathbf{M}^{1/2}\mathbf{e}^{n-1}\|$, $d = \|\mathbf{M}^{1/2}\mathbf{e}^n\|$ and $\rho = 1$ to get

$$2\|\mathbf{M}^{1/2}\mathbf{e}^{n-1}\| \|\mathbf{M}^{1/2}\mathbf{e}^n\| \le \langle \mathbf{M}\mathbf{e}^{n-1}, \mathbf{e}^{n-1} \rangle + \langle \mathbf{M}\mathbf{e}^n, \mathbf{e}^n \rangle \quad (A.4)$$

and another time for $c = \|\mathbf{r}^n\|_{-1}$, $d = \|\mathbf{e}^n\|_{G^*}$ and $\rho = \sqrt{\alpha_{\mathbf{A}_{\text{sym}}}}$ to obtain

$$2\|\mathbf{r}^n\|_{-1}\|\mathbf{e}^n\|_{G^*} \le \frac{1}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 + \alpha_{\mathbf{A}_{\text{sym}}} \|\mathbf{e}^n\|_{G^*}^2. \quad (A.5)$$

Now the definition of the coercivity constant (3.7) leads to

$$2\|\mathbf{r}^n\|_{-1} \|\mathbf{e}^n\|_{G^*} \le \frac{1}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2 + \langle \mathbf{A}_{\text{sym}} \mathbf{e}^n, \mathbf{e}^n \rangle. \quad (A.6)$$

Combining (A.2), (A.4) and (A.6) yields

$$\langle \mathbf{M}\mathbf{e}^n, \mathbf{e}^n \rangle + \Delta t \langle \mathbf{A}\mathbf{e}^n, \mathbf{e}^n \rangle \le \frac{1}{2} \langle \mathbf{M}\mathbf{e}^{n-1}, \mathbf{e}^{n-1} \rangle + \frac{1}{2} \langle \mathbf{M}\mathbf{e}^n, \mathbf{e}^n \rangle \\ + \frac{\Delta t}{2} \langle \mathbf{A}_{\text{sym}} \mathbf{e}^n, \mathbf{e}^n \rangle + \frac{\Delta t}{2\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2. \quad (A.7)$$

Since $\langle \mathbf{A}_{\text{skew}} \mathbf{e}^n, \mathbf{e}^n \rangle = 0$ and $\langle \mathbf{A}_{\text{sym}} \mathbf{e}^n, \mathbf{e}^n \rangle = \langle \mathbf{A}\mathbf{e}^n, \mathbf{e}^n \rangle$, we obtain

$$\langle \mathbf{M}\mathbf{e}^n, \mathbf{e}^n \rangle - \langle \mathbf{M}\mathbf{e}^{n-1}, \mathbf{e}^{n-1} \rangle + \Delta t \langle \mathbf{A}_{\text{sym}} \mathbf{e}^n, \mathbf{e}^n \rangle \le \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^n\|_{-1}^2. \quad (A.8)$$

Finally, we sum (A.8) over $\{1, \dots, n\}$ and consider that $\mathbf{e}^0 = 0$ to get

$$\langle \mathbf{M}\mathbf{e}^n, \mathbf{e}^n \rangle + \sum_{m=1}^{n} \Delta t \langle \mathbf{A}_{\text{sym}} \mathbf{e}^m, \mathbf{e}^m \rangle \le \sum_{m=1}^{n} \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \|\mathbf{r}^m\|_{-1}^2. \quad (A.9)$$

From (A.9), we have

$$\sum_{n=1}^{N} \langle \mathbf{M} \mathbf{e}^n, \mathbf{e}^n \rangle \le \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{n=1}^{N} (N+1-n) \|\mathbf{r}^n\|_{-1}^2 \le \frac{T}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2. \quad (A.10)$$

As a consequence,

$$\sum_{n=1}^{N} \langle \mathbf{M} \mathbf{e}^n, \mathbf{e}^n \rangle \le \frac{T}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2. \quad (A.11)$$

On the other hand, using (A.9), we can write

$$\sum_{m=1}^{n} \langle \mathbf{A}_{\text{sym}} \mathbf{e}^m, \mathbf{e}^m \rangle \le \frac{1}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=1}^{n} \|\mathbf{r}^m\|_{-1}^2, \quad \forall n \in \{1, \dots, N\}. \quad (A.12)$$

Now (A.11) with (A.12) for $n = N$ enable the following inequality

$$\sum_{m=1}^{N} [\langle \mathbf{M} \mathbf{e}^m, \mathbf{e}^m \rangle + \Delta t \langle \mathbf{A}_{\text{sym}} \mathbf{e}^m, \mathbf{e}^m \rangle] \le \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=1}^{N} \|\mathbf{r}^m\|_{-1}^2 \quad (A.13)$$

where it is possible to use (3.8) and write

$$\begin{aligned} \alpha_{\mathbf{G},\mathbf{LB}} \sum_{m=1}^{N} [\langle \mathbf{M} \mathbf{e}^m, \mathbf{e}^m \rangle + \Delta t \langle \mathbf{A}_{\text{sym}}^* \mathbf{e}^m, \mathbf{e}^m \rangle] &\le \sum_{m=1}^{N} [\langle \mathbf{M} \mathbf{e}^m, \mathbf{e}^m \rangle + \Delta t \langle \mathbf{A}_{\text{sym}} \mathbf{e}^m, \mathbf{e}^m \rangle] \\ &\le \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=1}^{N} \|\mathbf{r}^m\|_{-1}^2 \\ &\le \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym},\mathbf{LB}}}} \sum_{m=1}^{N} \|\mathbf{r}^m\|_{-1}^2 \end{aligned} \quad (A.14)$$

i.e. an upper bound of the error that is independent of the parameter $\xi$ as opposed to [21, Prop. 4.3]. $\square$
## A.2 Proof of Proposition 3.2

*Proof.* We start by writing

$$ \langle (M + \Delta t A^T) \boldsymbol{\varepsilon}^m, \boldsymbol{v} \rangle = \langle M \boldsymbol{\varepsilon}^{m+1}, \boldsymbol{v} \rangle - \Delta t \langle \boldsymbol{\varrho}^m, \boldsymbol{v} \rangle, \quad \forall \boldsymbol{v} \in \mathbb{R}^N, $$

and then take $\boldsymbol{v} = \boldsymbol{\varepsilon}^m$ to obtain

$$ \langle (M + \Delta t A^T) \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle = \langle M \boldsymbol{\varepsilon}^{m+1}, \boldsymbol{\varepsilon}^m \rangle - \Delta t \langle \boldsymbol{\varrho}^m, \boldsymbol{\varepsilon}^m \rangle. $$

We apply Cauchy-Schwarz inequality and use (3.19)

$$ \langle (M + \Delta t A^T) \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle \le \|M^{1/2} \boldsymbol{\varepsilon}^{m+1}\| \|M^{1/2} \boldsymbol{\varepsilon}^m\| + \Delta t \|\boldsymbol{\varrho}^m\|_{-1} \|\boldsymbol{\varepsilon}^m\|_{G^*}. \quad (A.15) $$

Similarly to the primal problem, we apply inequality (A.3) twice and get

$$ 2 \|M^{1/2} \boldsymbol{\varepsilon}^{m+1}\| \|M^{1/2} \boldsymbol{\varepsilon}^m\| \le \langle M \boldsymbol{\varepsilon}^{m+1}, \boldsymbol{\varepsilon}^{m+1} \rangle + \langle M \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle, \quad (A.16) $$

$$ 2 \|\boldsymbol{\varrho}^m\|_{-1} \|\boldsymbol{\varepsilon}^m\|_{G^*} \le \frac{1}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2 + \alpha_{A_{\text{sym}}} \|\boldsymbol{\varepsilon}^m\|_{G^*}^2. \quad (A.17) $$

Based on the definition of $\alpha_{A_{\text{sym}}}$, (A.17) becomes

$$ \begin{aligned} 2 \|\boldsymbol{\varrho}^m\|_{-1} \|\boldsymbol{\varepsilon}^m\|_{G^*} &\le \frac{1}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2 + \langle A_{\text{sym}} \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle \\ &= \frac{1}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2 + \langle A \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle, \quad \text{since } \langle A_{\text{skew}} \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle = 0 \\ &= \frac{1}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2 + \langle \boldsymbol{\varepsilon}^m, A^T \boldsymbol{\varepsilon}^m \rangle \\ &= \frac{1}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2 + \langle A^T \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle. \end{aligned} \quad (A.18) $$

Now, inequalities (A.15)-(A.18) lead to

$$ \langle M \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle - \langle M \boldsymbol{\varepsilon}^{m+1}, \boldsymbol{\varepsilon}^{m+1} \rangle + \Delta t \langle A^T \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle \le \frac{\Delta t}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2. \quad (A.19) $$

We finally sum (A.19) over $\{n, \dots, N-1\}$ and suppose that $\boldsymbol{\varepsilon}^N = 0$. We get

$$ \langle \boldsymbol{\varepsilon}^n, M \boldsymbol{\varepsilon}^n \rangle + \Delta t \sum_{m=n}^{N-1} \langle \boldsymbol{\varepsilon}^m, A^T \boldsymbol{\varepsilon}^m \rangle \le \frac{\Delta t}{\alpha_{A_{\text{sym}}}} \sum_{m=n}^{N-1} \|\boldsymbol{\varrho}^m\|_{-1}^2. \quad (A.20) $$

Inequality (A.20) holds true for all $n \in \{0, \dots, N-1\}$. So we can write

$$ \begin{aligned} \sum_{n=0}^{N-1} \langle \boldsymbol{\varepsilon}^n, M \boldsymbol{\varepsilon}^n \rangle &\le \frac{\Delta t}{\alpha_{A_{\text{sym}}}} \sum_{n=0}^{N-1} \sum_{m=n}^{N-1} \|\boldsymbol{\varrho}^m\|_{-1}^2 \\ &\le \frac{\Delta t}{\alpha_{A_{\text{sym}}}} \sum_{n=0}^{N-1} (n+1) \|\boldsymbol{\varrho}^n\|_{-1}^2 \\ &\le \frac{N \Delta t}{\alpha_{A_{\text{sym}}}} \sum_{n=0}^{N-1} \|\boldsymbol{\varrho}^n\|_{-1}^2 \\ &= \frac{T}{\alpha_{A_{\text{sym}}}} \sum_{n=0}^{N-1} \|\boldsymbol{\varrho}^n\|_{-1}^2. \end{aligned} \quad (A.21) $$

Then,

$$ \sum_{n=0}^{N-1} \langle \boldsymbol{\varepsilon}^n, M \boldsymbol{\varepsilon}^n \rangle \le \frac{T}{\alpha_{A_{\text{sym}}}} \sum_{n=0}^{N-1} \|\boldsymbol{\varrho}^n\|_{-1}^2. \quad (A.22) $$
## A.2 Proof of Proposition 3.2

*Proof.* We start by writing

$$ \langle (M + \Delta t A^T) \boldsymbol{\varepsilon}^m, \boldsymbol{v} \rangle = \langle M \boldsymbol{\varepsilon}^{m+1}, \boldsymbol{v} \rangle - \Delta t \langle \boldsymbol{\varrho}^m, \boldsymbol{v} \rangle, \quad \forall \boldsymbol{v} \in \mathbb{R}^N, $$

and then take $\boldsymbol{v} = \boldsymbol{\varepsilon}^m$ to obtain

$$ \langle (M + \Delta t A^T) \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle = \langle M \boldsymbol{\varepsilon}^{m+1}, \boldsymbol{\varepsilon}^m \rangle - \Delta t \langle \boldsymbol{\varrho}^m, \boldsymbol{\varepsilon}^m \rangle. $$

We apply Cauchy-Schwarz inequality and use (3.19)

$$ \langle (M + \Delta t A^T) \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle \le \|M^{1/2} \boldsymbol{\varepsilon}^{m+1}\| \|M^{1/2} \boldsymbol{\varepsilon}^m\| + \Delta t \|\boldsymbol{\varrho}^m\|_{-1} \|\boldsymbol{\varepsilon}^m\|_{G^*}. \quad (A.15) $$

Similarly to the primal problem, we apply inequality (A.3) twice and get

$$ 2 \|M^{1/2} \boldsymbol{\varepsilon}^{m+1}\| \|M^{1/2} \boldsymbol{\varepsilon}^m\| \le \langle M \boldsymbol{\varepsilon}^{m+1}, \boldsymbol{\varepsilon}^{m+1} \rangle + \langle M \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle, \quad (A.16) $$

$$ 2 \|\boldsymbol{\varrho}^m\|_{-1} \|\boldsymbol{\varepsilon}^m\|_{G^*} \le \frac{1}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2 + \alpha_{A_{\text{sym}}} \|\boldsymbol{\varepsilon}^m\|_{G^*}^2. \quad (A.17) $$

Based on the definition of $\alpha_{A_{\text{sym}}}$, (A.17) becomes

$$ \begin{aligned} 2 \|\boldsymbol{\varrho}^m\|_{-1} \|\boldsymbol{\varepsilon}^m\|_{G^*} &\le \frac{1}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2 + \langle A_{\text{sym}} \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle \\ &= \frac{1}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2 + \langle A \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle, \quad \text{since } \langle A_{\text{skew}} \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle = 0 \\ &= \frac{1}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2 + \langle \boldsymbol{\varepsilon}^m, A^T \boldsymbol{\varepsilon}^m \rangle \\ &= \frac{1}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2 + \langle A^T \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle. \end{aligned} \quad (A.18) $$

Now, inequalities (A.15)-(A.18) lead to

$$ \langle M \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle - \langle M \boldsymbol{\varepsilon}^{m+1}, \boldsymbol{\varepsilon}^{m+1} \rangle + \Delta t \langle A^T \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle \le \frac{\Delta t}{\alpha_{A_{\text{sym}}}} \|\boldsymbol{\varrho}^m\|_{-1}^2. \quad (A.19) $$

We finally sum (A.19) over $\{n, \dots, N-1\}$ and suppose that $\boldsymbol{\varepsilon}^N = 0$. We get

$$ \langle \boldsymbol{\varepsilon}^n, M \boldsymbol{\varepsilon}^n \rangle + \Delta t \sum_{m=n}^{N-1} \langle \boldsymbol{\varepsilon}^m, A^T \boldsymbol{\varepsilon}^m \rangle \le \frac{\Delta t}{\alpha_{A_{\text{sym}}}} \sum_{m=n}^{N-1} \|\boldsymbol{\varrho}^m\|_{-1}^2. \quad (A.20) $$

Inequality (A.20) holds true for all $n \in \{0, \dots, N-1\}$. So we can write

$$ \begin{aligned} \sum_{n=0}^{N-1} \langle \boldsymbol{\varepsilon}^n, M \boldsymbol{\varepsilon}^n \rangle &\le \frac{\Delta t}{\alpha_{A_{\text{sym}}}} \sum_{n=0}^{N-1} \sum_{m=n}^{N-1} \|\boldsymbol{\varrho}^m\|_{-1}^2 \\ &\le \frac{\Delta t}{\alpha_{A_{\text{sym}}}} \sum_{n=0}^{N-1} (n+1) \|\boldsymbol{\varrho}^n\|_{-1}^2 \\ &\le \frac{N \Delta t}{\alpha_{A_{\text{sym}}}} \sum_{n=0}^{N-1} \|\boldsymbol{\varrho}^n\|_{-1}^2 \\ &\le \frac{T}{\alpha_{A_{\text{sym}}}} \sum_{n=0}^{N-1} \|\boldsymbol{\varrho}^n\|_{-1}^2. \end{aligned} \quad (A.21) $$

Then,

$$ \sum_{n=0}^{N-1} \langle \boldsymbol{\varepsilon}^n, M \boldsymbol{\varepsilon}^n \rangle \le \frac{T}{\alpha_{A_{\text{sym}}}} \sum_{n=0}^{N-1} \|\boldsymbol{\varrho}^n\|_{-1}^2. \quad (A.22) $$
On the other side, again using (A.20), we have in particular for $n = 0$

$$ \Delta t \sum_{m=0}^{N-1} \langle \boldsymbol{\varepsilon}^m, \mathbf{A}^T \boldsymbol{\varepsilon}^m \rangle \leq \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=0}^{N-1} \| \boldsymbol{\varrho}^m \|_{-1}^2. \quad (A.23) $$

Combining (A.22) and (A.23), leads to

$$ \sum_{m=0}^{N-1} \left[ \langle \boldsymbol{\varepsilon}^m, \mathbf{M} \boldsymbol{\varepsilon}^m \rangle + \Delta t \langle \boldsymbol{\varepsilon}^m, \mathbf{A}^T \boldsymbol{\varepsilon}^m \rangle \right] \leq \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=0}^{N-1} \| \boldsymbol{\varrho}^m \|_{-1}^2. \quad (A.24) $$

Now, since

$$ \langle \boldsymbol{\varepsilon}^m, \mathbf{A}^T \boldsymbol{\varepsilon}^m \rangle = \langle \mathbf{A} \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle = \langle \mathbf{A}_{\text{sym}} \boldsymbol{\varepsilon}^m, \boldsymbol{\varepsilon}^m \rangle, $$

we get

$$ \sum_{m=0}^{N-1} \left[ \langle \boldsymbol{\varepsilon}^m, \mathbf{M} \boldsymbol{\varepsilon}^m \rangle + \Delta t \langle \boldsymbol{\varepsilon}^m, \mathbf{A}_{\text{sym}} \boldsymbol{\varepsilon}^m \rangle \right] \leq \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=0}^{N-1} \| \boldsymbol{\varrho}^m \|_{-1}^2. \quad (A.25) $$

Finally, we recall the definition of $\alpha_G$, which results in

$$ \sum_{m=0}^{N-1} \left[ \langle \boldsymbol{\varepsilon}^m, \mathbf{M} \boldsymbol{\varepsilon}^m \rangle + \Delta t \langle \boldsymbol{\varepsilon}^m, \mathbf{A}_{\text{sym}}^* \boldsymbol{\varepsilon}^m \rangle \right] \leq \frac{T + \Delta t}{\alpha_{G,\text{LB}} \alpha_{\mathbf{A}_{\text{sym}},\text{LB}}} \sum_{m=0}^{N-1} \| \boldsymbol{\varrho}^m \|_{-1}^2. \quad (A.26) $$

□

## A.3 Proof of Proposition 3.3

*Proof.* From (3.13), we have

$$ \langle \mathbf{M}(\psi_{M,n}^k - \psi_{M,n}^{k+1}) + \Delta t \mathbf{A}^T \psi_{M,n}^k, e^{k+1} \rangle = 0. $$

We then sum over $k = 0, \dots, n-1$, to obtain

$$ \langle \mathbf{M}(\psi_{M,n}^0 - \psi_{M,n}^1), e^1 \rangle + \langle \mathbf{M}(\psi_{M,n}^1 - \psi_{M,n}^2), e^2 \rangle + \dots + \langle \mathbf{M}(\psi_{M,n}^{n-1} - \psi_{M,n}^n), e^n \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{A}^T \psi_{M,n}^k, e^{k+1} \rangle = 0, $$

which gives us

$$ \langle \mathbf{M} \psi_{M,n}^0, e^1 \rangle - \langle \mathbf{M} \psi_{M,n}^1, e^1 \rangle + \langle \mathbf{M} \psi_{M,n}^1, e^2 \rangle - \langle \mathbf{M} \psi_{M,n}^2, e^2 \rangle + \dots \\ + \langle \mathbf{M} \psi_{M,n}^{n-1}, e^n \rangle - \langle \mathbf{M} \psi_{M,n}^n, e^n \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{A}^T \psi_{M,n}^k, e^{k+1} \rangle = 0, $$

leading to

$$ \sum_{k=0}^{n-1} \langle \mathbf{M} \psi_{M,n}^k, e^{k+1} \rangle - \sum_{k=1}^{n-1} \langle \mathbf{M} \psi_{M,n}^k, e^k \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{A}^T \psi_{M,n}^k, e^{k+1} \rangle = \langle \mathbf{M} \psi_{M,n}^n, e^n \rangle. $$

Since $\langle \mathbf{M} \psi_{M,n}^0, e^0 \rangle = 0$, the above equation becomes

$$ \sum_{k=0}^{n-1} \langle \mathbf{M} \psi_{M,n}^k, e^{k+1} - e^k \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{A}^T \psi_{M,n}^k, e^{k+1} \rangle = \langle \mathbf{M} \psi_{M,n}^n, e^n \rangle. $$

Using the final condition of the dual problem (3.13), we can write

$$ \langle \mathbf{M} \psi_{M,n}^n, e^n \rangle = - \langle \mathbf{l}, \mathbf{p}_M^n - \mathbf{p}^{\text{Npr},n} \rangle = \sum_{k=0}^{n-1} [ \langle \mathbf{M} \psi_{M,n}^k, e^{k+1} - e^k \rangle + \Delta t \langle \mathbf{A}^T \psi_{M,n}^k, e^{k+1} \rangle ]. \quad (A.27) $$
On the other side, again using (A.20), we have in particular for $n = 0$

$$ \Delta t \sum_{m=0}^{N-1} \langle \boldsymbol{\epsilon}^m, \mathbf{A}^T \boldsymbol{\epsilon}^m \rangle \le \frac{\Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=0}^{N-1} \| \boldsymbol{\varrho}^m \|_{-1}^2. \quad (A.23) $$

Combining (A.22) and (A.23), leads to

$$ \sum_{m=0}^{N-1} \left[ \langle \boldsymbol{\epsilon}^m, \mathbf{M} \boldsymbol{\epsilon}^m \rangle + \Delta t \langle \boldsymbol{\epsilon}^m, \mathbf{A}^T \boldsymbol{\epsilon}^m \rangle \right] \le \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=0}^{N-1} \| \boldsymbol{\varrho}^m \|_{-1}^2. \quad (A.24) $$

Now, since

$$ \langle \boldsymbol{\epsilon}^m, \mathbf{A}^T \boldsymbol{\epsilon}^m \rangle = \langle \mathbf{A} \boldsymbol{\epsilon}^m, \boldsymbol{\epsilon}^m \rangle = \langle \mathbf{A}_{\text{sym}} \boldsymbol{\epsilon}^m, \boldsymbol{\epsilon}^m \rangle, $$

we get

$$ \sum_{m=0}^{N-1} \left[ \langle \boldsymbol{\epsilon}^m, \mathbf{M} \boldsymbol{\epsilon}^m \rangle + \Delta t \langle \boldsymbol{\epsilon}^m, \mathbf{A}_{\text{sym}} \boldsymbol{\epsilon}^m \rangle \right] \le \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{m=0}^{N-1} \| \boldsymbol{\varrho}^m \|_{-1}^2. \quad (A.25) $$

Finally, we recall the definition of $\alpha_G$, which results in

$$ \sum_{m=0}^{N-1} \left[ \langle \boldsymbol{\epsilon}^m, \mathbf{M} \boldsymbol{\epsilon}^m \rangle + \Delta t \langle \boldsymbol{\epsilon}^m, \mathbf{A}_{\text{sym}}^* \boldsymbol{\epsilon}^m \rangle \right] \le \frac{T + \Delta t}{\alpha_{G,\text{LB}} \alpha_{\mathbf{A}_{\text{sym}},\text{LB}}} \sum_{m=0}^{N-1} \| \boldsymbol{\varrho}^m \|_{-1}^2. \quad (A.26) $$

□

## A.3 Proof of Proposition 3.3

*Proof.* From (3.13), we have

$$ \langle \mathbf{M}(\psi_{M,n}^k - \psi_{M,n}^{k+1}) + \Delta t \mathbf{A}^T \psi_{M,n}^k, \mathbf{e}^{k+1} \rangle = 0. $$

We then sum over $k = 0, \dots, n-1$, to obtain

$$ \langle \mathbf{M}(\psi_{M,n}^0 - \psi_{M,n}^1), \mathbf{e}^1 \rangle + \langle \mathbf{M}(\psi_{M,n}^1 - \psi_{M,n}^2), \mathbf{e}^2 \rangle + \dots + \langle \mathbf{M}(\psi_{M,n}^{n-1} - \psi_{M,n}^n), \mathbf{e}^n \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{A}^T \psi_{M,n}^k, \mathbf{e}^{k+1} \rangle = 0, $$

which gives us

$$ \langle \mathbf{M} \psi_{M,n}^0, \mathbf{e}^1 \rangle - \langle \mathbf{M} \psi_{M,n}^1, \mathbf{e}^1 \rangle + \langle \mathbf{M} \psi_{M,n}^1, \mathbf{e}^2 \rangle - \langle \mathbf{M} \psi_{M,n}^2, \mathbf{e}^2 \rangle + \dots \\ + \langle \mathbf{M} \psi_{M,n}^{n-1}, \mathbf{e}^n \rangle - \langle \mathbf{M} \psi_{M,n}^n, \mathbf{e}^n \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{A}^T \psi_{M,n}^k, \mathbf{e}^{k+1} \rangle = 0, $$

leading to

$$ \sum_{k=0}^{n-1} \langle \mathbf{M} \psi_{M,n}^k, \mathbf{e}^{k+1} \rangle - \sum_{k=1}^{n-1} \langle \mathbf{M} \psi_{M,n}^k, \mathbf{e}^k \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{A}^T \psi_{M,n}^k, \mathbf{e}^{k+1} \rangle = \langle \mathbf{M} \psi_{M,n}^n, \mathbf{e}^n \rangle. $$

Since $\langle \mathbf{M} \psi_{M,n}^0, \mathbf{e}^0 \rangle = 0$, the above equation becomes

$$ \sum_{k=0}^{n-1} \langle \mathbf{M} \psi_{M,n}^k, \mathbf{e}^{k+1} - \mathbf{e}^k \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{A}^T \psi_{M,n}^k, \mathbf{e}^{k+1} \rangle = \langle \mathbf{M} \psi_{M,n}^n, \mathbf{e}^n \rangle. $$

Using the final condition of the dual problem (3.13), we can write

$$ \langle \mathbf{M} \psi_{M,n}^n, \mathbf{e}^n \rangle = - \langle \mathbf{l}, \mathbf{p}_M^n - \mathbf{p}^{\text{Npr},n} \rangle = \sum_{k=0}^{n-1} \left[ \langle \mathbf{M} \psi_{M,n}^k, \mathbf{e}^{k+1} - \mathbf{e}^k \rangle + \Delta t \langle \mathbf{A}^T \psi_{M,n}^k, \mathbf{e}^{k+1} \rangle \right]. \quad (A.27) $$
Equation (A.27) can be rewritten as

$$
\begin{align}
\langle l, p_{\mathcal{M}}^n - p^{\mathcal{N}_{\text{pr},n}} \rangle &= - \sum_{k=0}^{n-1} \langle (M + \Delta t A) p_{\mathcal{M}}^{k+1}, \psi_{\mathcal{M},n}^k \rangle + \sum_{k=0}^{n-1} \langle (M + \Delta t A) p^{\mathcal{N}_{\text{pr},k+1}}, \psi_{\mathcal{M},n}^k \rangle \nonumber \\
&\quad + \sum_{k=0}^{n-1} \langle M p_{\mathcal{M}}^k, \psi_{\mathcal{M},n}^k \rangle - \sum_{k=0}^{n-1} \langle M p^{\mathcal{N}_{\text{pr},k}}, \psi_{\mathcal{M},n}^k \rangle \nonumber \\
&= -\Delta t \sum_{k=0}^{n-1} \langle b, \psi_{\mathcal{M},n}^k \rangle - \sum_{k=0}^{n-1} \langle M p^{\mathcal{N}_{\text{pr},k}}, \psi_{\mathcal{M},n}^k \rangle + \sum_{k=0}^{n-1} \langle (M + \Delta t A) p^{\mathcal{N}_{\text{pr},k+1}}, \psi_{\mathcal{M},n}^k \rangle \nonumber \\
&\quad + \sum_{k=0}^{n-1} \langle (M + \Delta t A) p^{\mathcal{N}_{\text{pr},k+1}}, \psi_{n}^{\mathcal{N}_{\text{du},k}} \rangle - \sum_{k=0}^{n-1} \langle (M + \Delta t A) p^{\mathcal{N}_{\text{pr},k+1}}, \psi_{n}^{\mathcal{N}_{\text{du},k}} \rangle \nonumber \\
&\quad + \sum_{k=0}^{n-1} \langle M p^{\mathcal{N}_{\text{pr},k}}, \psi_{n}^{\mathcal{N}_{\text{du},k}} \rangle - \sum_{k=0}^{n-1} \langle M p^{\mathcal{N}_{\text{pr},k}}, \psi_{n}^{\mathcal{N}_{\text{du},k}} \rangle + \Delta t \sum_{k=0}^{n-1} \langle b, \psi_{n}^{\mathcal{N}_{\text{du},k}} \rangle - \Delta t \sum_{k=0}^{n-1} \langle b, \psi_{n}^{\mathcal{N}_{\text{du},k}} \rangle \nonumber \\
&= \Delta t \sum_{k=0}^{n-1} \langle r^{k+1}, \psi_{\mathcal{M},n}^k - \psi_{n}^{\mathcal{N}_{\text{du},k}} \rangle + \Delta t \sum_{k=0}^{n-1} \langle r^{k+1}, \psi_{n}^{\mathcal{N}_{\text{du},k}} \rangle \nonumber \\
&= \Delta t \sum_{k=0}^{n-1} \langle r^{k+1}, \Psi_{\mathcal{M}}^{N-n+k} - \Psi^{N-n+k} \rangle + \Delta t \sum_{k=0}^{n-1} \langle r^{k+1}, \Psi^{N-n+k} \rangle \nonumber \\
&= \Delta t \sum_{k=0}^{n-1} \langle r^{k+1}, \varepsilon^{N-n+k} \rangle + \Delta t \sum_{k=0}^{n-1} \langle r^{k+1}, \Psi^{N-n+k} \rangle. \tag{A.28}
\end{align}
$$

• First choice. The error bound is evaluated according to

$$
|s^n - s^{\mathcal{N}_s, n}| = \sum_{k=0}^{n-1} \Delta t |\langle \mathbf{r}^{k+1}, \varepsilon^{N-n+k} \rangle|. \qquad (A.29)
$$

We use (3.6) and Cauchy Schwarz inequality to obtain

$$
|s^n - s^{\mathcal{N}_s, n}| \le \left( \sum_{k=0}^{n-1} \Delta t \|\mathbf{r}^{k+1}\|_{-1}^2 \right)^{1/2} \left( \sum_{k=0}^{n-1} \Delta t \|\varepsilon^{N-n+k}\|_{\mathcal{G}^*}^2 \right)^{1/2}. \quad (A.30)
$$

Inequality (A.30) is valid for $n = N$. Hence, we can write

$$
|s^N - s^{\mathcal{N}_s, N}| \le \left( \sum_{k=0}^{N-1} \Delta t \|\mathbf{r}^{k+1}\|_{-1}^2 \right)^{1/2} \left( \sum_{k=0}^{N-1} \Delta t \|\boldsymbol{\varepsilon}^k\|_{\mathcal{G}^*}^2 \right)^{1/2}. \quad (A.31)
$$

We have from the definition of $\alpha_G$ that

$$
\alpha_G \|\boldsymbol{\varepsilon}^k\|_{\mathcal{G}^*}^2 \le \langle (M + \Delta t \mathbf{A}_{\text{sym}}) \boldsymbol{\varepsilon}^k, \boldsymbol{\varepsilon}^k \rangle. \quad (A.32)
$$

We then sum over $\{0, \dots, N-1\}$ and use (A.24) to obtain

$$
\alpha_G \sum_{k=0}^{N-1} \|\boldsymbol{\varepsilon}^k\|_{\mathcal{G}^*}^2 \le \sum_{k=0}^{N-1} \langle (M + \Delta t \mathbf{A}_{\text{sym}}) \boldsymbol{\varepsilon}^k, \boldsymbol{\varepsilon}^k \rangle \le \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{k=0}^{N-1} \|\boldsymbol{\varrho}^k\|_{-1}^2. \quad (A.33)
$$

Therefore,

$$
|s^N - s^{\mathcal{N}_s, N}| \le \Delta t \left( \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2} \Delta_{\text{du}}^N := \Delta_s^N. \qquad (A.34)
$$
Equation (A.27) can be rewritten as

$$
\begin{align}
\langle l, \mathbf{p}_{\mathcal{M}}^n - \mathbf{p}^{\mathcal{N}_{\text{pr}},n} \rangle &= -\sum_{k=0}^{n-1} \langle (\mathbf{M} + \Delta t \mathbf{A}) \mathbf{p}_{\mathcal{M}}^{k+1}, \psi_{\mathcal{M},n}^k \rangle + \sum_{k=0}^{n-1} \langle (\mathbf{M} + \Delta t \mathbf{A}) \mathbf{p}^{\mathcal{N}_{\text{pr}},k+1}, \psi_{\mathcal{M},n}^k \rangle \nonumber \\
&\quad + \sum_{k=0}^{n-1} \langle \mathbf{M} \mathbf{p}_{\mathcal{M}}^k, \psi_{\mathcal{M},n}^k \rangle - \sum_{k=0}^{n-1} \langle \mathbf{M} \mathbf{p}^{\mathcal{N}_{\text{pr}},k}, \psi_{\mathcal{M},n}^k \rangle \nonumber \\
&= -\Delta t \sum_{k=0}^{n-1} \langle \mathbf{b}, \psi_{\mathcal{M},n}^k \rangle - \sum_{k=0}^{n-1} \langle \mathbf{M} \mathbf{p}^{\mathcal{N}_{\text{pr}},k}, \psi_{\mathcal{M},n}^k \rangle + \sum_{k=0}^{n-1} \langle (\mathbf{M} + \Delta t \mathbf{A}) \mathbf{p}^{\mathcal{N}_{\text{pr}},k+1}, \psi_{\mathcal{M},n}^k \rangle \nonumber \\
&\quad + \sum_{k=0}^{n-1} \langle (\mathbf{M} + \Delta t \mathbf{A}) \mathbf{p}^{\mathcal{N}_{\text{pr}},k+1}, \psi_{n}^{\mathcal{N}_{\text{du}},k} \rangle - \sum_{k=0}^{n-1} \langle (\mathbf{M} + \Delta t \mathbf{A}) \mathbf{p}^{\mathcal{N}_{\text{pr}},k+1}, \psi_{n}^{\mathcal{N}_{\text{du}},k} \rangle \nonumber \\
&\quad + \sum_{k=0}^{n-1} \langle \mathbf{M} \mathbf{p}^{\mathcal{N}_{\text{pr}},k}, \psi_{n}^{\mathcal{N}_{\text{du}},k} \rangle - \sum_{k=0}^{n-1} \langle \mathbf{M} \mathbf{p}^{\mathcal{N}_{\text{pr}},k}, \psi_{n}^{\mathcal{N}_{\text{du}},k} \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{b}, \psi_{n}^{\mathcal{N}_{\text{du}},k} \rangle - \Delta t \sum_{k=0}^{n-1} \langle \mathbf{b}, \psi_{n}^{\mathcal{N}_{\text{du}},k} \rangle \nonumber \\
&= \Delta t \sum_{k=0}^{n-1} \langle \mathbf{r}^{k+1}, \psi_{\mathcal{M},n}^k - \psi_{n}^{\mathcal{N}_{\text{du}},k} \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{r}^{k+1}, \psi_{n}^{\mathcal{N}_{\text{du}},k} \rangle \nonumber \\
&= \Delta t \sum_{k=0}^{n-1} \langle \mathbf{r}^{k+1}, \Psi_{\mathcal{M}}^{N-n+k} - \Psi^{N_{\text{du}},N-n+k} \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{r}^{k+1}, \Psi^{N_{\text{du}},N-n+k} \rangle \nonumber \\
&= \Delta t \sum_{k=0}^{n-1} \langle \mathbf{r}^{k+1}, \varepsilon^{N-n+k} \rangle + \Delta t \sum_{k=0}^{n-1} \langle \mathbf{r}^{k+1}, \Psi^{N_{\text{du}},N-n+k} \rangle. \tag{A.28}
\end{align}
$$

• First choice. The error bound is evaluated according to

$$
|s^n - s^{\mathcal{N}_s, n}| = \sum_{k=0}^{n-1} \Delta t |\langle \mathbf{r}^{k+1}, \varepsilon^{N-n+k} \rangle|. \qquad (A.29)
$$

We use (3.6) and Cauchy Schwarz inequality to obtain

$$
|s^n - s^{\mathcal{N}_s, n}| \le \left( \sum_{k=0}^{n-1} \Delta t \|\mathbf{r}^{k+1}\|_{-1}^2 \right)^{1/2} \left( \sum_{k=0}^{n-1} \Delta t \|\varepsilon^{N-n+k}\|_{\mathcal{G}^*}^2 \right)^{1/2}. \quad (A.30)
$$

Inequality (A.30) is valid for $n = N$. Hence, we can write

$$
|s^N - s^{\mathcal{N}_s, N}| \le \left( \sum_{k=0}^{N-1} \Delta t \|\mathbf{r}^{k+1}\|_{-1}^2 \right)^{1/2} \left( \sum_{k=0}^{N-1} \Delta t \|\boldsymbol{\varepsilon}^k\|_{\mathcal{G}^*}^2 \right)^{1/2}. \quad (A.31)
$$

We have from the definition of $\alpha_G$ that

$$
\alpha_G \|\boldsymbol{\varepsilon}^k\|_{\mathcal{G}^*}^2 \le \langle (\mathbf{M} + \Delta t \mathbf{A}_{\text{sym}}) \boldsymbol{\varepsilon}^k, \boldsymbol{\varepsilon}^k \rangle. \quad (A.32)
$$

We then sum over $\{0, \dots, N-1\}$ and use (A.24) to obtain

$$
\alpha_G \sum_{k=0}^{N-1} \|\boldsymbol{\varepsilon}^k\|_{\mathcal{G}^*}^2 \le \sum_{k=0}^{N-1} \langle (\mathbf{M} + \Delta t \mathbf{A}_{\text{sym}}) \boldsymbol{\varepsilon}^k, \boldsymbol{\varepsilon}^k \rangle \le \frac{T + \Delta t}{\alpha_{\mathbf{A}_{\text{sym}}}} \sum_{k=0}^{N-1} \|\boldsymbol{\varrho}^k\|_{-1}^2. \quad (A.33)
$$

Therefore,

$$
|s^N - s^{\mathcal{N}_s, N}| \le \Delta t \left( \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2} \Delta_{\text{du}}^N := \Delta_s^N. \qquad (A.34)
$$
• Second choice. The error bound at $n = N$ is evaluated according to

$$ \begin{aligned} |s^N - \hat{s}^{N_s, N}| &\le \Delta t \sum_{n=0}^{N-1} |\langle \mathbf{r}^{n+1}, \varepsilon^{N-N+n} \rangle| + \Delta t \sum_{n=0}^{N-1} |\langle \mathbf{r}^{n+1}, \mathbf{\Psi}^{N_{\text{du}}, N-N+n} \rangle| \\ &\le \Delta t \left( \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2} \Delta_{\text{du}}^N + \Delta t \sum_{n=0}^{N-1} |\langle \mathbf{r}^{n+1}, \mathbf{\Psi}^{N_{\text{du}}, n} \rangle| =: \tilde{\Delta}_s^N. \end{aligned} \quad (A.35) $$

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

## B Successive constraint method

Let $\Xi$ be a set of parameter values. For each $\xi \in \Xi$, the successive constraint method (SCM) consists in finding an upper bound $\alpha_{\text{UB}}(\xi)$ and a lower bound $\alpha_{\text{LB}}(\xi)$ of the coercivity constant $\alpha(\xi)$ through an offline-online strategy. The SCM relies on the affine decomposition assumption (3.25), which enables us to express $\alpha(\xi)$ as

$$ \alpha(\xi) = \inf_{\mathbf{v} \in \mathbb{R}^N} \sum_{d=1}^{D_a} \Theta_d^a(\xi) \frac{\mathbf{v}^T \mathbf{A}_d \mathbf{v}}{\|\mathbf{v}\|_{G^*}^2} = \inf_{\mathbf{v} \in \mathbb{R}^N} \sum_{d=1}^{D_a} \Theta_d^a(\xi) w_d. \quad (B.1) $$

To define the lower bound $\alpha_{\text{LB}}(\xi)$, we express (B.1) as a minimization problem

$$ \alpha(\xi) = \inf_{\mathbf{w} \in \mathcal{W}} \mathcal{J}(\xi, \mathbf{w}), \quad (B.2) $$

where the set $\mathcal{W}$ is defined as

$$ \mathcal{W} := \left\{ \mathbf{w} = (w_1, \dots, w_{D_a}) \in \mathbb{R}^{D_a} \mid \exists \mathbf{v} \in \mathbb{R}^N \text{ s.t. } w_d = \frac{\mathbf{v}^T \mathbf{A}_d \mathbf{v}}{\|\mathbf{v}\|_{G^*}^2}, 1 \le d \le D_a \right\}, $$

and the objective function is given by

$$ \mathcal{J}: \Xi \times \mathbb{R}^{D_a} \to \mathbb{R} $$

$$ (\xi, \mathbf{w}) \mapsto \mathcal{J}(\xi, \mathbf{w}) = \sum_{d=1}^{D_a} \Theta_d^a(\xi) w_d. $$

The idea of the SCM is based on creating two sets $\mathcal{W}_{\text{LB}}$ and $\mathcal{W}_{\text{UB}}$, such that $\mathcal{W}_{\text{UB}} \subset \mathcal{W} \subset \mathcal{W}_{\text{LB}}$, where we perform the minimization over these two sets and define

$$ \alpha_{\text{LB}}(\xi) = \min_{\mathbf{w} \in \mathcal{W}_{\text{LB}}} \mathcal{J}(\xi, \mathbf{w}) \quad \text{and} \quad \alpha_{\text{UB}}(\xi) = \min_{\mathbf{w} \in \mathcal{W}_{\text{UB}}} \mathcal{J}(\xi, \mathbf{w}). $$

**Definition of $\mathcal{W}_{\text{UB}}$.** We introduce the subset of parameter values $\Xi_M \subset \Xi$ obtained using a greedy algorithm (see Algorithm 5). The construction of $\Xi_M$ requires a training set $\Xi_{\text{training}}$ and a fixed tolerance $0 \le \text{tol} \le 1$ that controls the relative gap between the lower and upper bounds.

For all $1 \le j \le M$ and for each $\xi_j \in \Xi_M$,

1. we assemble $\mathbf{A}(\xi_j) = \sum_{d=1}^{D_a} \Theta_d^a(\xi_j) \mathbf{A}_d$,

2. we solve the generalized eigenvalue problem

$$ \mathbf{A}(\xi_j)\mathbf{y} = \lambda \mathbf{G}^*(\xi^*)\mathbf{y}, \quad (B.3) $$

and extract the smallest eigenvalue $\alpha^j$ and its corresponding eigenvector $\mathbf{v}^j$,

3. we compute the vector $\mathbf{w}^j \in \mathbb{R}^{D_a}$ such that

$$ (\mathbf{w}^j)_d = \frac{(\mathbf{v}^j)^T \mathbf{A}_d \mathbf{v}^j}{\|\mathbf{v}^j\|_{G^*}^2} $$

4. we define the set

$$ \mathcal{W}_{\text{UB}} = \{\mathbf{w}^j \mid 1 \le j \le M\}, $$

and compute the upper bound

$$ \alpha_{\text{UB}}(\xi) = \arg \min_{\mathbf{w} \in \mathcal{W}_{\text{UB}}} \mathcal{J}(\xi, \mathbf{w}). $$
• Second choice. The error bound at $n = N$ is evaluated according to

$$ \begin{aligned} |s^N - \tilde{s}^{N_s, N}| &\le \Delta t \sum_{n=0}^{N-1} |\langle \mathbf{r}^{n+1}, \varepsilon^{N-N+n} \rangle| + \Delta t \sum_{n=0}^{N-1} |\langle \mathbf{r}^{n+1}, \mathbf{\Psi}^{N_{\text{du}}, N-N+n} \rangle| \\ &\le \Delta t \left( \sum_{n=1}^{N} \|\mathbf{r}^n\|_{-1}^2 \right)^{1/2} \Delta_{\text{du}}^N + \Delta t \sum_{n=0}^{N-1} |\langle \mathbf{r}^{n+1}, \mathbf{\Psi}^{N_{\text{du}}, n} \rangle| =: \tilde{\Delta}_s^N. \end{aligned} \quad (A.35) $$

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

## B Successive constraint method

Let $\Xi$ be a set of parameter values. For each $\xi \in \Xi$, the successive constraint method (SCM) consists in finding an upper bound $\alpha_{\text{UB}}(\xi)$ and a lower bound $\alpha_{\text{LB}}(\xi)$ of the coercivity constant $\alpha(\xi)$ through an offline-online strategy. The SCM relies on the affine decomposition assumption (3.25), which enables us to express $\alpha(\xi)$ as

$$ \alpha(\xi) = \inf_{\mathbf{v} \in \mathbb{R}^N} \sum_{d=1}^{D_a} \Theta_d^a(\xi) \frac{\mathbf{v}^T \mathbf{A}_d \mathbf{v}}{\|\mathbf{v}\|_{G^*}^2} = \inf_{\mathbf{v} \in \mathbb{R}^N} \sum_{d=1}^{D_a} \Theta_d^a(\xi) w_d. \quad (B.1) $$

To define the lower bound $\alpha_{\text{LB}}(\xi)$, we express (B.1) as a minimization problem

$$ \alpha(\xi) = \inf_{\mathbf{w} \in \mathcal{W}} \mathcal{J}(\xi, \mathbf{w}), \quad (B.2) $$

where the set $\mathcal{W}$ is defined as

$$ \mathcal{W} := \left\{ \mathbf{w} = (w_1, \dots, w_{D_a}) \in \mathbb{R}^{D_a} \mid \exists \mathbf{v} \in \mathbb{R}^N \text{ s.t. } w_d = \frac{\mathbf{v}^T \mathbf{A}_d \mathbf{v}}{\|\mathbf{v}\|_{G^*}^2}, 1 \le d \le D_a \right\}, $$

and the objective function is given by

$$ \mathcal{J}: \Xi \times \mathbb{R}^{D_a} \to \mathbb{R} $$

$$ (\xi, \mathbf{w}) \mapsto \mathcal{J}(\xi, \mathbf{w}) = \sum_{d=1}^{D_a} \Theta_d^a(\xi) w_d. $$

The idea of the SCM is based on creating two sets $\mathcal{W}_{\text{LB}}$ and $\mathcal{W}_{\text{UB}}$, such that $\mathcal{W}_{\text{UB}} \subset \mathcal{W} \subset \mathcal{W}_{\text{LB}}$, where we perform the minimization over these two sets and define

$$ \alpha_{\text{LB}}(\xi) = \min_{\mathbf{w} \in \mathcal{W}_{\text{LB}}} \mathcal{J}(\xi, \mathbf{w}) \quad \text{and} \quad \alpha_{\text{UB}}(\xi) = \min_{\mathbf{w} \in \mathcal{W}_{\text{UB}}} \mathcal{J}(\xi, \mathbf{w}). $$

**Definition of $\mathcal{W}_{\text{UB}}$.** We introduce the subset of parameter values $\Xi_M \subset \Xi$ obtained using a greedy algorithm (see Algorithm 5). The construction of $\Xi_M$ requires a training set $\Xi_{\text{training}}$ and a fixed tolerance $0 \le \text{tol} \le 1$ that controls the relative gap between the lower and upper bounds.

For all $1 \le j \le M$ and for each $\xi_j \in \Xi_M$,

1. we assemble $\mathbf{A}(\xi_j) = \sum_{d=1}^{D_a} \Theta_d^a(\xi_j) \mathbf{A}_d$,

2. we solve the generalized eigenvalue problem

$$ \mathbf{A}(\xi_j)\mathbf{y} = \lambda\mathbf{G}^*(\xi^*)\mathbf{y}, \quad (B.3) $$

and extract the smallest eigenvalue $\alpha^j$ and its corresponding eigenvector $\mathbf{v}^j$,

3. we compute the vector $\mathbf{w}^j \in \mathbb{R}^{D_a}$ such that

$$ (\mathbf{w}^j)_d = \frac{(\mathbf{v}^j)^T \mathbf{A}_d \mathbf{v}^j}{\|\mathbf{v}^j\|_{G^*}^2} $$

4. we define the set

$$ \mathcal{W}_{\text{UB}} = \{\mathbf{w}^j \mid 1 \le j \le M\}, $$

and compute the upper bound

$$ \alpha_{\text{UB}}(\xi) = \arg \min_{\mathbf{w} \in \mathcal{W}_{\text{UB}}} \mathcal{J}(\xi, \mathbf{w}). $$
**Algorithm 5 Construction of Ξ<sub>M</sub>**

**Input:** tol, Ξ.

1. Choose arbitrary $\xi_1 \in \Xi$.

2. Set $j = 1$ and $\Xi_j = \{\xi_1\}$.

3. Compute $\eta_j(\xi) = \frac{\alpha_{UB}(\xi) - \alpha_{LB}(\xi)}{\alpha_{UB}(\xi)}$.

4. **while** $\max_{\xi \in \Xi} \eta_j(\xi) > \text{tol}$ **do**

5. Compute $\xi_{j+1} = \arg \max_{\xi \in \Xi} \eta_j(\xi)$.

6. Set $\Xi_{j+1} = \Xi_j \cup \{\xi_{j+1}\}$.

7. $j \leftarrow j + 1$.

8. $\eta_j(\xi) = \frac{\alpha_{UB}(\xi) - \alpha_{LB}(\xi)}{\alpha_{UB}(\xi)}$.

9. **end while**

**Definition of W<sub>LB</sub>.** First, we need to introduce the constraint interval

$$ \mathcal{B} = \prod_{d=1}^{D_a} \left[ \inf_{\mathbf{v} \in \mathbb{R}^N} \frac{\mathbf{v}^T \mathbf{A}_d \mathbf{v}}{\|\mathbf{v}\|_{G^*}^2}, \sup_{\mathbf{v} \in \mathbb{R}^N} \frac{\mathbf{v}^T \mathbf{A}_d \mathbf{v}}{\|\mathbf{v}\|_{G^*}^2} \right], $$

obtained by computing, once at the beginning of the SCM, the smallest and largest eigenvalues of a problem similar to (B.3) and obtained by replacing $\mathbf{A}(\xi_j)$ by $\mathbf{A}_d$. We define the set

$$ \mathcal{W}_{\text{LB}}^{j}(\xi) = \left\{ \mathbf{w} \in \mathcal{B} \mid \begin{array}{l} \mathcal{J}(\xi'; \mathbf{w}) \ge \alpha(\xi'), \quad \forall \xi' \in P_{M_1}(\xi; \Xi_j); \\ \mathcal{J}(\xi'; \mathbf{w}) \ge \alpha_{\text{LB}}^{j-1}(\xi'), \quad \forall \xi' \in P_{M_2}(\xi; \Xi \setminus \Xi_j) \end{array} \right\}, $$

where $P_M(\xi; \mathbb{D}) := \{\text{M closest points to } \xi \text{ in } \mathbb{D}\}$.

## C Empirical interpolation method

The efficiency of the RB method relies on the affine decomposition (3.25) proposed in Section 3.5. However, this decomposition is not always available. But the empirical interpolation method (EIM) can provide one to approximate, in our case, $\hat{v}(\xi)$ with an affine sum. Given a family of parameter-dependent vectors $\mathcal{T} = \{\hat{v}(\xi) \in \mathbb{R}^F; \xi \in \Xi_{\text{training}}\}$, the EIM aims at finding an approximation to the elements of $\mathcal{T}$ through an operator $\mathcal{I}_{M_{\text{EIM}}}$ that interpolates the vector $\hat{v}(\xi)$ at some selected points. Using a greedy process, we construct the set of vectors $\{\tilde{v}^1, \dots, \tilde{v}^{M_{\text{EIM}}}\}$ and the interpolation points $\{x_1, \dots, x_{M_{\text{EIM}}}\}$ such that

$$ \mathcal{I}_{M_{\text{EIM}}}[\hat{v}(\xi)] \approx \sum_{i=1}^{M_{\text{EIM}}} \Theta_i(\xi) \tilde{v}^i, \qquad (C.1a) $$

where $\Theta_d(\xi) \in \mathbb{R}$ and $\tilde{v}^i \in \mathbb{R}^F$, $1 \le i \le M_{\text{EIM}}$, do not depend on $\xi$.

To begin the procedure, we randomly choose $\xi_1$ from $\Xi_{\text{training}}$ and set $\hat{v}^1 = \hat{v}(\xi_1)$. The first interpolation point is chosen such that

$$ x_1 = \arg \max_{1 \le j \le F} |\hat{v}_j^1|, $$

where $\hat{v}_j^1$ is the j-th element of $\hat{v}^1$. We then initialize the first basis function as

$$ \tilde{v}^1 = \hat{v}^1 / \hat{v}_{j_1}^1, $$

with $1 \le j_1 \le F$, the index corresponding to the selected point $x_1$. At the m-th step, $m = 2, \dots, M_{\text{EIM}}-1$, given the set of interpolations points $\{x_1, \dots, x_{M_{\text{EIM}}-1}\}$ and the set of basis elements $\{\tilde{v}^1, \dots, \tilde{v}^{M_{\text{EIM}}-1}\}$, we select the next snapshot as the worst approximated one by the current interpolant. To do so, we first write the m equations stating the equality between the current EIM approximation and a vector $\hat{v}(\xi)$ at the current m interpolation points. This
**Algorithm 5 Construction of Ξ<sub>M</sub>**

**Input:** tol, Ξ.

1. Choose arbitrary $\xi_1 \in \Xi$.

2. Set $j = 1$ and $\Xi_j = \{\xi_1\}$.

3. Compute $\eta_j(\xi) = \frac{\alpha_{UB}(\xi)-\alpha_{LB}(\xi)}{\alpha_{UB}(\xi)}$.

4. **while** $\max_{\xi \in \Xi} \eta_j(\xi) > \text{tol}$ **do**

5. Compute $\xi_{j+1} = \arg \max_{\xi \in \Xi} \eta_j(\xi)$.

6. Set $\Xi_{j+1} = \Xi_j \cup \{\xi_{j+1}\}$.

7. $j \leftarrow j + 1$.

8. $\eta_j(\xi) = \frac{\alpha_{UB}(\xi)-\alpha_{LB}(\xi)}{\alpha_{UB}(\xi)}$.

9. **end while**

**Definition of W<sub>LB</sub>.** First, we need to introduce the constraint interval

$$ \mathcal{B} = \prod_{d=1}^{D_a} \left[ \inf_{v \in \mathbb{R}^N} \frac{v^T A_d v}{\|v\|_{G^*}^2}, \sup_{v \in \mathbb{R}^N} \frac{v^T A_d v}{\|v\|_{G^*}^2} \right], $$

obtained by computing, once at the beginning of the SCM, the smallest and largest eigenvalues of a problem similar to (B.3) and obtained by replacing $A(\xi_j)$ by $A_d$. We define the set

$$ W_{\text{LB}}^{j}(\xi) = \left\{ w \in \mathcal{B} \mid \begin{array}{l} \mathcal{J}(\xi'; w) \ge \alpha(\xi'), \quad \forall \xi' \in P_{M_1}(\xi; \Xi_j); \\ \mathcal{J}(\xi'; w) \ge \alpha_{\text{LB}}^{j-1}(\xi'), \quad \forall \xi' \in P_{M_2}(\xi; \Xi \setminus \Xi_j) \end{array} \right\}, $$

where $P_M(\xi; D) := \{\text{M closest points to } \xi \text{ in } D\}$.

## C Empirical interpolation method

The efficiency of the RB method relies on the affine decomposition (3.25) proposed in Section 3.5. However, this decomposition is not always available. But the empirical interpolation method (EIM) can provide one to approximate, in our case, $\hat{v}(\xi)$ with an affine sum. Given a family of parameter-dependent vectors $\mathcal{T} = \{\hat{v}(\xi) \in \mathbb{R}^F; \xi \in \Xi_{\text{training}}\}$, the EIM aims at finding an approximation to the elements of $\mathcal{T}$ through an operator $\mathcal{I}_{M_{\text{EIM}}}$ that interpolates the vector $\hat{v}(\xi)$ at some selected points. Using a greedy process, we construct the set of vectors $\{\tilde{v}^1, \dots, \tilde{v}^{M_{\text{EIM}}}\}$ and the interpolation points $\{x_1, \dots, x_{M_{\text{EIM}}}\}$ such that

$$ \mathcal{I}_{M_{\text{EIM}}}[\hat{v}(\xi)] \approx \sum_{i=1}^{M_{\text{EIM}}} \Theta_i(\xi) \tilde{v}^i, \qquad (C.1a) $$

where $\Theta_d(\xi) \in \mathbb{R}$ and $\tilde{v}^i \in \mathbb{R}^F$, $1 \le i \le M_{\text{EIM}}$, do not depend on $\xi$.

To begin the procedure, we randomly choose $\xi_1$ from $\Xi_{\text{training}}$ and set $\hat{v}^1 = \hat{v}(\xi_1)$. The first interpolation point is chosen such that

$$ x_1 = \arg \max_{1 \le j \le F} |\hat{v}_j^1|, $$

where $\hat{v}_j^1$ is the j-th element of $\hat{v}^1$. We then initialize the first basis function as

$$ \tilde{v}^1 = \hat{v}^1 / \hat{v}_{j_1}^1, $$

with $1 \le j_1 \le F$, the index corresponding to the selected point $x_1$. At the m-th step, $m = 2, \dots, M_{\text{EIM}}-1$, given the set of interpolations points $\{x_1, \dots, x_{M_{\text{EIM}}-1}\}$ and the set of basis elements $\{\tilde{v}^1, \dots, \tilde{v}^{M_{\text{EIM}}-1}\}$, we select the next snapshot as the worst approximated one by the current interpolant. To do so, we first write the m equations stating the equality between the current EIM approximation and a vector $\hat{v}(\xi)$ at the current m interpolation points. This
leads to the following lower triangular linear system:

$$ \begin{bmatrix} 1 & 0 & \cdots & 0 & 0 \\ \tilde{\mathbf{v}}_{j_2}^1 & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & & \vdots & \vdots \\ \tilde{\mathbf{v}}_{j_{M_{\text{EIM}}-1}}^1 & \tilde{\mathbf{v}}_{j_{M_{\text{EIM}}-1}}^2 & \cdots & 1 & 0 \\ \tilde{\mathbf{v}}_{j_{M_{\text{EIM}}}}^1 & \tilde{\mathbf{v}}_{j_{M_{\text{EIM}}}}^2 & \cdots & \tilde{\mathbf{v}}_{j_{M_{\text{EIM}}}}^{M_{\text{EIM}}} & 1 \end{bmatrix} \begin{bmatrix} \Theta_1 \\ \Theta_2 \\ \vdots \\ \Theta_{M_{\text{EIM}}-1} \\ \Theta_{M_{\text{EIM}}} \end{bmatrix} (\xi) = \begin{bmatrix} \hat{\mathbf{v}}_{j_1} \\ \hat{\mathbf{v}}_{j_2} \\ \vdots \\ \hat{\mathbf{v}}_{j_{M_{\text{EIM}}-1}} \\ \hat{\mathbf{v}}_{j_{M_{\text{EIM}}}} \end{bmatrix} (\xi). $$

We choose

$$ \xi_{m+1} = \arg \max_{\xi \in \Xi_{\text{training}}} \| \hat{\mathbf{v}}(\xi) - \mathcal{I}_m[\hat{\mathbf{v}}(\xi)] \|_{L^\infty}. \qquad (C.2a) $$

The $(m+1)$-th interpolation point is then defined as

$$ x_{m+1} = \arg \max_{1 \le j \le F} |r_j^{m+1}| $$

with $r^{m+1} = \hat{\mathbf{v}}(\xi_{m+1}) - \mathcal{I}_m[\hat{\mathbf{v}}(\xi_{m+1})]$ and the corresponding basis vector is taken as

$$ \tilde{\mathbf{v}}^{m+1} = \mathbf{r}^{m+1} / r_{j_{m+1}}^{m+1}. $$

We repeat this procedure until a given tolerance $\epsilon_{\text{EIM}} > 0$ is reached, i.e.

$$ \max_{\xi \in \Xi_{\text{training}}} \| \hat{\mathbf{v}}(\xi) - \mathcal{I}_m[\hat{\mathbf{v}}(\xi)] \|_{L^\infty} < \epsilon_{\text{EIM}}. $$
leads to the following lower triangular linear system:

$$ \begin{bmatrix} 1 & 0 & \cdots & 0 & 0 \\ \tilde{\mathbf{v}}_{j_2}^1 & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \tilde{\mathbf{v}}_{j_{M_{\text{EIM}}-1}}^1 & \tilde{\mathbf{v}}_{j_{M_{\text{EIM}}-1}}^2 & \cdots & 1 & 0 \\ \tilde{\mathbf{v}}_{j_{M_{\text{EIM}}}}^1 & \tilde{\mathbf{v}}_{j_{M_{\text{EIM}}}}^2 & \cdots & \tilde{\mathbf{v}}_{j_{M_{\text{EIM}}}}^{M_{\text{EIM}}} & 1 \end{bmatrix} \begin{bmatrix} \Theta_1 \\ \Theta_2 \\ \vdots \\ \Theta_{M_{\text{EIM}}-1} \\ \Theta_{M_{\text{EIM}}} \end{bmatrix} (\xi) = \begin{bmatrix} \hat{\mathbf{v}}_{j_1} \\ \hat{\mathbf{v}}_{j_2} \\ \vdots \\ \hat{\mathbf{v}}_{j_{M_{\text{EIM}}-1}} \\ \hat{\mathbf{v}}_{j_{M_{\text{EIM}}}} \end{bmatrix} (\xi). $$

We choose

$$ \xi_{m+1} = \arg \max_{\xi \in \Xi_{\text{training}}} \| \hat{\mathbf{v}}(\xi) - \mathcal{I}_m[\hat{\mathbf{v}}(\xi)] \|_{L^\infty}. \quad (C.2a) $$

The $(m+1)$-th interpolation point is then defined as

$$ x_{m+1} = \arg \max_{1 \le j \le F} |r_j^{m+1}| $$

with $r^{m+1} = \hat{\mathbf{v}}(\xi_{m+1}) - \mathcal{I}_m[\hat{\mathbf{v}}(\xi_{m+1})]$ and the corresponding basis vector is taken as

$$ \tilde{\mathbf{v}}^{m+1} = \mathbf{r}^{m+1}/r_{j_{m+1}}^{m+1}. $$

We repeat this procedure until a given tolerance $\epsilon_{\text{EIM}} > 0$ is reached, i.e.

$$ \max_{\xi \in \Xi_{\text{training}}} \| \hat{\mathbf{v}}(\xi) - \mathcal{I}_m[\hat{\mathbf{v}}(\xi)] \|_{L^\infty} < \epsilon_{\text{EIM}}. $$
References

[1] I. AAVATSMARK, T. BARKVE, Ø. BØE, AND T. MANNSETH, Discretization on non-orthogonal, curvilinear grids for multi-phase flow, ECMOR IV - 4th European Conference on the Mathematics of Oil Recovery, (1994), https://doi.org/10.3997/2214-4609.201411179.

[2] L. AGÉLAS, R. EYMARD, AND R. HERBIN, *A nine-point finite volume scheme for the simulation of diffusion in heterogeneous media*, Comptes Rendus Mathématique, 347 (2009), pp. 673-676, https://doi.org/10.1016/j.crma.2009.03.013.

[3] L. BEIRÃO DA VEIGA, F. BREZZI, A. CANGIANI, G. MANZINI, L. D. MARINI, AND A. RUSSO, *Basic principles of virtual element methods*, Mathematical Models and Methods in Applied Sciences, 23 (2013), pp. 199–214, https://doi.org/10.1142/S0218202512500492.

[4] M. J. BERGER AND J. OLIGER, *Adaptive mesh refinement for hyperbolic partial differential equations*, J. Comput. Phys., 53 (1984), pp. 484–512, https://doi.org/10.1016/0021-9991(84)90073-1.

[5] Y. BRENIER AND J. JAFFRÉ, *Upstream differencing for multiphase flow in reservoir simulation*, SIAM Journal on Numerical Analysis, 28 (1991), pp. 685–696, https://doi.org/10.1137/0728036.

[6] F. BREZZI, K. LIPNIKOV, AND V. SIMONCINI, *A family of mimetic finite difference methods on polygonal and polyhedral meshes*, Mathematical Models and Methods in Applied Sciences, 15 (2005), pp. 1533–1551, https://doi.org/10.1142/S0218202505000832.

[7] A. BUHR, C. ENGWER, M. OHLBERGER, AND S. RAVE, *A numerically stable a posteriori error estimator for reduced basis approximations of elliptic equations*, 2014, https://arxiv.org/abs/1407.8005.

[8] C. CANUTO, T. TONN, AND K. URBAN, *A posteriori error analysis of the reduced basis method for nonaffine parametrized nonlinear pdes*, SIAM Journal on Numerical Analysis, 47 (2009), pp. 2001–2022, https://doi.org/10.1137/080724812.

[9] F. CASENAVE, *Accurate a posteriori error evaluation in the reduced basis method*, Comptes Rendus Mathématique, 350 (2012), pp. 539–542, https://doi.org/10.1016/j.crma.2012.05.012.

[10] S. CHATURANTABUT AND D. C. SORENSEN, *Nonlinear model reduction via discrete empirical interpolation*, SIAM Journal of Scientific Computing, 32 (2010), pp. 2737–2764, https://doi.org/10.1137/090766498.

[11] G. CHAVENT AND J. JAFFRÉ, *Mathematical models and finite elements for reservoir simulation: single phase, multiphase and multicomponent flows through porous media*, vol. 17, Elsevier, 1986.

[12] Y. CHEN, J. S. HESTHAVEN, Y. MADAY, AND J. RODRÍGUEZ, *A monotonic evaluation of lower bounds for inf-sup stability constants in the frame of reduced basis approximations*, Comptes Rendus Mathématique, 346 (2008), pp. 1295–1300, https://doi.org/10.1016/j.crma.2008.10.012.

[13] Z. CHEN, *Reservoir simulation: mathematical techniques in oil recovery*, vol. 77 of CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, 2007.

[14] Z. CHEN AND R. E. EWING, *Degenerate two-phase incompressible flow iii*, Numerische Mathematik, 90 (2001), pp. 215–240, https://doi.org/10.1007/s002110100291.

[15] R. EYMARD, T. GALLOUËT, AND R. HERBIN, *Finite volume methods*, Handbook of numerical analysis, 7 (2000), pp. 713–7108, https://doi.org/10.1016/S1570-8659(00)07005-8.

[16] R. EYMARD, T. GALLOUËT, AND R. HERBIN, *Discretisation of heterogeneous and anisotropic diffusion problems on general nonconforming meshes*, *sushi: a scheme using stabilisation and hybrid interfaces*, IMA Journal of Numerical Analysis, 30 (2009), pp. 1009–1043, https://doi.org/10.1093/imanum/drn084.

[17] R. EYMARD, C. GUICHARD, R. HERBIN, AND R. MASSON, *Vertex-centred discretization of multiphase compositional darcy flows on general meshes*, Computational Geosciences, 16 (2012), pp. 987–1005, https://doi.org/10.1007/s10596-012-9299-x.

[18] M. A. GREPL, *Reduced basis approximation and a posteriori error estimation for parabolic partial differential equations*, PhD thesis, MIT, 2005, http://hdl.handle.net/1721.1/32387.

[19] M. A. GREPL, Y. MADAY, N. C. NGUYEN, AND A. T. PATERA, *Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations*, ESAIM: Mathematical Modelling and Numerical Analysis, 41 (2007), pp. 575–605, https://doi.org/DOI:10.1051/m2an:2007031.

[20] M. A. GREPL AND A. T. PATERA, *A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations*, ESAIM: Mathematical Modelling and Numerical Analysis, 39 (2005), pp. 157–181, https://doi.org/10.1051/m2an:2005006.

[21] B. HAASDONK AND M. OHLBERGER, *Reduced basis method for finite volume approximations of parametrized linear evolution equations*, ESAIM: Mathematical Modelling and Numerical Analysis, 42 (2008), pp. 277–302, https://doi.org/10.1051/m2an:2008001.