# Adaptive space-time model order reduction with dual-weighted residual (MORé DWR) error control for poroelasticity

Hendrik Fischer, Julian Roth, Ludovic Chamoin, Amélie Fau, Mary Wheeler, Thomas Wick

► To cite this version:

Hendrik Fischer, Julian Roth, Ludovic Chamoin, Amélie Fau, Mary Wheeler, et al., Adaptive space-time model order reduction with dual-weighted residual (MORé DWR) error control for poroelasticity. *Advanced Modeling and Simulation in Engineering Sciences*, 2024, 11 (1), pp.9. 10.1186/s40323-024-00262-6. [hal-04552013](https://hal.science/hal-04552013)

HAL Id: hal-04552013

https://hal.science/hal-04552013

Submitted on 18 Apr 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

Open Access

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

Check for
updates

# Adaptive space-time model order reduction with dual-weighted residual (MORe DWR) error control for poroelasticity

Hendrik Fischer<fnref n="1" />,<fnref n="2" />, Julian Roth<fnref n="1" />,<fnref n="2" />, Ludovic Chamoin<fnref n="2" />, Amélie Fau<fnref n="2" />, Mary Wheeler<fnref n="3" /> and Thomas Wick<fnref n="1" />,<fnref n="2" />*

## Abstract

In this work, the space-time MORe DWR (Model Order Reduction with Dual-Weighted Residual error estimates) framework is extended and further developed for single-phase flow problems in porous media. Specifically, our problem statement is the Biot system which consists of vector-valued displacements (geomechanics) coupled to a Darcy flow pressure equation. The MORe DWR method introduces a goal-oriented adaptive incremental proper orthogonal decomposition (POD) based-reduced-order model (ROM). The error in the reduced goal functional is estimated during the simulation, and the POD basis is enriched on-the-fly if the estimate exceeds a given threshold. This results in a reduction of the total number of full-order-model solves for the simulation of the porous medium, a robust estimation of the quantity of interest and well-suited reduced bases for the problem at hand. We apply a space-time Galerkin discretization with Taylor-Hood elements in space and a discontinuous Galerkin method with piecewise constant functions in time. The latter is well-known to be similar to the backward Euler scheme. We demonstrate the efficiency of our method on the well-known two-dimensional Mandel benchmark and a three-dimensional footing problem.

## Introduction

Porous media problems have long-standing applications in subsurface modeling, ground-water flow, hydraulic fracturing, geothermal energy recovery, and nuclear waste storage. The resulting mathematical models yield coupled systems of partial differential equations (PDEs) from which one of the most well-known is the so-called Biot problem [18–21,71]. Developments over the last two decades have led to more complicated models, due to additional physics that can be incorporated. Nonstationary, nonlinear coupled PDE systems are often obtained. They are also sometimes subject to inequality constraints as in multiphysics phase-field fracture in porous media [75,78]. Despite advances in numerical solvers of iterative or multigrid type, high-performance parallel

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

<footnote>*Correspondence:
thomas.wick@ifam.uni-
hannover.de</footnote>

<footnote n="1">¹Leibniz Universität Hannover,
Institut für Angewandte
Mathematik, AG
Wissenschaftliches Rechnen,
Welfengarten 1, 30167
Hannover, Germany</footnote>

<footnote n="2">²Université Paris-Saclay,
CentraleSupélec, ENS
Paris-Saclay, CNRS, LMPS,
Laboratoire de Mécanique
Paris-Saclay, 91190
Gif-sur-Yvette, France</footnote>

<footnote n="3">³The University of Texas at
Austin, Oden Institute, Austin, TX
78712, USA</footnote>
RESEARCH ARTICLE

Open Access

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

Check for
updates

# Adaptive space-time model order reduction with dual-weighted residual (MORe DWR) error control for poroelasticity

Hendrik Fischer<fnref n="1" />,<fnref n="2" />, Julian Roth<fnref n="1" />,<fnref n="2" />, Ludovic Chamoin<fnref n="2" />, Amélie Fau<fnref n="2" />, Mary Wheeler<fnref n="3" /> and Thomas Wick<fnref n="1" />,<fnref n="2" />*

## Abstract

In this work, the space-time MORe DWR (Model Order Reduction with Dual-Weighted Residual error estimates) framework is extended and further developed for single-phase flow problems in porous media. Specifically, our problem statement is the Biot system which consists of vector-valued displacements (geomechanics) coupled to a Darcy flow pressure equation. The MORe DWR method introduces a goal-oriented adaptive incremental proper orthogonal decomposition (POD) based-reduced-order model (ROM). The error in the reduced goal functional is estimated during the simulation, and the POD basis is enriched on-the-fly if the estimate exceeds a given threshold. This results in a reduction of the total number of full-order-model solves for the simulation of the porous medium, a robust estimation of the quantity of interest and well-suited reduced bases for the problem at hand. We apply a space-time Galerkin discretization with Taylor-Hood elements in space and a discontinuous Galerkin method with piecewise constant functions in time. The latter is well-known to be similar to the backward Euler scheme. We demonstrate the efficiency of our method on the well-known two-dimensional Mandel benchmark and a three-dimensional footing problem.

## Introduction

Porous media problems have long-standing applications in subsurface modeling, ground-water flow, hydraulic fracturing, geothermal energy recovery, and nuclear waste storage. The resulting mathematical models yield coupled systems of partial differential equations (PDEs) from which one of the most well-known is the so-called Biot problem [18–21,71]. Developments over the last two decades have led to more complicated models, due to additional physics that can be incorporated. Nonstationary, nonlinear coupled PDE systems are often obtained. They are also sometimes subject to inequality constraints as in multiphysics phase-field fracture in porous media [75,78]. Despite advances in numerical solvers of iterative or multigrid type, high-performance parallel

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

<footnote>*Correspondence:
thomas.wick@ifam.uni-
hannover.de</footnote>

<footnote n="1">¹Leibniz Universität Hannover,
Institut für Angewandte
Mathematik, AG
Wissenschaftliches Rechnen,
Welfengarten 1, 30167
Hannover, Germany</footnote>

<footnote n="2">²Université Paris-Saclay,
CentraleSupélec, ENS
Paris-Saclay, CNRS, LMPS,
Laboratoire de Mécanique
Paris-Saclay, 91190
Gif-sur-Yvette, France</footnote>

<footnote n="3">³The University of Texas at
Austin, Oden Institute, Austin, TX
78712, USA</footnote>
lel computing, and improved hardware performance, the computational cost for solving such multiphysics problems remains high. Specific examples in the field of poroelasticity are [5, 9, 10, 22, 42, 44, 53, 54, 76, 77].

Model order reduction (MOR) and reduced-order modeling (ROM) techniques [14, 15, 30, 37, 47, 56, 57, 65, 68, 80] provide one possible solution to significantly reduce the computational cost. Therein, the problem is split into two phases, an offline phase (solving a costly original high-fidelity model) to build a representative reduced basis and an online phase in which the reduced-order model is solved fast. This splitting comes with an additional cost, which pays off when the original model needs to be solved several (hundreds and more) times as for example in uncertainty quantification, inverse modeling, and optimal control. However, providing a robust and problem-specific reduced basis for nonlinear and general problems is a challenge.

Other model-order reduction approaches are based on solving the Galerkin problem with a separation of variables [2]. For nonlinear problems, some robust methods have been proposed, but they require a highly intrusive and specific numerical framework [17], or less intrusive and more flexible approaches, but without a mathematical proof of convergence for nonlinear cases [16]. Thus, there is a need for reduced-order numerical strategies which are flexible, general and for which the error accompanying the reduced solution space can be controlled for any quantity of interest. Additionally, a goal-oriented version of the proper generalized decomposition (PGD) has been applied to stationary problems in [43]. In contrast to our work, an a priori computed dual solution was used to control the error of the primal reduced-order solution.

Thanks to our own recent advances in space-time modeling and goal-oriented space-time error control for nonstationary, coupled problems [64], an incremental error-controlled proper orthogonal decomposition (POD) based ROM method, i.e. MORe DWR (Model Order Reduction with Dual-Weighted Residual error estimates), was suggested for the heat and elastodynamics equations in [32]. Therein, satisfying results for adaptively refining the ROM basis according to some distributed-in-time goal functionals were obtained.

In the current work, the Biot system of poroelasticity, namely vector-valued displacements and a Darcy-type pressure equation, is considered and modeled in a space-time fashion. Specifically, the temporal discretization is based on discontinuous Galerkin (dG) finite elements, while the spatial discretization employs classical continuous Galerkin (cG) finite elements. This serves as a basis for space-time error-controlled adaptivity by applying the dual-weighted residual method (DWR) [8, 12] to our ROM concepts. Importantly, the ROM error is computed with respect to a goal functional, which is often motivated from technical quantities of interest arising in physics, engineering and practical applications. The ROM updates are performed in an incremental fashion as recently suggested in [46] and do not require a full FOM (full-order model) run in advance. Thus, the ROM is updated on the fly using FOM solutions estimated only for one temporal element. Here, the singular value decomposition (SVD) for computing the leading components reduces to a truncated version and provides the snapshot matrix. An important point in coupled problems resorts to considering a monolithic version and computing one single SVD or, alternatively, splitting the SVD into its PDE components. In this work, we aim for the latter and split the displacements and pressure into separate SVDs in order to facilitate computing different POD sizes for the individual solution components. These develop-
lel computing, and improved hardware performance, the computational cost for solving such multiphysics problems remains high. Specific examples in the field of poroelasticity are [5, 9, 10, 22, 42, 44, 53, 54, 76, 77].

Model order reduction (MOR) and reduced-order modeling (ROM) techniques [14, 15, 30, 37, 47, 56, 57, 65, 68, 80] provide one possible solution to significantly reduce the computational cost. Therein, the problem is split into two phases, an offline phase (solving a costly original high-fidelity model) to build a representative reduced basis and an online phase in which the reduced-order model is solved fast. This splitting comes with an additional cost, which pays off when the original model needs to be solved several (hundreds and more) times as for example in uncertainty quantification, inverse modeling, and optimal control. However, providing a robust and problem-specific reduced basis for nonlinear and general problems is a challenge.

Other model-order reduction approaches are based on solving the Galerkin problem with a separation of variables [2]. For nonlinear problems, some robust methods have been proposed, but they require a highly intrusive and specific numerical framework [17], or less intrusive and more flexible approaches, but without a mathematical proof of convergence for nonlinear cases [16]. Thus, there is a need for reduced-order numerical strategies which are flexible, general and for which the error accompanying the reduced solution space can be controlled for any quantity of interest. Additionally, a goal-oriented version of the proper generalized decomposition (PGD) has been applied to stationary problems in [43]. In contrast to our work, an a priori computed dual solution was used to control the error of the primal reduced-order solution.

Thanks to our own recent advances in space-time modeling and goal-oriented space-time error control for nonstationary, coupled problems [64], an incremental error-controlled proper orthogonal decomposition (POD) based ROM method, i.e. MORe DWR (Model Order Reduction with Dual-Weighted Residual error estimates), was suggested for the heat and elastodynamics equations in [32]. Therein, satisfying results for adaptively refining the ROM basis according to some distributed-in-time goal functionals were obtained.

In the current work, the Biot system of poroelasticity, namely vector-valued displacements and a Darcy-type pressure equation, is considered and modeled in a space-time fashion. Specifically, the temporal discretization is based on discontinuous Galerkin (dG) finite elements, while the spatial discretization employs classical continuous Galerkin (cG) finite elements. This serves as a basis for space-time error-controlled adaptivity by applying the dual-weighted residual method (DWR) [8, 12] to our ROM concepts. Importantly, the ROM error is computed with respect to a goal functional, which is often motivated from technical quantities of interest arising in physics, engineering and practical applications. The ROM updates are performed in an incremental fashion as recently suggested in [46] and do not require a full FOM (full-order model) run in advance. Thus, the ROM is updated on the fly using FOM solutions estimated only for one temporal element. Here, the singular value decomposition (SVD) for computing the leading components reduces to a truncated version and provides the snapshot matrix. An important point in coupled problems resorts to considering a monolithic version and computing one single SVD or, alternatively, splitting the SVD into its PDE components. In this work, we aim for the latter and split the displacements and pressure into separate SVDs in order to facilitate computing different POD sizes for the individual solution components. These develop-
ments result into a final MORe DWR algorithm, which is tested for two-dimensional and three-dimensional configurations. First, the well-known Mandel benchmark is addressed. In a second numerical test, a poroelastic three-dimensional block is considered. Therein, on the one hand, typical cost factors such as speedup, FOM solves and ROM sizes are investigated. On the other hand, standard quantities of a posteriori error control such as effectivity indices and indicator indices are computed as well, demonstrating the performance of our overall proposed concepts.

The outline of this paper is as follows. In Sect. 2 the strong form problem statement is presented first. Then, a space-time formulation is derived. The dual-weighted residual method for goal-oriented error control is exposed, which requires a linear, backward-in-time running, adjoint problem. Next, in Sect. 3, the MORe DWR (Model Order Reduction with Dual-Weighted Residual error estimates) concepts are explained. Goal-oriented error control for refining the ROM basis is also discussed. In Sect. 4 two numerical examples are investigated. The first test is the two-dimensional Mandel benchmark. The second test is a three-dimensional footing problem.

## Problem formulation

In this work, we model a porous medium by coupling a vector-valued displacement equation (geomechanics) to Darcy flow in a poro-elastic medium. In the problem description, $I := (0, T)$ denotes the temporal domain and $\Omega \subset \mathbb{R}^d$ a sufficiently smooth spatial domain with spatial dimension $d \in \{2, 3\}$.

## Strong form of poroelasticity and function spaces

The governing equations for poroelasticity read [28,49]: Find pressure $p: \tilde{\Omega} \times \tilde{I} \to \mathbb{R}$ and displacement $u: \tilde{\Omega} \times \tilde{I} \to \mathbb{R}^d$ such that

$$ \begin{aligned} \partial_t(c p + \alpha(\nabla_x \cdot u)) - \frac{1}{\nu}\nabla_x \cdot (K \nabla_x p) &= 0 & \text{in } \Omega \times I, \\ -\nabla_x \cdot \sigma(u) + \alpha \nabla_x p &= 0 & \text{in } \Omega \times I, \end{aligned} \quad (1) $$

with the isotropic stress tensor

$$ \sigma(u) := \mu(\nabla_x u + (\nabla_x u)^T) + \lambda(\nabla_x u)I, $$

and the (constrained specific) storage coefficient $c \ge c^* > 0$, which is assumed strictly positive in this paper, the Biot-Willis constant $\alpha \in [0, 1]$, the permeability tensor $K$, the fluid viscosity $\nu$ and the Lamé parameters $\lambda, \mu > 0$. We notice that the constant $c$ may depend on space, i.e., $c(x)$, and is linked to the compressibility $M > 0$. When $c$ tends to zero, numerical instabilities may arise (e.g., [58,60]) for which mixed methods or enriched Galerkin or discontinuous Galerkin finite elements in the numerical discretization should be employed [48,50,55,59,69]. This coupled system of equations is also known as the Biot system [18–21,71]. A rigorous mathematical analysis of this poroelasticity problem can be found in [67].

Assuming homogeneous Dirichlet boundary conditions for the displacement on $\Gamma_D$ and inhomogeneous Neumann/traction boundary conditions
ments result into a final MORe DWR algorithm, which is tested for two-dimensional and three-dimensional configurations. First, the well-known Mandel benchmark is addressed. In a second numerical test, a poroelastic three-dimensional block is considered. Therein, on the one hand, typical cost factors such as speedup, FOM solves and ROM sizes are investigated. On the other hand, standard quantities of a posteriori error control such as effectivity indices and indicator indices are computed as well, demonstrating the performance of our overall proposed concepts.

The outline of this paper is as follows. In Sect. 2 the strong form problem statement is presented first. Then, a space-time formulation is derived. The dual-weighted residual method for goal-oriented error control is exposed, which requires a linear, backward-in-time running, adjoint problem. Next, in Sect. 3, the MORe DWR (Model Order Reduction with Dual-Weighted Residual error estimates) concepts are explained. Goal-oriented error control for refining the ROM basis is also discussed. In Sect. 4 two numerical examples are investigated. The first test is the two-dimensional Mandel benchmark. The second test is a three-dimensional footing problem.

## Problem formulation

In this work, we model a porous medium by coupling a vector-valued displacement equation (geomechanics) to Darcy flow in a poro-elastic medium. In the problem description, $I := (0, T)$ denotes the temporal domain and $\Omega \subset \mathbb{R}^d$ a sufficiently smooth spatial domain with spatial dimension $d \in \{2, 3\}$.

## Strong form of poroelasticity and function spaces

The governing equations for poroelasticity read [28,49]: Find pressure $p: \tilde{\Omega} \times \tilde{I} \to \mathbb{R}$ and displacement $u: \tilde{\Omega} \times \tilde{I} \to \mathbb{R}^d$ such that

$$ \begin{aligned} \partial_t(c p + \alpha(\nabla_x \cdot u)) - \frac{1}{\nu}\nabla_x \cdot (K \nabla_x p) &= 0 && \text{in } \Omega \times I, \\ -\nabla_x \cdot \sigma(u) + \alpha \nabla_x p &= 0 && \text{in } \Omega \times I, \end{aligned} \tag{1} $$

with the isotropic stress tensor

$$ \sigma(u) := \mu(\nabla_x u + (\nabla_x u)^T) + \lambda(\nabla_x u)I, $$

and the (constrained specific) storage coefficient $c \ge c^* > 0$, which is assumed strictly positive in this paper, the Biot-Willis constant $\alpha \in [0, 1]$, the permeability tensor $K$, the fluid viscosity $\nu$ and the Lamé parameters $\lambda, \mu > 0$. We notice that the constant $c$ may depend on space, i.e., $c(x)$, and is linked to the compressibility $M > 0$. When $c$ tends to zero, numerical instabilities may arise (e.g., [58,60]) for which mixed methods or enriched Galerkin or discontinuous Galerkin finite elements in the numerical discretization should be employed [48,50,55,59,69]. This coupled system of equations is also known as the Biot system [18–21,71]. A rigorous mathematical analysis of this poroelasticity problem can be found in [67].

Assuming homogeneous Dirichlet boundary conditions for the displacement on $\Gamma_D$ and inhomogeneous Neumann/traction boundary conditions
$$ \sigma(u) \cdot n = t $$

on $\Gamma_N = \partial\Omega \setminus \Gamma_D$, the spatial function spaces are given by

$$ V_u(\Omega) = (H_{0,\Gamma_D}^1(\Omega))^d, $$

$$ V_p(\Omega) = H^1(\Omega), $$

$$ V(\Omega) = V_u(\Omega) \times V_p(\Omega). $$

All boundary conditions, including pressure conditions as well, are carefully specified in Sect. 4 and the function spaces are redefined w.r.t. the modified boundary conditions, too.

### Weak space-time form of the primal problem

We define the space-time function space $X(I, V(\Omega))$ as

$$ X(I, V(\Omega)) := L^2(I, V(\Omega)) \cap H^1(I, V^*(\Omega)) $$

with the dual space of $V(\Omega)$ being denoted as $V^*(\Omega) = L(V(\Omega), \mathbb{R})$. Furthermore, we will use the notation

$$ (f,g) := (f,g)_{L^2(\Omega)} := \int_{\Omega} f \cdot g \, dx, \quad \langle (f,g) \rangle := (f,g)_{L^2(I,L^2(\Omega))} := \int_{I} (f,g) \, dt, $$

$$ \langle f,g \rangle := \langle f,g \rangle_{L^2(\Gamma)} := \int_{\Gamma} f \cdot g \, ds, \quad \langle (f,g) \rangle := (f,g)_{L^2(I,L^2(\Gamma))} := \int_{I} \langle f,g \rangle \, dt. $$

In this notation, $f \cdot g$ represents the Euclidean inner product if $f$ and $g$ are scalar- or vector-valued and represents the Frobenius inner product if $f$ and $g$ are matrices.

We can now derive the space-time variational formulation for the (primal) poroelasticity problem. By integration by parts, the variational formulation reads: Find $U := \{u, p\} \in X(I, V(\Omega))$ such that

$$ c(\langle \partial_t p, \phi^p \rangle) + \alpha(\langle \partial_t (\nabla_x \cdot u), \phi^p \rangle) + \frac{K}{\nu}(\langle \nabla_x p, \nabla_x \phi^p \rangle) = 0, $$

$$ (\langle \sigma(u), \nabla_x \phi^u \rangle - \alpha \langle pI, \nabla_x \phi^u \rangle) + \alpha \langle pn, \phi^u \rangle_{\Gamma_N \times I} = \langle t, \phi^u \rangle_{\Gamma_N \times I} $$

$$ \forall \Phi := \begin{pmatrix} \phi^u \\ \phi^p \end{pmatrix} \in X(I, V(\Omega)). $$

Here, $n$ denotes the normal vector on the Neumann boundary $\Gamma_N$.

For the time discretization, let $\mathcal{T}_k := \{I_m \mid 1 \le m \le M\}$ with $I_m := (t_{m-1}, t_m)$ and

$$ 0 := t_0 < t_1 < \dots < t_{M-1} < t_M =: T $$

be a partitioning of time i.e., $\tilde{I} = [0, T] = \cup_{m=1}^M \tilde{I}_m$. Then, in an intermediate step similar to spatial discontinuous Galerkin (dG) method, we can define a broken function space that allows for discontinuities at the temporal grid points:

$$ \tilde{X}(\mathcal{T}_k, V(\Omega)) := \left\{ \begin{pmatrix} u \\ p \end{pmatrix} \in L^2(I, L^2(\Omega)^{d+1}), \begin{pmatrix} u \\ p \end{pmatrix} \Big|_{I_m} \in X(I_m, V(\Omega)) \forall I_m \in \mathcal{T}_k \right\}. $$
$$ \sigma(u) \cdot n = t $$

on $\Gamma_N = \partial\Omega \setminus \Gamma_D$, the spatial function spaces are given by

$$ V_u(\Omega) = (H_{0,\Gamma_D}^1(\Omega))^d, $$

$$ V_p(\Omega) = H^1(\Omega), $$

$$ V(\Omega) = V_u(\Omega) \times V_p(\Omega). $$

All boundary conditions, including pressure conditions as well, are carefully specified in Sect. 4 and the function spaces are redefined w.r.t. the modified boundary conditions, too.

### Weak space-time form of the primal problem

We define the space-time function space $X(I, V(\Omega))$ as

$$ X(I, V(\Omega)) := L^2(I, V(\Omega)) \cap H^1(I, V^*(\Omega)) $$

with the dual space of $V(\Omega)$ being denoted as $V^*(\Omega) = L(V(\Omega), \mathbb{R})$. Furthermore, we will use the notation

$$ (f,g) := (f,g)_{L^2(\Omega)} := \int_{\Omega} f \cdot g \, dx, \quad \langle (f,g) \rangle := (f,g)_{L^2(I,L^2(\Omega))} := \int_{I} (f,g) \, dt, $$

$$ \langle f,g \rangle := \langle f,g \rangle_{L^2(\Gamma)} := \int_{\Gamma} f \cdot g \, ds, \quad \langle (f,g) \rangle := (f,g)_{L^2(I,L^2(\Gamma))} := \int_{I} \langle f,g \rangle \, dt. $$

In this notation, $f \cdot g$ represents the Euclidean inner product if $f$ and $g$ are scalar- or vector-valued and represents the Frobenius inner product if $f$ and $g$ are matrices.

We can now derive the space-time variational formulation for the (primal) poroelasticity problem. By integration by parts, the variational formulation reads: Find $U := \{u, p\} \in X(I, V(\Omega))$ such that

$$ c(\langle \partial_t p, \phi^p \rangle) + \alpha(\langle \partial_t (\nabla_x \cdot u), \phi^p \rangle) + \frac{K}{\nu}(\langle \nabla_x p, \nabla_x \phi^p \rangle) = 0, $$

$$ (\langle \sigma(u), \nabla_x \phi^u \rangle - \alpha(\langle pI, \nabla_x \phi^u \rangle) + \alpha(\langle pn, \phi^u \rangle))_{\Gamma_N \times I} = (\langle t, \phi^u \rangle)_{\Gamma_N \times I} $$

$$ \forall \Phi := \begin{pmatrix} \phi^u \\ \phi^p \end{pmatrix} \in X(I, V(\Omega)). $$

Here, $n$ denotes the normal vector on the Neumann boundary $\Gamma_N$.

For the time discretization, let $\mathcal{T}_k := \{I_m \mid 1 \le m \le M\}$ with $I_m := (t_{m-1}, t_m)$ and

$$ 0 := t_0 < t_1 < \dots < t_{M-1} < t_M := T $$

be a partitioning of time i.e., $\tilde{I} = [0, T] = \cup_{m=1}^M \tilde{I}_m$. Then, in an intermediate step similar to spatial discontinuous Galerkin (dG) method, we can define a broken function space that allows for discontinuities at the temporal grid points:

$$ \tilde{X}(\mathcal{T}_k, V(\Omega)) := \left\{ \begin{pmatrix} u \\ p \end{pmatrix} \in L^2(I, L^2(\Omega)^{d+1}), \begin{pmatrix} u \\ p \end{pmatrix} \Big|_{I_m} \in X(I_m, V(\Omega)) \forall I_m \in \mathcal{T}_k \right\}. $$
Due to these discontinuities, the limits of a function $f$ at time $t_m$ from above and from below are

$$f_m^{\pm} := \lim_{\epsilon \searrow 0} f(t_m \pm \epsilon),$$

and the jump of the function value of $f$ at time $t_m$ is

$$[f]_m := f_m^+ - f_m^-.$$

Accounting for the discontinuities in time, we thus need to solve the problem:

Find $U \in \tilde{X}(\mathcal{T}_k, V(\Omega))$ such that

$$A(U)(\Phi) = F(\Phi) \quad \forall \Phi \in \tilde{X}(\mathcal{T}_k, V(\Omega)). \qquad (2)$$

The bilinear form and right-hand side read

$$A(U)(\Phi) = \sum_{m=1}^{M} \int_{I_m} c(\partial_t p, \phi^p) + \alpha(\partial_t(\nabla_x \cdot u), \phi^p) + \frac{K}{\nu}(\nabla_x p, \nabla_x \phi^p) \, dt \\ + \sum_{m=1}^{M} \int_{I_m} (\sigma(u), \nabla_x \phi^u) - \alpha(pI, \nabla_x \phi^u) + \alpha(pn, \phi^u)_{\Gamma_N} \, dt \\ + \sum_{m=1}^{M-1} \alpha([\nabla_x \cdot u]_{m\phi}, \phi_m^{p,+}) + \alpha(\nabla_x \cdot u_0^+, \phi_0^{p,+}) \\ + \sum_{m=1}^{M-1} c([p]_{m\phi}, \phi_m^{p,+}) + c(p_0^+, \phi_0^{p,+})$$

and

$$F(\Phi) := (\langle t, \phi^u \rangle)_{\Gamma_N \times I} + \alpha(\nabla_x \cdot u^0, \phi_0^{p,+}) + c(p^0, \phi_0^{p,+}),$$

where

$$U := (u, p), \quad \Phi := (\phi^u, \phi^p).$$

**Remark 2.1** This space-time formulation can also be easily extended to a multirate-in-time setting with different timestep sizes for displacement and pressure [63]; a rigorous analysis starting directly with the backward Euler time discretization and mixed spaces for flow and conformal Galerkin for geomechanics is provided in [3].

For the projection of the problem onto the POD vectors, the block structure of the linear system is exploited to build separate bases for displacement and pressure. We can decompose the bilinear form into the four following blocks

$$A_1(u)(\phi^u) = \sum_{m=1}^{M} \int_{I_m} (\sigma(u), \nabla_x \phi^u) \, dt,$$

$$A_2(p)(\phi^p) = \sum_{m=1}^{M} \int_{I_m} c(\partial_t p, \phi^p) + \frac{K}{\nu}(\nabla_x p, \nabla_x \phi^p) \, dt$$
Due to these discontinuities, the limits of a function $f$ at time $t_m$ from above and from below are

$$f_m^{\pm} := \lim_{\epsilon \searrow 0} f(t_m \pm \epsilon),$$

and the jump of the function value of $f$ at time $t_m$ is

$$[f]_m := f_m^+ - f_m^-.$$

Accounting for the discontinuities in time, we thus need to solve the problem:

Find $U \in \tilde{X}(\mathcal{T}_k, V(\Omega))$ such that

$$A(U)(\Phi) = F(\Phi) \quad \forall \Phi \in \tilde{X}(\mathcal{T}_k, V(\Omega)). \qquad (2)$$

The bilinear form and right-hand side read

$$A(U)(\Phi) = \sum_{m=1}^{M} \int_{I_m} c(\partial_t p, \phi^p) + \alpha(\partial_t(\nabla_x \cdot u), \phi^p) + \frac{K}{\nu}(\nabla_x p, \nabla_x \phi^p) \, dt \\ + \sum_{m=1}^{M} \int_{I_m} (\sigma(u), \nabla_x \phi^u) - \alpha(pI, \nabla_x \phi^u) + \alpha(pn, \phi^u)_{\Gamma_N} \, dt \\ + \sum_{m=1}^{M-1} \alpha([\nabla_x \cdot u]_{m\phi}, \phi_m^{p,+}) + \alpha(\nabla_x \cdot u_0^+, \phi_0^{p,+}) \\ + \sum_{m=1}^{M-1} c([p]_{m\phi}, \phi_m^{p,+}) + c(p_0^+, \phi_0^{p,+})$$

and

$$F(\Phi) := (\langle t, \phi^u \rangle)_{\Gamma_N \times I} + \alpha(\nabla_x \cdot u^0, \phi_0^{p,+}) + c(p^0, \phi_0^{p,+}),$$

where

$$U := (u, p), \quad \Phi := (\phi^u, \phi^p).$$

**Remark 2.1** This space-time formulation can also be easily extended to a multirate-in-time setting with different timestep sizes for displacement and pressure [63]; a rigorous analysis starting directly with the backward Euler time discretization and mixed spaces for flow and conformal Galerkin for geomechanics is provided in [3].

For the projection of the problem onto the POD vectors, the block structure of the linear system is exploited to build separate bases for displacement and pressure. We can decompose the bilinear form into the four following blocks

$$A_1(u)(\phi^u) = \sum_{m=1}^{M} \int_{I_m} (\sigma(u), \nabla_x \phi^u) \, dt,$$

$$A_2(p)(\phi^p) = \sum_{m=1}^{M} \int_{I_m} c(\partial_t p, \phi^p) + \frac{K}{\nu}(\nabla_x p, \nabla_x \phi^p) \, dt$$
$$ 
\begin{aligned}
&+ \sum_{m=1}^{M-1} c([p]_{m\nu}, \phi_m^{p,+}) + c(p_0^+, \phi_0^{p,+}), \\
B_1(p)(\phi^u) &= \sum_{m=1}^{M} \int_{I_m} -\alpha(pI, \nabla_x\phi^u) + \alpha(pn, \phi^u)_{\Gamma_N} \, dt, \\
B_2(u)(\phi^p) &= \sum_{m=1}^{M} \int_{I_m} \alpha(\partial_t(\nabla_x \cdot u), \phi^p) \, dt \\
&\quad + \sum_{m=1}^{M-1} \alpha([\nabla_x \cdot u]_{m\nu}, \phi_m^{p,+}) + \alpha(\nabla_x \cdot u_0^+, \phi_0^{p,+}).
\end{aligned}
$$

### dG(0), i.e., backward Euler, time discretization of the primal problem

Using the space-time formulation of the poroelasticity problem (2), a dG(0) time discretization is derived by using discontinuous, piecewise-constant finite elements in time; see e.g., [79]. This results into a backward Euler scheme. Then, on each temporal element $I_m = (t_{m-1}, t_m)$, we have temporally constant functions $(u_m, p_m) \in V(\Omega)$ such that

$$ u|_{I_m} := u_m, \quad p|_{I_m} := p_m. $$

Using piecewise-constant functions $f$ in time, we can insert the relations

$$ \partial_t f|_{I_m} = 0, \quad [f]_m = f_{m+1} - f_m $$

and arrive at the time-stepping scheme: Find $U_m := (u_m, p_m) \in V(\Omega)$ such that

$$ 
\begin{aligned}
c(p_m - p_{m-1}, \phi^p) + \alpha(\nabla_x \cdot u_m - \nabla_x \cdot u_{m-1}, \phi^p) + k \frac{K}{\nu} (\nabla_x p_m, \nabla_x \phi^p) \\
+ k(\sigma(u_m), \nabla_x \phi^u) - k\alpha(p_m I, \nabla_x \phi^u) + k\alpha(p_m n, \phi^u)_{\Gamma_N} = k(t, \phi^u)_{\Gamma_N} \\
\forall \Phi = (\phi^u, \phi^p) \in V(\Omega).
\end{aligned}
$$

Here, $k := t_m - t_{m-1}$ denotes a (constant) timestep size.

**Remark 2.2** Instead of starting from the space-time formulation, we can also directly derive the backward Euler time discretization from the strong formulation of poroelasticity (1). The main difference is that now we do not scale the mechanics equation by the timestep size $k$. This is motivated by the displacement equation being quasi-static. Hence, in traditional time discretizations the time integral is omitted for this equation. For a consistent mathematical description, the previous time-stepping scheme fits better in the space-time setting, but in the actual implementation we can neglect the time step size $k$ in the mechanics equation. We then have: Find $U_m := (u_m, p_m) \in V(\Omega)$ such that

$$ 
\begin{aligned}
c(p_m - p_{m-1}, \phi^p) + \alpha(\nabla_x \cdot u_m - \nabla_x \cdot u_{m-1}, \phi^p) + k \frac{K}{\nu} (\nabla_x p_m, \nabla_x \phi^p) \\
+ (\sigma(u_m), \nabla_x \phi^u) - \alpha(p_m I, \nabla_x \phi^u) + \alpha(p_m n, \phi^u)_{\Gamma_N} = \langle t, \phi^u \rangle_{\Gamma_N} \\
\forall \Phi = (\phi^u, \phi^p) \in V(\Omega).
\end{aligned}
$$

We will use this formulation for our numerical tests since it is equivalent to the dG(0) time discretization.
$$ 
\begin{aligned}
& + \sum_{m=1}^{M-1} c([p]_m, \phi_m^{p,+}) + c(p_0^+, \phi_0^{p,+}), \\
B_1(p)(\phi^u) &= \sum_{m=1}^{M} \int_{I_m} -\alpha(pI, \nabla_x \phi^u) + \alpha(pn, \phi^u)_{\Gamma_N} \, dt, \\
B_2(u)(\phi^p) &= \sum_{m=1}^{M} \int_{I_m} \alpha(\partial_t(\nabla_x \cdot u), \phi^p) \, dt \\
&\quad + \sum_{m=1}^{M-1} \alpha([\nabla_x \cdot u]_m, \phi_m^{p,+}) + \alpha(\nabla_x \cdot u_0^+, \phi_0^{p,+}).
\end{aligned}
$$

### dG(0), i.e., backward Euler, time discretization of the primal problem

Using the space-time formulation of the poroelasticity problem (2), a dG(0) time discretization is derived by using discontinuous, piecewise-constant finite elements in time; see e.g., [79]. This results into a backward Euler scheme. Then, on each temporal element $I_m = (t_{m-1}, t_m)$, we have temporally constant functions $(u_m, p_m) \in V(\Omega)$ such that

$$ u|_{I_m} := u_m, \quad p|_{I_m} := p_m. $$

Using piecewise-constant functions $f$ in time, we can insert the relations

$$ \partial_t f|_{I_m} = 0, \quad [f]_m = f_{m+1} - f_m $$

and arrive at the time-stepping scheme: Find $U_m := (u_m, p_m) \in V(\Omega)$ such that

$$ 
\begin{aligned}
c(p_m - p_{m-1}, \phi^p) + \alpha(\nabla_x \cdot u_m - \nabla_x \cdot u_{m-1}, \phi^p) + k \frac{K}{\nu} (\nabla_x p_m, \nabla_x \phi^p) \\
+ k(\sigma(u_m), \nabla_x \phi^u) - k\alpha(p_m I, \nabla_x \phi^u) + k\alpha(p_m n, \phi^u)_{\Gamma_N} = k(t, \phi^u)_{\Gamma_N} \\
\forall \Phi = (\phi^u, \phi^p) \in V(\Omega).
\end{aligned}
$$

Here, $k := t_m - t_{m-1}$ denotes a (constant) timestep size.

**Remark 2.2** Instead of starting from the space-time formulation, we can also directly derive the backward Euler time discretization from the strong formulation of poroelasticity (1). The main difference is that now we do not scale the mechanics equation by the timestep size $k$. This is motivated by the displacement equation being quasi-static. Hence, in traditional time discretizations the time integral is omitted for this equation. For a consistent mathematical description, the previous time-stepping scheme fits better in the space-time setting, but in the actual implementation we can neglect the time step size $k$ in the mechanics equation. We then have: Find $U_m := (u_m, p_m) \in V(\Omega)$ such that

$$ 
\begin{aligned}
c(p_m - p_{m-1}, \phi^p) + \alpha(\nabla_x \cdot u_m - \nabla_x \cdot u_{m-1}, \phi^p) + k \frac{K}{\nu} (\nabla_x p_m, \nabla_x \phi^p) \\
+ (\sigma(u_m), \nabla_x \phi^u) - \alpha(p_m I, \nabla_x \phi^u) + \alpha(p_m n, \phi^u)_{\Gamma_N} = (t, \phi^u)_{\Gamma_N} \\
\forall \Phi = (\phi^u, \phi^p) \in V(\Omega).
\end{aligned}
$$

We will use this formulation for our numerical tests since it is equivalent to the dG(0) time discretization.
## Weak space-time form of the adjoint problem

The MORE DWR method [32] measures the adjoint sensitivity of the primal reduced-order solution with respect to some quantities of interest. Let a goal functional $J: \tilde{X}(\mathcal{T}_k, V(\Omega)) \to \mathbb{R}$ of the form

$$J(U) = \int_{0}^{T} J_{1}(U(t)) \, dt + J_{2}(U(T)), \qquad (3)$$

be given, which represents some physical quantity of interest (QoI). The adjoint equations for porous media are derived, following [79]. Since the problem is linear, we just need to switch trial and test functions in the bilinear form and replace the primal with the adjoint solution. The forcing function on the right-hand side is being substituted by the goal functional $J$, which is justified by the Lagrangian formalism, c.f. Sect. 3. The space-time formulation of the adjoint problem thus reads: Find $Z \in \tilde{X}(\mathcal{T}_k, V(\Omega))$ such that

$$A(\Phi)(Z) = J'_{U}(U)(\Phi) \quad \forall \Phi \in \tilde{X}(\mathcal{T}_k, V(\Omega)).$$

More concretely, the bilinear form reads

$$A(\Phi)(Z) = \sum_{m=1}^{M} \int_{I_m} c(\partial_t \phi^p, z^p) + \alpha(\partial_t (\nabla_x \cdot \phi^u), z^p) + \frac{K}{\nu} (\nabla_x \phi^p, \nabla_x z^p) \, dt \\ + \sum_{m=1}^{M} \int_{I_m} (\sigma(\phi^u), \nabla_x z^u) - \alpha(\phi^p I, \nabla_x z^u) + \alpha(\phi^p n, z^u)_{\Gamma_N} \, dt \\ + \sum_{m=1}^{M-1} \alpha([\nabla_x \cdot \phi^u]_{mv}, z_m^{p,+}) + \alpha(\nabla_x \cdot \phi_0^{u,+}, z_0^{p,+}) \\ + \sum_{m=1}^{M-1} c([\phi^p]_{mv}, z_m^{p,+}) + c(\phi_0^{p,+}, z_0^{p,+}),$$

where

$$Z := (z^u, z^p), \quad \Phi := (\phi^u, \phi^p).$$

The time derivatives are moved from the test function to the adjoint solution by integration by parts in time e.g., on the time-continuous level we use integration by parts for

$$\int_{0}^{T} \alpha(\partial_t (\nabla_x \cdot \phi^u), z^p) = - \int_{0}^{T} \alpha(\nabla_x \cdot \phi^u, \partial_t z^p) \, dt + \alpha(\nabla_x \cdot \phi^u(T), z^p(T)) \\ - \alpha(\nabla_x \cdot \phi^u(0), z^p(0)).$$

The adjoint problem runs backward in time and starts with a final time condition, which depends on the quantity of interest. In particular, we have

$$\alpha(\nabla_x \cdot \phi^u(T), z^p(T)) + c(\phi^p(T), z^p(T)) = J'_{2,U}(U(T))(\Phi(T))$$
## Weak space-time form of the adjoint problem

The MORE DWR method [32] measures the adjoint sensitivity of the primal reduced-order solution with respect to some quantities of interest. Let a goal functional $J: \tilde{X}(\mathcal{T}_k, V(\Omega)) \to \mathbb{R}$ of the form

$$J(U) = \int_{0}^{T} J_{1}(U(t)) \, dt + J_{2}(U(T)), \qquad (3)$$

be given, which represents some physical quantity of interest (QoI). The adjoint equations for porous media are derived, following [79]. Since the problem is linear, we just need to switch trial and test functions in the bilinear form and replace the primal with the adjoint solution. The forcing function on the right-hand side is being substituted by the goal functional $J$, which is justified by the Lagrangian formalism, c.f. Sect. 3. The space-time formulation of the adjoint problem thus reads: Find $Z \in \tilde{X}(\mathcal{T}_k, V(\Omega))$ such that

$$A(\Phi)(Z) = J'_{U}(U)(\Phi) \quad \forall \Phi \in \tilde{X}(\mathcal{T}_k, V(\Omega)).$$

More concretely, the bilinear form reads

$$A(\Phi)(Z) = \sum_{m=1}^{M} \int_{I_m} c(\partial_t \phi^p, z^p) + \alpha(\partial_t (\nabla_x \cdot \phi^u), z^p) + \frac{K}{\nu} (\nabla_x \phi^p, \nabla_x z^p) \, dt \\ + \sum_{m=1}^{M} \int_{I_m} (\sigma(\phi^u), \nabla_x z^u) - \alpha(\phi^p I, \nabla_x z^u) + \alpha \langle \phi^p n, z^u \rangle_{\Gamma_N} \, dt \\ + \sum_{m=1}^{M-1} \alpha([\nabla_x \cdot \phi^u]_{mv}, z_m^{p,+}) + \alpha(\nabla_x \cdot \phi_0^{u,+}, z_0^{p,+}) \\ + \sum_{m=1}^{M-1} c([\phi^p]_{mv}, z_m^{p,+}) + c(\phi_0^{p,+}, z_0^{p,+}),$$

where

$$Z := (z^u, z^p), \quad \Phi := (\phi^u, \phi^p).$$

The time derivatives are moved from the test function to the adjoint solution by integration by parts in time e.g., on the time-continuous level we use integration by parts for

$$\int_{0}^{T} \alpha(\partial_{t}(\nabla_{x} \cdot \phi^{u}), z^{p}) = - \int_{0}^{T} \alpha(\nabla_{x} \cdot \phi^{u}, \partial_{t} z^{p}) \, dt + \alpha(\nabla_{x} \cdot \phi^{u}(T), z^{p}(T)) \\ - \alpha(\nabla_{x} \cdot \phi^{u}(0), z^{p}(0)).$$

The adjoint problem runs backward in time and starts with a final time condition, which depends on the quantity of interest. In particular, we have

$$\alpha(\nabla_x \cdot \phi^u(T), z^p(T)) + c(\phi^p(T), z^p(T)) = J'_{2,U}(U(T))(\Phi(T))$$
i.e., Z(T) = 0 for $J_2 = 0$. The integration by parts can be applied to each temporal element, which leads exemplarily to

$$
\begin{align*}
& \sum_{m=1}^{M} \int_{I_m} \alpha(\partial_t(\nabla_x \cdot \phi^u), z^p) \, dt + \sum_{m=1}^{M-1} \alpha([\nabla_x \cdot \phi^u]_{m\nu} z_m^{p,+}) + \alpha(\nabla_x \cdot \phi_0^{u,+}, z_0^{p,+}) \\
&= -\sum_{m=1}^{M} \int_{I_m} \alpha(\nabla_x \cdot \phi^u, \partial_t z^p) \, dt + \sum_{m=1}^{M-1} \alpha([\nabla_x \cdot \phi^u]_{m\nu} z_m^{p,+}) + \alpha(\nabla_x \cdot \phi_0^{u,+}, z_0^{p,+}) \\
&\quad + \sum_{m=1}^{M} \alpha(\nabla_x \cdot \phi_m^{u,-}, z_m^{p,-}) - \alpha(\nabla_x \cdot \phi_{m-1}^{u,+}, z_{m-1}^{p,+}) \\
&= -\sum_{m=1}^{M} \int_{I_m} \alpha(\nabla_x \cdot \phi^u, \partial_t z^p) \, dt + \sum_{m=1}^{M-1} \alpha(\alpha(\overline{\nabla_x \cdot \phi_m^{u,+}}, z_m^{p,+}) - \alpha(\nabla_x \cdot \phi_m^{u,-}, z_m^{p,+}) \\
&\quad + \alpha(\nabla_x \cdot \phi_m^{u,-}, z_m^{p,-}) - \alpha(\overline{\nabla_x \cdot \phi_m^{u,+}}, z_m^{p,+})) \\
&\quad + \alpha(\nabla_x \cdot \phi_M^{u,-}, z_M^{p,-}) - \alpha(\overline{\nabla_x \cdot \phi_0^{u,+}}, z_0^{p,+}) + \alpha(\overline{\nabla_x \cdot \phi_0^{u,+}}, z_0^{p,+}) \\
&= -\sum_{m=1}^{M} \int_{I_m} \alpha(\nabla_x \cdot \phi^u, \partial_t z^p) \, dt - \sum_{m=1}^{M-1} \alpha(\nabla_x \cdot \phi_m^{u,-}, [z^p]_m) + \alpha(\nabla_x \cdot \phi_M^{u,-}, z_M^{p,-}).
\end{align*}
$$

Overall, we can rewrite the left-hand side of the adjoint problem as

$$
\begin{align}
A(\Phi)(Z) ={}& \sum_{m=1}^{M} \int_{I_m} -c(\phi^p, \partial_t z^p) - \alpha(\nabla_x \cdot \phi^u, \partial_t z^p) + \frac{K}{\nu} (\nabla_x \phi^p, \nabla_x z^p) \, dt \nonumber \\
& + \sum_{m=1}^{M} \int_{I_m} (\sigma(\phi^u), \nabla_x z^u) - \alpha(\phi^p I, \nabla_x z^u) + \alpha \langle \phi^p n, z^u \rangle_{\Gamma_N} \, dt \nonumber \\
& - \sum_{m=1}^{M-1} \alpha(\nabla_x \cdot \phi_m^{u,-}, [z^p]_m) + \alpha(\nabla_x \cdot \phi_M^{u,-}, z_M^{p,-}) \nonumber \\
& - \sum_{m=1}^{M-1} c(\phi_m^{p,-}, [z^p]_m) + c(\phi_M^{p,-}, z_M^{p,-}). \tag{4}
\end{align}
$$

dG(0), i.e., backward Euler, time discretization of the adjoint problem

A dG(0) time discretization is again used for the adjoint problem. Then on each temporal element $I_{m+1} = (t_m, t_{m+1})$, we have temporally constant functions $(z_m^u, z_m^p) \in V(\Omega)$ such that

$$
z^u|_{I_{m+1}} =: z_m^u \qquad z^p|_{I_{m+1}} =: z_m^p.
$$

Using piecewise-constant functions $f$ in time, the relations

$$
\partial_t f|_{I_m} = 0, \quad [-f]_m = -(f_{m+1} - f_m) = f_m - f_{m+1}
$$

hold and give the adjoint time-stepping scheme: Find $Z_m := (z_m^u, z_m^p) \in V(\Omega)$ such that

$$
\begin{align*}
& c(z_m^p - z_{m+1}^p, \phi^p) + \alpha(z_m^p - z_{m+1}^p, \nabla_x \cdot \phi^u) + k \frac{K}{\nu} (\nabla_x z_m^p, \nabla_x \phi^p) \\
& \quad + k(\nabla_x z_m^u, \sigma(\phi^u)) - k\alpha(\nabla_x z_m^u, \phi^p I) + k\alpha \langle z_m^u, \phi^p n \rangle_{\Gamma_N} = J'_{U}(U)(\Phi) \\
& \forall \Phi = (\phi^u, \phi^p) \in V(\Omega).
\end{align*}
$$
i.e., Z(T) = 0 for $J_2 = 0$. The integration by parts can be applied to each temporal element, which leads exemplarily to

$$
\begin{aligned}
& \sum_{m=1}^{M} \int_{I_m} \alpha(\partial_t(\nabla_x \cdot \phi^u), z^p) \, dt + \sum_{m=1}^{M-1} \alpha([\nabla_x \cdot \phi^u]_{m\nu} z_m^{p+}) + \alpha(\nabla_x \cdot \phi_0^{u+}, z_0^{p+}) \\
& = -\sum_{m=1}^{M} \int_{I_m} \alpha(\nabla_x \cdot \phi^u, \partial_t z^p) \, dt + \sum_{m=1}^{M-1} \alpha([\nabla_x \cdot \phi^u]_{m\nu} z_m^{p+}) + \alpha(\nabla_x \cdot \phi_0^{u+}, z_0^{p+}) \\
& \qquad + \sum_{m=1}^{M} \alpha(\nabla_x \cdot \phi_m^{u-}, z_m^{p-}) - \alpha(\nabla_x \cdot \phi_{m-1}^{u+}, z_{m-1}^{p+}) \\
& = -\sum_{m=1}^{M} \int_{I_m} \alpha(\nabla_x \cdot \phi^u, \partial_t z^p) \, dt + \sum_{m=1}^{M-1} \alpha(\alpha(\overline{\nabla_x \cdot \phi_m^{u+}, z_m^{p+}}) - \alpha(\nabla_x \cdot \phi_m^{u-}, z_m^{p+}) \\
& \qquad + \alpha(\nabla_x \cdot \phi_m^{u-}, z_m^{p-}) - \alpha(\overline{\nabla_x \cdot \phi_m^{u+}, z_m^{p+}}) \\
& \qquad + \alpha(\nabla_x \cdot \phi_M^{u-}, z_M^{p-}) - \alpha(\overline{\nabla_x \cdot \phi_0^{u+}, z_0^{p+}}) + \alpha(\overline{\nabla_x \cdot \phi_0^{u+}, z_0^{p+}}) \\
& = -\sum_{m=1}^{M} \int_{I_m} \alpha(\nabla_x \cdot \phi^u, \partial_t z^p) \, dt - \sum_{m=1}^{M-1} \alpha(\nabla_x \cdot \phi_m^{u-}, [z^p]_m) + \alpha(\nabla_x \cdot \phi_M^{u-}, z_M^{p-}).
\end{aligned}
$$

Overall, we can rewrite the left-hand side of the adjoint problem as

$$
\begin{aligned}
A(\Phi)(Z) = & \sum_{m=1}^{M} \int_{I_m} -c(\phi^p, \partial_t z^p) - \alpha(\nabla_x \cdot \phi^u, \partial_t z^p) + \frac{K}{\nu}(\nabla_x \phi^p, \nabla_x z^p) \, dt \\
& + \sum_{m=1}^{M} \int_{I_m} (\sigma(\phi^u), \nabla_x z^u) - \alpha(\phi^p I, \nabla_x z^u) + \alpha \langle \phi^p n, z^u \rangle_{\Gamma_N} \, dt \\
& - \sum_{m=1}^{M-1} \alpha(\nabla_x \cdot \phi_m^{u-}, [z^p]_m) + \alpha(\nabla_x \cdot \phi_M^{u-}, z_M^{p-}) \\
& - \sum_{m=1}^{M-1} c(\phi_m^{p-}, [z^p]_m) + c(\phi_M^{p-}, z_M^{p-}).
\end{aligned}
\tag{4}
$$

dG(0), i.e., backward Euler, time discretization of the adjoint problem

A dG(0) time discretization is again used for the adjoint problem. Then on each temporal element $I_{m+1} = (t_m, t_{m+1})$, we have temporally constant functions $(z_m^u, z_m^p) \in V(\Omega)$ such that

$$z^u|_{I_{m+1}} =: z_m^u \qquad z^p|_{I_{m+1}} =: z_m^p.$$

Using piecewise-constant functions $f$ in time, the relations

$$\partial_t f|_{I_m} = 0, \quad [-f]_m = -(f_{m+1} - f_m) = f_m - f_{m+1}$$

hold and give the adjoint time-stepping scheme: Find $Z_m := (z_m^u, z_m^p) \in V(\Omega)$ such that

$$
\begin{aligned}
& c(z_m^p - z_{m+1}^p, \phi^p) + \alpha(z_m^p - z_{m+1}^p, \nabla_x \cdot \phi^u) + k \frac{K}{\nu} (\nabla_x z_m^p, \nabla_x \phi^p) \\
& \qquad + k(\nabla_x z_m^u, \sigma(\phi^u)) - k\alpha(\nabla_x z_m^u, \phi^p I) + k\alpha \langle z_m^u, \phi^p n \rangle_{\Gamma_N} = J'_{U}(U)(\Phi) \\
& \forall \Phi = (\phi^u, \phi^p) \in V(\Omega).
\end{aligned}
$$
**Remark 2.3** To be consistent with the backward Euler time discretization in Remark 2.2, we again get rid of the timestep size *k* in displacement related terms and obtain

$$
\begin{align*}
c(z_m^p - z_{m+1}^p, \phi^p) + \alpha(z_m^p - z_{m+1}^p, \nabla_x \cdot \phi^u) + k \frac{K}{\nu} (\nabla_x z_m^p, \nabla_x \phi^p) \\
+ (\nabla_x z_m^u, \sigma(\phi^u)) - \alpha(\nabla_x z_m^u, \phi^p I) + \alpha \langle z_m^u, \phi^p n \rangle_{\Gamma_N} &= J'_{U}(U)(\Phi) \\
\forall \Phi = (\phi^u, \phi^p) \in V(\Omega).
\end{align*}
$$

## MORe DWR: model order reduction with dual-weighted residual error estimates

The MORe DWR approach aims at evaluating the quantities of interest using a reduced-order model, while guaranteeing the error due to this approximation. Thus, the difference between the reduced-order-model (ROM) solution $U^{\text{ROM}} := U_{kh}^{\text{ROM}} \in X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\text{ROM}})$ and the full-order-model (FOM) solution $U^{\text{FOM}} := U_{kh}^{\text{FOM}} \in X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\text{FOM}})$ with $V_h^{\text{ROM}} \subset V_h^{\text{FOM}} =: V_h$ is controlled by employing a dual-weighted residual method [8,11,12]. To this end, we obtain an optimization problem in which the model error between FOM and ROM, both measured in terms of some goal functional $J(\cdot)$, shall be minimized:

$$J(U^{\text{FOM}}) - J(U^{\text{ROM}}) \qquad (5)$$

subject to the constraint that the variational formulation of the poroelasticity problem (2) is satisfied by $U^{\text{FOM}}$ i.e., $A(U^{\text{FOM}})(\Phi) = F(\Phi)$ for all test functions $\Phi \in X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\text{FOM}})$. For more information on space-time error control, we refer the reader to [64,66,70,79]. We focus in this work on the enrichment of the reduced basis depending on the temporal evolution of the quantities of interest. The reduced basis is refined in a goal-oriented way to accurately and efficiently compute the solution over the whole temporal domain. In principle coarsening would also be possible, but is not the objective in this work. For coarsening, we would need to follow the work of Meyer and Matthies [52].

## Space-time dual-weighted residual method

The two Lagrange functionals for the constrained optimization problem (5) are defined as

$$
\begin{align*}
\mathcal{L}_{\square} : & X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\square}) \times X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\square}) \rightarrow \mathbb{R}, \\
& (U^{\square}, Z^{\square}) \mapsto J(U^{\square}) - A(U^{\square})(Z^{\square}) + F(Z^{\square})
\end{align*}
$$

with $\square \in \{\text{FOM, ROM}\}$. The stationary points $(U^{\text{FOM}}, Z^{\text{FOM}})$ and $(U^{\text{ROM}}, Z^{\text{ROM}})$ of the Lagrange functionals $\mathcal{L}_{\text{FOM}}$ and $\mathcal{L}_{\text{ROM}}$ need to satisfy the Karush-Kuhn-Tucker first-order optimality conditions.

## Primal problem

Firstly, the stationary points are solution to the following primal problems

$$
\begin{align*}
\mathcal{L}'_{\square,Z} (U^{\square}, Z^{\square}) (\delta Z^{\square}) &= -A(U^{\square}) (\delta Z^{\square}) + F(\delta Z^{\square}) = 0 \\
\forall \delta Z^{\square} \in X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\square}), \quad \square \in \{\text{FOM, ROM}\}.
\end{align*}
$$
Remark 2.3 To be consistent with the backward Euler time discretization in Remark 2.2, we again get rid of the timestep size k in displacement related terms and obtain

$$c(z_m^p - z_{m+1}^p, \phi^p) + \alpha(z_m^p - z_{m+1}^p, \nabla_x \cdot \phi^u) + k \frac{K}{\nu} (\nabla_x z_m^p, \nabla_x \phi^p) \\ + (\nabla_x z_m^u, \sigma(\phi^u)) - \alpha(\nabla_x z_m^u, \phi^p I) + \alpha \langle z_m^u, \phi^p n \rangle_{\Gamma_N} = J'_{U}(U)(\Phi) \\ \forall \Phi = (\phi^u, \phi^p) \in V(\Omega).$$

## MORe DWR: model order reduction with dual-weighted residual error estimates

The MORe DWR approach aims at evaluating the quantities of interest using a reduced-order model, while guaranteeing the error due to this approximation. Thus, the difference between the reduced-order-model (ROM) solution $U^{\text{ROM}} := U_{kh}^{\text{ROM}} \in X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\text{ROM}})$ and the full-order-model (FOM) solution $U^{\text{FOM}} := U_{kh}^{\text{FOM}} \in X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\text{FOM}})$ with $V_h^{\text{ROM}} \subset V_h^{\text{FOM}} =: V_h$ is controlled by employing a dual-weighted residual method [8,11,12]. To this end, we obtain an optimization problem in which the model error between FOM and ROM, both measured in terms of some goal functional $J(\cdot)$, shall be minimized:

$$J(U^{\text{FOM}}) - J(U^{\text{ROM}}) \qquad (5)$$

subject to the constraint that the variational formulation of the poroelasticity problem (2) is satisfied by $U^{\text{FOM}}$ i.e., $A(U^{\text{FOM}})(\Phi) = F(\Phi)$ for all test functions $\Phi \in X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\text{FOM}})$. For more information on space-time error control, we refer the reader to [64,66,70,79]. We focus in this work on the enrichment of the reduced basis depending on the temporal evolution of the quantities of interest. The reduced basis is refined in a goal-oriented way to accurately and efficiently compute the solution over the whole temporal domain. In principle coarsening would also be possible, but is not the objective in this work. For coarsening, we would need to follow the work of Meyer and Matthies [52].

## Space-time dual-weighted residual method

The two Lagrange functionals for the constrained optimization problem (5) are defined as

$$\mathcal{L}_{\square, Z}: X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\square}) \times X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\square}) \to \mathbb{R}, \\ (U^{\square}, Z^{\square}) \mapsto J(U^{\square}) - A(U^{\square})(Z^{\square}) + F(Z^{\square})$$

with $\square \in \{\text{FOM, ROM}\}$. The stationary points $(U^{\text{FOM}}, Z^{\text{FOM}})$ and $(U^{\text{ROM}}, Z^{\text{ROM}})$ of the Lagrange functionals $\mathcal{L}_{\text{FOM}}$ and $\mathcal{L}_{\text{ROM}}$ need to satisfy the Karush-Kuhn-Tucker first-order optimality conditions.

## Primal problem

Firstly, the stationary points are solution to the following primal problems

$$\mathcal{L}'_{\square,Z}(U^{\square}, Z^{\square})(\delta Z^{\square}) = -A(U^{\square})(\delta Z^{\square}) + F(\delta Z^{\square}) = 0 \\ \forall \delta Z^{\square} \in X_k^{\text{dG}(0)}(\mathcal{T}_k, V_h^{\square}), \quad \square \in \{\text{FOM, ROM}\}.$$
The corresponding solutions denoted $U^{\text{FOM}}$ and $U^{\text{ROM}}$ are the primal solutions. We observe that the primal solution can be obtained by solving the original problem, i.e., the quasi-static poroelasticity problem, forward in time as presented in Sect. 2.3.

### Adjoint problem

Secondly, the stationary points must also satisfy the following adjoint (or dual) problems

$$ 
\begin{aligned}
\mathcal{L}'_{\square,U}(U^{\square}, Z^{\square})(\delta U^{\square}) &= J'_{U}(U^{\square})(\delta U^{\square}) - A(\delta U^{\square})(Z^{\square}) = 0 \\
\forall \delta U^{\square} \in X_k^{\text{dG}(0)}(T_k, V_h^{\square}), \quad \square \in \{\text{FOM}, \text{ROM}\},
\end{aligned}
$$

which give the adjoint solutions $Z^{\text{FOM}}$ and $Z^{\text{ROM}}$. Note that here we already use that the poroelasticity equations are linear, which ensures that

$$ A'_{U}(U)(\delta U, Z) = A(\delta U)(Z). $$

Hence, the adjoint solution is obtained by solving

$$ A(\delta U)(Z) = J'_{U}(U)(\delta U), \qquad (6) $$

which is backward in time as discussed in Sect. 2.5.

**Remark 3.1** For linear goal functionals, the right-hand side of the adjoint problem (6) reduces to

$$ J'_{U}(U)(\delta U) = J(\delta U) $$

and the adjoint problem does not depend on the primal solution anymore.

### Error estimator

As a first result, we have

**Proposition 3.1** The computable ROM error estimator is given by

$$ J(U^{\text{FOM}}) \approx \eta := -A(U^{\text{ROM}})(Z^{\text{FOM}}) + F(Z^{\text{ROM}}). \qquad (7) $$

*Proof* The proof follows from Theorem 4.2 from [32]. It holds

$$ 
\begin{aligned}
J(U^{\text{FOM}}) - J(U^{\text{ROM}}) &= -A(U^{\text{ROM}})(Z^{\text{FOM}}) + F(Z^{\text{FOM}}) \\
&= -A(U^{\text{ROM}})(Z^{\text{ROM}}) + F(Z^{\text{ROM}}) \\
&\quad -A(U^{\text{ROM}})(Z^{\text{FOM}} - Z^{\text{ROM}}) + F(Z^{\text{FOM}} - Z^{\text{ROM}}) \\
&\approx -A(U^{\text{ROM}})(Z^{\text{ROM}}) + F(Z^{\text{ROM}}) =: \eta.
\end{aligned}
\qquad (8) 
$$

□

The previous result is assembled on each temporal element $I_m$ separately to localize the error in time. More concretely, for poroelasticity the temporally-localized error estimator reads:
The corresponding solutions denoted $U^{\text{FOM}}$ and $U^{\text{ROM}}$ are the primal solutions. We observe that the primal solution can be obtained by solving the original problem, i.e., the quasi-static poroelasticity problem, forward in time as presented in Sect. 2.3.

### Adjoint problem

Secondly, the stationary points must also satisfy the following adjoint (or dual) problems

$$ 
\begin{aligned}
\mathcal{L}'_{\square,U}(U^{\square}, Z^{\square})(\delta U^{\square}) &= J'_{U}(U^{\square})(\delta U^{\square}) - A(\delta U^{\square})(Z^{\square}) = 0 \\
\forall \delta U^{\square} \in X_k^{\text{dG}(0)}(T_k, V_h^{\square}), \quad \square \in \{\text{FOM}, \text{ROM}\},
\end{aligned}
$$

which give the adjoint solutions $Z^{\text{FOM}}$ and $Z^{\text{ROM}}$. Note that here we already use that the poroelasticity equations are linear, which ensures that

$$ A'_{U}(U)(\delta U, Z) = A(\delta U)(Z). $$

Hence, the adjoint solution is obtained by solving

$$ A(\delta U)(Z) = J'_{U}(U)(\delta U), \qquad (6) $$

which is backward in time as discussed in Sect. 2.5.

**Remark 3.1** For linear goal functionals, the right-hand side of the adjoint problem (6) reduces to

$$ J'_{U}(U)(\delta U) = J(\delta U) $$

and the adjoint problem does not depend on the primal solution anymore.

### Error estimator

As a first result, we have

**Proposition 3.1** The computable ROM error estimator is given by

$$ J(U^{\text{FOM}}) \approx \eta := -A(U^{\text{ROM}})(Z^{\text{FOM}}) + F(Z^{\text{ROM}}). \qquad (7) $$

*Proof* The proof follows from Theorem 4.2 from [32]. It holds

$$ 
\begin{aligned}
J(U^{\text{FOM}}) - J(U^{\text{ROM}}) &= -A(U^{\text{ROM}})(Z^{\text{FOM}}) + F(Z^{\text{FOM}}) \\
&= -A(U^{\text{ROM}})(Z^{\text{ROM}}) + F(Z^{\text{ROM}}) \\
&\quad -A(U^{\text{ROM}})(Z^{\text{FOM}} - Z^{\text{ROM}}) + F(Z^{\text{FOM}} - Z^{\text{ROM}}) \\
&\approx -A(U^{\text{ROM}})(Z^{\text{ROM}}) + F(Z^{\text{ROM}}) =: \eta.
\end{aligned}
\qquad (8) 
$$

□

The previous result is assembled on each temporal element $I_m$ separately to localize the error in time. More concretely, for poroelasticity the temporally-localized error estimator reads:
**Proposition 3.2** Given the primal problem (2), the goal functional (3), and the corresponding adjoint problem (6) with the adjoint bilinear form (4), and η from Proposition 3.1, the temporally-localized error estimator is given by

$$ \eta_m := \eta|_{I_m} = -c(p_m - p_{m-1}, z_{m-1}^p) - \alpha(\nabla_x \cdot u_m - \nabla_x \cdot u_{m-1}, z_{m-1}^p) \\ - k \cdot \frac{K}{\nu}(\nabla_x p_m, \nabla_x z_{m-1}^p) - (\sigma(u_m), \nabla_x z_{m-1}^u) \\ + \alpha(p_m I, \nabla_x z_{m-1}^u) - \alpha \langle p_m n, z_{m-1}^u \rangle_{\Gamma_N} + \langle t, z_{m-1}^u \rangle_{\Gamma_N}, $$

for $m = 1, \dots, M$. Moreover, from (7), we have

$$ J(U^{\text{FOM}}) - J(U^{\text{ROM}}) \approx \eta := \sum_{m=1}^{M} \eta_m. $$

The superscripts indicating reduced-order-model solutions for both the primal and dual solutions are omitted for a clearer presentation. The (reciprocal) effectivity index [6], which is the ratio between the estimated and the true errors, i.e.,

$$ I_{\text{eff}} := \left| \frac{J(U^{\text{FOM}}) - J(U^{\text{ROM}})}{\eta} \right| \qquad (9) $$

is used to measure the quality of our error estimator. We desire $I_{\text{eff}} \approx 1$ since then the error estimator can reliably predict the reduced-order-modeling error. We refer to [31] for two-sided proofs of discretization errors measured in goal functionals using saturation assumptions. We also observe that the effectivity index $I_{\text{eff}}$ is close to one in the numerical tests in Sect. 4. Finally, the quality of the adaptive enrichment is measured in terms of the indicator index [62]:

$$ I_{\text{ind}} := \frac{|J(U^{\text{FOM}}) - J(U^{\text{ROM}})|}{\sum_{m=1}^{M} |\eta_m|}. \qquad (10) $$

## Error estimator-based ROM updates

In this section, we give a compact summary of our novel approach of a goal-oriented incremental reduced-order model from [32] and discuss its extension to porous media settings. In the MORE DWR method, we marry a reduced-order model with a dual-weighted-residual-based error estimator and an incremental version of the POD algorithm. The goal is then to use the error estimator to identify when the solution behavior is not captured accurately by the reduced basis, such that it is incrementally enriched on-the-fly with new full-order-model snapshots. The approach is particularly cheap because the FOM snapshots are computed only for one temporal element.

In more detail, we apply our findings on error control of Sect. 3.1 to a backward Euler reduced-order model of poroelasticity. Further, an incremental basis generation is mandatory for the method to reduce computational operations and thus to be fast, which we realize by means of an incremental SVD. The incremental SVD is presented in Sect. 3.3.1. In this context, we also introduce the incremental POD as a trimmed version of the incremental SVD. Subsequently, the overall MORE DWR framework is depicted in Sect. 3.4, where all the ingredients are assembled and the final algorithm is presented.
**Proposition 3.2** Given the primal problem (2), the goal functional (3), and the corresponding adjoint problem (6) with the adjoint bilinear form (4), and η from Proposition 3.1, the temporally-localized error estimator is given by

$$ \eta_m := \eta|_{I_m} = -c(p_m - p_{m-1}, z_{m-1}^p) - \alpha(\nabla_x \cdot u_m - \nabla_x \cdot u_{m-1}, z_{m-1}^p) \\ - k \cdot \frac{K}{\nu}(\nabla_x p_m, \nabla_x z_{m-1}^p) - (\sigma(u_m), \nabla_x z_{m-1}^u) \\ + \alpha(p_m I, \nabla_x z_{m-1}^u) - \alpha \langle p_m n, z_{m-1}^u \rangle_{\Gamma_N} + \langle t, z_{m-1}^u \rangle_{\Gamma_N}, $$

for $m = 1, \dots, M$. Moreover, from (7), we have

$$ J(U^{\text{FOM}}) - J(U^{\text{ROM}}) \approx \eta := \sum_{m=1}^{M} \eta_m. $$

The superscripts indicating reduced-order-model solutions for both the primal and dual solutions are omitted for a clearer presentation. The (reciprocal) effectivity index [6], which is the ratio between the estimated and the true errors, i.e.,

$$ I_{\text{eff}} := \left| \frac{J(U^{\text{FOM}}) - J(U^{\text{ROM}})}{\eta} \right| \qquad (9) $$

is used to measure the quality of our error estimator. We desire $I_{\text{eff}} \approx 1$ since then the error estimator can reliably predict the reduced-order-modeling error. We refer to [31] for two-sided proofs of discretization errors measured in goal functionals using saturation assumptions. We also observe that the effectivity index $I_{\text{eff}}$ is close to one in the numerical tests in Sect. 4. Finally, the quality of the adaptive enrichment is measured in terms of the indicator index [62]:

$$ I_{\text{ind}} := \frac{|J(U^{\text{FOM}}) - J(U^{\text{ROM}})|}{\sum_{m=1}^{M} |\eta_m|}. \qquad (10) $$

### Error estimator-based ROM updates

In this section, we give a compact summary of our novel approach of a goal-oriented incremental reduced-order model from [32] and discuss its extension to porous media settings. In the MORe DWR method, we marry a reduced-order model with a dual-weighted-residual-based error estimator and an incremental version of the POD algorithm. The goal is then to use the error estimator to identify when the solution behavior is not captured accurately by the reduced basis, such that it is incrementally enriched on-the-fly with new full-order-model snapshots. The approach is particularly cheap because the FOM snapshots are computed only for one temporal element.

In more detail, we apply our findings on error control of Sect. 3.1 to a backward Euler reduced-order model of poroelasticity. Further, an incremental basis generation is mandatory for the method to reduce computational operations and thus to be fast, which we realize by means of an incremental SVD. The incremental SVD is presented in Sect. 3.3.1. In this context, we also introduce the incremental POD as a trimmed version of the incremental SVD. Subsequently, the overall MORe DWR framework is depicted in Sect. 3.4, where all the ingredients are assembled and the final algorithm is presented.
In summary, our novel approach avoids a computationally heavy offline phase and directly solves the reduced model. Our approach can be thought of as a replacement for traditional full-order-model simulations, since we highly reduce the computational cost, while controlling the error between the full-order (FOM) and the reduced-order model (ROM).

### Incremental proper orthogonal decomposition

The proper orthogonal decomposition (POD) [7, 15, 25, 27, 30, 39, 45, 47, 61, 68, 80] can be employed to project the differential equation onto a lower-dimensional solution manifold. In traditional applications of POD, one first collects full-order solution snapshots

$$U(t_0), U(t_1), \dots, U(t_M),$$

and finds the POD basis by computing the singular value decomposition

$$Y := [U(t_0) U(t_1) \dots U(t_M)] = \Psi \Sigma \Phi^T$$

with orthogonal matrices $\Psi \in \mathbb{R}^{n \times d}$, $\Phi \in \mathbb{R}^{(M+1) \times d}$ and $\Sigma = \text{diag}(\sigma_1, \dots, \sigma_d) \in \mathbb{R}^{d \times d}$. Here, $n$ is the number of spatial degrees of freedom, $M + 1$ is the number of FOM snapshots, and $d \le \min(n, M + 1)$. The SVD is then being truncated to the size $N \ll d$ using a retained energy criterion, cf. [36, 38, 47] e.g.,

$$\varepsilon(N) = \frac{\sum_{i=1}^{N} \sigma_i^2}{\sum_{i=1}^{d} \sigma_i^2} = \frac{\sum_{i=1}^{N} \sigma_i^2}{\sum_{i=1}^{M+1} \|U_i\|_{\mathbb{R}^n}^2} < 99.9\%. \quad (11)$$

For coupled problems with several PDEs, like poroelasticity, we could again perform the SVD of the snapshot matrix $Y$ consisting of the entire solution vectors $U(t_i) = \begin{pmatrix} u(t_i) \\ p(t_i) \end{pmatrix}$, but we choose to create separate SVDs for displacement $u(t_i)$ and pressure $p(t_i)$ such that POD bases with various sizes can be used for the different solution components.

Instead of building repeatedly a reduced basis from scratch, we suggest updating an already existing truncated SVD (tSVD) or solely its left-singular (POD) vectors according to modifications of the snapshot matrix without recomputing the whole tSVD or requiring access to the snapshot matrix [32,46]. The POD becomes incremental by appending additional snapshots to the initial snapshot matrix. This difference between the classical SVD and incremental SVD (iSVD) is illustrated in Fig. 1. In this context, we rely on the general approach of an additive rank-$b$ modification of the SVD, mainly developed by [23,24] and applied to the model-order reduction of fluid flows in [46].

Let $Y \in \mathbb{R}^{n \times m}$ be a given snapshot matrix that includes $\tilde{m} > 0$ snapshots. Usually, $\tilde{m}$ is equal or connected to the number of already computed time steps. Further, we have the rank-$N$ tSVD $\Psi \Sigma \Phi^T$ of the matrix $Y$. Additionally, let $b \in \mathbb{N}$ newly computed snapshots $\{U_1, \dots, U_b\}$ be stored in the bunch matrix

$$B = [u_1 \dots u_b] \in \mathbb{R}^{n \times b}. \quad (12)$$

We aim to compute the tSVD that is updated by the information contained in the bunch matrix $B$ according to

$$\tilde{\Psi} \tilde{\Sigma} \tilde{\Phi}^T = \tilde{Y} = [Y B]$$
In summary, our novel approach avoids a computationally heavy offline phase and directly solves the reduced model. Our approach can be thought of as a replacement for traditional full-order-model simulations, since we highly reduce the computational cost, while controlling the error between the full-order (FOM) and the reduced-order model (ROM).

### Incremental proper orthogonal decomposition

The proper orthogonal decomposition (POD) [7, 15, 25, 27, 30, 39, 45, 47, 61, 68, 80] can be employed to project the differential equation onto a lower-dimensional solution manifold. In traditional applications of POD, one first collects full-order solution snapshots

$$U(t_0), U(t_1), \dots, U(t_M),$$

and finds the POD basis by computing the singular value decomposition

$$Y := [U(t_0) U(t_1) \dots U(t_M)] = \Psi \Sigma \Phi^T$$

with orthogonal matrices $\Psi \in \mathbb{R}^{n \times d}$, $\Phi \in \mathbb{R}^{(M+1) \times d}$ and $\Sigma = \text{diag}(\sigma_1, \dots, \sigma_d) \in \mathbb{R}^{d \times d}$. Here, $n$ is the number of spatial degrees of freedom, $M + 1$ is the number of FOM snapshots, and $d \le \min(n, M + 1)$. The SVD is then being truncated to the size $N \ll d$ using a retained energy criterion, cf. [36, 38, 47] e.g.,

$$\varepsilon(N) = \frac{\sum_{i=1}^{N} \sigma_i^2}{\sum_{i=1}^{d} \sigma_i^2} = \frac{\sum_{i=1}^{N} \sigma_i^2}{\sum_{i=1}^{M+1} \|U_i\|_{\mathbb{R}^n}^2} < 99.9\%. \quad (11)$$

For coupled problems with several PDEs, like poroelasticity, we could again perform the SVD of the snapshot matrix $Y$ consisting of the entire solution vectors $U(t_i) = \begin{pmatrix} u(t_i) \\ p(t_i) \end{pmatrix}$, but we choose to create separate SVDs for displacement $u(t_i)$ and pressure $p(t_i)$ such that POD bases with various sizes can be used for the different solution components.

Instead of building repeatedly a reduced basis from scratch, we suggest updating an already existing truncated SVD (tSVD) or solely its left-singular (POD) vectors according to modifications of the snapshot matrix without recomputing the whole tSVD or requiring access to the snapshot matrix [32, 46]. The POD becomes incremental by appending additional snapshots to the initial snapshot matrix. This difference between the classical SVD and incremental SVD (iSVD) is illustrated in Fig. 1. In this context, we rely on the general approach of an additive rank-$b$ modification of the SVD, mainly developed by [23, 24] and applied to the model-order reduction of fluid flows in [46].

Let $Y \in \mathbb{R}^{n \times m}$ be a given snapshot matrix that includes $\tilde{m} > 0$ snapshots. Usually, $\tilde{m}$ is equal or connected to the number of already computed time steps. Further, we have the rank-$N$ tSVD $\Psi \Sigma \Phi^T$ of the matrix $Y$. Additionally, let $b \in \mathbb{N}$ newly computed snapshots $\{U_1, \dots, U_b\}$ be stored in the bunch matrix

$$B = [u_1 \dots u_b] \in \mathbb{R}^{n \times b}. \quad (12)$$

We aim to compute the tSVD that is updated by the information contained in the bunch matrix $B$ according to

$$\tilde{\Psi} \tilde{\Sigma} \tilde{\Phi}^T = \tilde{Y} = [Y B]$$
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Methodology of SVD (top) and incremental SVD (bottom)</figcaption></figure>

without explicitly recomputing $Y$ or $\tilde{Y}$ due to performance and memory reasons which was the original motivation of Brand's work on the incremental SVD, cf. [23, 24]. We skip the technical derivations of the incremental SVD in this work and instead summarize the incremental POD routine in Algorithm 1. A prototypical implementation of the incremental POD in Python can be found at https://github.com/Hendrik240298/Incremental_POD and for a mathematical deep dive into the theory behind the incremental POD, we refer the reader to [46][Section 2.2] and [32][Section 4.2.1].

**Algorithm 1 Incremental POD update**

**Input:** Reduced basis matrix $\Psi_N \in \mathbb{R}^{n \times N}$, singular value vector $\Sigma = [\sigma_1, \dots, \sigma_N] \in \mathbb{R}^N$, bunch matrix $B \in \mathbb{R}^{n \times b}$, and energy threshold $\varepsilon \in [0, 1]$.

**Output:** Reduced basis matrix $\tilde{\Psi}_N \in \mathbb{R}^{n \times \tilde{N}}$, singular value vector $\tilde{\Sigma} = [\tilde{\sigma}_1, \dots, \tilde{\sigma}_{\tilde{N}}]$

1: $H = \Psi_N^T B$

2: $P = B - \Psi_N H$

3: $[Q_P, R_P] = QR(P)$

4: $Q = [\Psi_N Q_P]$

5: $F = \begin{bmatrix} \text{diag}(\Sigma) & H \\ 0 & R_P \end{bmatrix}$

6: **if** Q not orthogonal **then**

7: $[Q, R] = QR(Q)$

8: $F = RF$

9: $[\Psi', \Sigma'] = SVD(F)$

10: $\tilde{N} = \min \left\{ N \in \mathbb{N} \mid \varepsilon(N) \ge \varepsilon, 1 \le N \le d \right\}$

11: $\tilde{\Sigma} = \text{diag}(\Sigma')(1 : \tilde{N})$

12: $\tilde{\Psi}_N = Q\Psi'(\cdot, 1 : \tilde{N})$

**Goal-oriented error-controlled incremental ROM**

This section outlines how goal-oriented error control (Sect. 3.1) drives the incremental reduced-order model (Sect. 3.3.1).

**The MORe DWR algorithm**

The overall procedure is illustrated in Fig. 2.
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Methodology of SVD (top) and incremental SVD (bottom)</figcaption></figure>

without explicitly recomputing $Y$ or $\tilde{Y}$ due to performance and memory reasons which was the original motivation of Brand's work on the incremental SVD, cf. [23, 24]. We skip the technical derivations of the incremental SVD in this work and instead summarize the incremental POD routine in Algorithm 1. A prototypical implementation of the incremental POD in Python can be found at https://github.com/Hendrik240298/Incremental_POD and for a mathematical deep dive into the theory behind the incremental POD, we refer the reader to [46][Section 2.2] and [32][Section 4.2.1].

**Algorithm 1 Incremental POD update**

**Input:** Reduced basis matrix $\Psi_N \in \mathbb{R}^{n \times N}$, singular value vector $\Sigma = [\sigma_1, \dots, \sigma_N] \in \mathbb{R}^N$, bunch matrix $B \in \mathbb{R}^{n \times b}$, and energy threshold $\varepsilon \in [0, 1]$.

**Output:** Reduced basis matrix $\tilde{\Psi}_N \in \mathbb{R}^{n \times \tilde{N}}$, singular value vector $\tilde{\Sigma} = [\tilde{\sigma}_1, \dots, \tilde{\sigma}_{\tilde{N}}]$

1: $H = \Psi_N^T B$

2: $P = B - \Psi_N H$

3: $[Q_P, R_P] = QR(P)$

4: $Q = [\Psi_N Q_P]$

5: $F = \begin{bmatrix} \text{diag}(\Sigma) H \\ 0 & R_P \end{bmatrix}$

6: **if** Q not orthogonal **then**

7: $[Q, R] = QR(Q)$

8: $F = RF$

9: $[\Psi', \Sigma'] = SVD(F)$

10: $\tilde{N} = \min \left\{ N \in \mathbb{N} \mid \varepsilon(N) \ge \varepsilon, 1 \le N \le d \right\}$

11: $\tilde{\Sigma} = \text{diag}(\Sigma')(1 : \tilde{N})$

12: $\tilde{\Psi}_N = Q\Psi'(:, 1 : \tilde{N})$

**Goal-oriented error-controlled incremental ROM**

This section outlines how goal-oriented error control (Sect. 3.1) drives the incremental reduced-order model (Sect. 3.3.1).

**The MORe DWR algorithm**

The overall procedure is illustrated in Fig. 2.
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Schematic representation of MORE DWR algorithm</figcaption></figure>

The aim is to solve the reduced-order model and adaptively enrich the reduced basis by means of the iPOD with full-order solutions until a given estimated error tolerance is reached for the chosen goal functional. The approach is designed to work without any prior knowledge or exploration of the solution manifold, while also attempting to minimize the number of full-order solves. The error associated with the reduced solution is estimated on each temporal element. The new full-order solution used for the basis enrichment is computed on the temporal element with the largest error. So, two full-order solves are conducted for each enrichment. We observe that this procedure is perfectly compatible with the adaptive basis selection based on DWR estimates presented by Meyer and Matthies in [52] to reduce the dimension of the reduced space. Thus, if incorporated it would be possible to either enrich or downsize the reduced basis according to the problem statement. The resulting approach is detailed in Algorithm 2. The reduced primal and dual solutions are used to provide a fast error estimation. Thus, the incremental ROM has both a primal and a dual reduced basis as input in the form of the reduced basis matrix.

**Initialization**

If no prior information is available, the bases are initialized by computing the primal and dual solution snapshots on the first temporal element. If prior reduced bases are available e.g., from a previous simulation on a different parameter configuration or some (very cheap) coarse-grid solution in the context of a multigrid idea, prior bases are used as an initial guess and are then altered by the goal-oriented adaptation. Thus, the MORE DWR approach is totally appropriate for reduced-order modeling of parameterized problems. The method can be an efficient substitute for the full-order model, as it builds a reduced basis tailored to a quantity of interest with a minimum of FOM solves.

Nevertheless, in the scenario without prior knowledge, we observed an issue when exclusively the first dual FOM solution ($t = T$) is utilized for the basis generation. We identified a significant discrepancy between the error estimates and the real errors, which is demonstrated and further described in Fig. 8. To mitigate this issue, further assessments utilizing the final dual snapshots for the basis enrichment at the first temporal elements ($t = 0$) were conducted, which consistently improved the accuracy of the error estimates.
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Schematic representation of MORE DWR algorithm</figcaption></figure>

The aim is to solve the reduced-order model and adaptively enrich the reduced basis by means of the iPOD with full-order solutions until a given estimated error tolerance is reached for the chosen goal functional. The approach is designed to work without any prior knowledge or exploration of the solution manifold, while also attempting to minimize the number of full-order solves. The error associated with the reduced solution is estimated on each temporal element. The new full-order solution used for the basis enrichment is computed on the temporal element with the largest error. So, two full-order solves are conducted for each enrichment. We observe that this procedure is perfectly compatible with the adaptive basis selection based on DWR estimates presented by Meyer and Matthies in [52] to reduce the dimension of the reduced space. Thus, if incorporated it would be possible to either enrich or downsize the reduced basis according to the problem statement. The resulting approach is detailed in Algorithm 2. The reduced primal and dual solutions are used to provide a fast error estimation. Thus, the incremental ROM has both a primal and a dual reduced basis as input in the form of the reduced basis matrix.

**Initialization**

If no prior information is available, the bases are initialized by computing the primal and dual solution snapshots on the first temporal element. If prior reduced bases are available e.g., from a previous simulation on a different parameter configuration or some (very cheap) coarse-grid solution in the context of a multigrid idea, prior bases are used as an initial guess and are then altered by the goal-oriented adaptation. Thus, the MORE DWR approach is totally appropriate for reduced-order modeling of parameterized problems. The method can be an efficient substitute for the full-order model, as it builds a reduced basis tailored to a quantity of interest with a minimum of FOM solves.

Nevertheless, in the scenario without prior knowledge, we observed an issue when exclusively the first dual FOM solution ($t = T$) is utilized for the basis generation. We identified a significant discrepancy between the error estimates and the real errors, which is demonstrated and further described in Fig. 8. To mitigate this issue, further assessments utilizing the final dual snapshots for the basis enrichment at the first temporal elements ($t = 0$) were conducted, which consistently improved the accuracy of the error estimates.
Nonetheless, a practical constraint is posed, as these snapshots are not readily accessible without completing the whole dual FOM solution process. Thus, a computationally efficient surrogate to bypass the need for the full-order solution process is required. For this, the dual FOM is substituted by the dual ROM, which is already computed during the enrichment process and only the last dual solution steps are conducted by the dual FOM. Given that the inconsistency in the estimate predominantly emerges during the initial enrichment iterations, we have optimized the process by confining the additional dual ROM enrichment solely to these early iterations for efficiency.

### Adaptive algorithm

At each MORe DWR iteration, the primal and dual reduced bases are enriched for the temporal element $m_{max}$ corresponding to the highest relative error comparing the local relative error

$$ \eta_m^{rel} = \frac{\eta_m}{J(U^{ROM}) + \sum_{l=1}^{M} \eta_l} \qquad (13) $$

associated with any temporal element $m$. Because the full-order solutions are not available, the error is normalized according to

$$ J(U^{FOM}) \approx J(U^{ROM}) + \eta = J(U^{ROM}) + \sum_{l=1}^{M} \eta_l. \qquad (14) $$

The adaptive MORe DWR algorithm is stopped according to the global relative error estimator

$$ \eta^{rel} = \frac{\eta}{J(U^{FOM})} \approx \frac{\sum_{m=1}^{M} \eta_m}{J(U^{ROM}) + \sum_{m=1}^{M} \eta_m} < TOL^{rel} \qquad (15) $$

with $TOL^{rel}$ a user-chosen threshold. This ensures that the relative error over the entire temporal domain is below a given tolerance, whereas in [32] we ensured that the local relative error is below a given tolerance. However, enforcing the accuracy locally is a more greedy approach, which leads to more full-order-model solves.

**Remark 3.2** The additional enrichment of the dual space to provide an accurate estimation of the error is omitted in Algorithm 2 and in Fig. 2 for the sake of clarity. This dual enrichment takes place after line 12 in Algorithm 2 and after enriching the primal and dual reduced bases in Fig. 2. The supplementary enrichment is based on the FOM solution of the dual problem for a few time steps, corresponding to the few last time steps of the dual problem i.e., the few first time steps of the primal problem. As the dual problem is backward in time, it starts from the ROM solution at $t = t_l$ with $l \in \mathbb{N}$ and the solution is computed for $t_{l-1}, t_{l-2}, \dots, t_0$. The dual reduced basis is then enriched with these additional snapshots. This procedure is applied only in the few first MORe DWR iterations prior to the update of the reduced system components and the error estimator.

### Numerical tests

The MORe DWR framework is numerically substantiated on two different problem configurations. The Mandel problem [1, 26, 29, 51, 73] (see also [33, 35, 40, 50, 78] and more recently also for nonlinear poroelasticity [73]) is first considered as a benchmark for two-dimensional poroelasticity and a linear goal functional. In Mandel's problem one
Nonetheless, a practical constraint is posed, as these snapshots are not readily accessible without completing the whole dual FOM solution process. Thus, a computationally efficient surrogate to bypass the need for the full-order solution process is required. For this, the dual FOM is substituted by the dual ROM, which is already computed during the enrichment process and only the last dual solution steps are conducted by the dual FOM. Given that the inconsistency in the estimate predominantly emerges during the initial enrichment iterations, we have optimized the process by confining the additional dual ROM enrichment solely to these early iterations for efficiency.

### Adaptive algorithm

At each MORe DWR iteration, the primal and dual reduced bases are enriched for the temporal element $m_{max}$ corresponding to the highest relative error comparing the local relative error

$$ \eta_m^{rel} = \frac{\eta_m}{J(U^{ROM}) + \sum_{l=1}^{M} \eta_l} \qquad (13) $$

associated with any temporal element $m$. Because the full-order solutions are not available, the error is normalized according to

$$ J(U^{FOM}) \approx J(U^{ROM}) + \eta = J(U^{ROM}) + \sum_{l=1}^{M} \eta_l. \qquad (14) $$

The adaptive MORe DWR algorithm is stopped according to the global relative error estimator

$$ \eta^{rel} = \frac{\eta}{J(U^{FOM})} \approx \frac{\sum_{m=1}^{M} \eta_m}{J(U^{ROM}) + \sum_{m=1}^{M} \eta_m} < TOL^{rel} \qquad (15) $$

with $TOL^{rel}$ a user-chosen threshold. This ensures that the relative error over the entire temporal domain is below a given tolerance, whereas in [32] we ensured that the local relative error is below a given tolerance. However, enforcing the accuracy locally is a more greedy approach, which leads to more full-order-model solves.

**Remark 3.2** The additional enrichment of the dual space to provide an accurate estimation of the error is omitted in Algorithm 2 and in Fig. 2 for the sake of clarity. This dual enrichment takes place after line 12 in Algorithm 2 and after enriching the primal and dual reduced bases in Fig. 2. The supplementary enrichment is based on the FOM solution of the dual problem for a few time steps, corresponding to the few last time steps of the dual problem i.e., the few first time steps of the primal problem. As the dual problem is backward in time, it starts from the ROM solution at $t = t_l$ with $l \in \mathbb{N}$ and the solution is computed for $t_{l-1}, t_{l-2}, \dots, t_0$. The dual reduced basis is then enriched with these additional snapshots. This procedure is applied only in the few first MORe DWR iterations prior to the update of the reduced system components and the error estimator.

### Numerical tests

The MORe DWR framework is numerically substantiated on two different problem configurations. The Mandel problem [1, 26, 29, 51, 73] (see also [33, 35, 40, 50, 78] and more recently also for nonlinear poroelasticity [73]) is first considered as a benchmark for two-dimensional poroelasticity and a linear goal functional. In Mandel's problem one
**Algorithm 2** MORe DWR algorithm

**Input:** Initial condition $U_0 := U(t_0)$, primal and dual reduced basis matrices $(\Psi_{N_p^u}^{p_u}, \Psi_{N_p^p}^{p_p})$ and $(\Psi_{N_d^u}^{d_u}, \Psi_{N_d^p}^{d_p})$, energy iPOD threshold $\epsilon \in [0, 1]$, and error tolerance tol > 0.

**Output:** Primal and dual reduced basis matrices $(\Psi_{N_p^u}^{p_u}, \Psi_{N_p^p}^{p_p})$ and $(\Psi_{N_d^u}^{d_u}, \Psi_{N_d^p}^{d_p})$ and reduced primal solutions $U_m$ for all $1 \le m \le M$.

1. **while** $\eta^{\text{rel}} > \text{tol do}$

2. **for** $m = 1, 2, \dots, M$ **do**

▷ Primal ROM

3. Solve reduced primal system on $I_m$ for $U_m^{\text{ROM}}$

4. **for** $m = M, M-1, \dots, 1$ **do**

▷ Dual ROM

5. Solve reduced dual system on $I_m$ for $Z_{m-1}^{\text{ROM}}$

6. Compute error estimate: $\eta^{\text{rel}}$

7. **if** $\eta^{\text{rel}} > \text{tol then}$

8. Find temporal element with maximum error: $m_{\text{max}} = \arg \max_{1 \le m \le L} \left| \frac{\eta_m}{J(U^{\text{ROM}}) + \sum_{l=1}^{M} \eta_l} \right|$

9. Solve primal full-order system on $I_{m_{\text{max}}}$ for $U_{m_{\text{max}}}^{\text{FOM}} = (u_{m_{\text{max}}}^{\text{FOM}}, p_{m_{\text{max}}}^{\text{FOM}})$

10. Update primal reduced basis:

$$ \Psi_{N_p^u}^{p_u} = \text{iPOD}(\Psi_{N_p^u}^{p_u}, \Sigma_{N_p^u}, [u_{m_{\text{max}}}^{\text{FOM}}], \epsilon) $$

$$ \Psi_{N_p^p}^{p_p} = \text{iPOD}(\Psi_{N_p^p}^{p_p}, \Sigma_{N_p^p}, [p_{m_{\text{max}}}^{\text{FOM}}], \epsilon) $$

11. Solve dual full-order system on $I_{m_{\text{max}}}$ for $Z_{m_{\text{max}}-1}^{\text{ROM}} = (z_{m_{\text{max}}-1}^{u, \text{FOM}}, z_{m_{\text{max}}-1}^{p, \text{FOM}})$

12. Update dual reduced basis:

$$ \Psi_{N_d^u}^{d_u} = \text{iPOD}(\Psi_{N_d^u}^{d_u}, \Sigma_{N_d^u}, [z_{m_{\text{max}}-1}^{u, \text{FOM}}], \epsilon) $$

$$ \Psi_{N_d^p}^{d_p} = \text{iPOD}(\Psi_{N_d^p}^{d_p}, \Sigma_{N_d^p}, [z_{m_{\text{max}}-1}^{p, \text{FOM}}], \epsilon) $$

13. Update reduced system components and error estimator

observes the so-called Mandel-Cryer effect [29,72] of a non-monotonic pressure evolution: first increasing pressure, followed by decreasing pressure. The second numerical test is a three-dimensional footing problem inspired by [34]. Our computations have been performed on an AMD EPYC 7H12 64-Core Processor. The FEM codes have been written in FEniCS [4] and the reduced-order modeling has been performed with NumPy [41] and SciPy [74].

**Mandel problem in 2D**

Let $\Omega := (0 \text{ m}, 100 \text{ m}) \times (0 \text{ m}, 20 \text{ m})$ and $I := (0 \text{ s}, 5\,000\,000 \text{ s}) \approx (0,58 \text{ d})$ with boundaries as shown in Fig. 3.

The initial and boundary conditions are given by

$$ p(\cdot, 0) = p^0(\cdot) = 0 \quad \text{in } \Omega \times \{0\}, $$

$$ u(\cdot, 0) = u^0(\cdot) = 0 \quad \text{in } \Omega \times \{0\}, $$

$$ \frac{K}{\nu} \nabla_x p \cdot n = 0 \quad \text{on } \partial\Omega \setminus \Gamma_{\text{right}} \times I, \text{ (No flow condition, homogeneous Neumann)} $$

$$ \sigma(u) \cdot n = -\bar{t}e_y \quad \text{on } \Gamma_{\text{top}} \times I, \text{ (Traction condition, inhomogeneous Neumann)} $$

$$ p = 0 \quad \text{on } \Gamma_{\text{right}} \times I, \text{ (Constant zero pressure, homogeneous Dirichlet)} $$

$$ \sigma(u) \cdot n = 0 \quad \text{on } \Gamma_{\text{right}} \times I, \text{ (Traction free, homogeneous Neumann)} $$

$$ u_y = 0 \text{ and } \partial_y u_x = 0 \quad \text{on } \Gamma_{\text{bottom}} \times I, \text{ (Confined conditions, mixed Dirichlet/Neumann)} $$
**Algorithm 2** MORe DWR algorithm

**Input:** Initial condition $U_0 := U(t_0)$, primal and dual reduced basis matrices $(\Psi_{N_p^u}^{p_u}, \Psi_{N_p^p}^{p_p})$ and $(\Psi_{N_d^u}^{d_u}, \Psi_{N_d^p}^{d_p})$, energy iPOD threshold $\epsilon \in [0, 1]$, and error tolerance $\text{tol} > 0$.

**Output:** Primal and dual reduced basis matrices $(\Psi_{N_p^u}^{p_u}, \Psi_{N_p^p}^{p_p})$ and $(\Psi_{N_d^u}^{d_u}, \Psi_{N_d^p}^{d_p})$ and reduced primal solutions $U_m$ for all $1 \le m \le M$.

1: **while** $\eta^{\text{rel}} > \text{tol do}$

2: **for** $m = 1, 2, \dots, M$ **do**

▷ Primal ROM

3: Solve reduced primal system on $I_m$ for $U_m^{\text{ROM}}$

4: **for** $m = M, M-1, \dots, 1$ **do**

▷ Dual ROM

5: Solve reduced dual system on $I_m$ for $Z_{m-1}^{\text{ROM}}$

6: Compute error estimate: $\eta^{\text{rel}}$

7: **if** $\eta^{\text{rel}} > \text{tol then}$

8: Find temporal element with maximum error: $m_{\text{max}} = \arg\max_{1 \le m \le L} \left| \frac{\eta_m}{J(U^{\text{ROM}}) + \sum_{l=1}^{M} \eta_l} \right|$

9: Solve primal full-order system on $I_{m_{\text{max}}}$ for $U_{m_{\text{max}}}^{\text{FOM}} = (u_{m_{\text{max}}}^{\text{FOM}}, p_{m_{\text{max}}}^{\text{FOM}})$

10: Update primal reduced basis:

$$ \Psi_{N_p^u}^{p_u} = \text{iPOD}(\Psi_{N_p^u}^{p_u}, \Sigma_{N_p^u}, [u_{m_{\text{max}}}^{\text{FOM}}], \epsilon) $$

$$ \Psi_{N_p^p}^{p_p} = \text{iPOD}(\Psi_{N_p^p}^{p_p}, \Sigma_{N_p^p}, [p_{m_{\text{max}}}^{\text{FOM}}], \epsilon) $$

11: Solve dual full-order system on $I_{m_{\text{max}}}$ for $Z_{m_{\text{max}}-1}^{\text{ROM}} = (z_{m_{\text{max}}-1}^{u, \text{FOM}}, z_{m_{\text{max}}-1}^{p, \text{FOM}})$

12: Update dual reduced basis:

$$ \Psi_{N_d^u}^{d_u} = \text{iPOD}(\Psi_{N_d^u}^{d_u}, \Sigma_{N_d^u}, [z_{m_{\text{max}}-1}^{u, \text{FOM}}], \epsilon) $$

$$ \Psi_{N_d^p}^{d_p} = \text{iPOD}(\Psi_{N_d^p}^{d_p}, \Sigma_{N_d^p}, [z_{m_{\text{max}}-1}^{p, \text{FOM}}], \epsilon) $$

13: Update reduced system components and error estimator

observes the so-called Mandel-Cryer effect [29,72] of a non-monotonic pressure evolution: first increasing pressure, followed by decreasing pressure. The second numerical test is a three-dimensional footing problem inspired by [34]. Our computations have been performed on an AMD EPYC 7H12 64-Core Processor. The FEM codes have been written in FEniCS [4] and the reduced-order modeling has been performed with NumPy [41] and SciPy [74].

**Mandel problem in 2D**

Let $\Omega := (0 \text{ m}, 100 \text{ m}) \times (0 \text{ m}, 20 \text{ m})$ and $I := (0 \text{ s}, 5\,000\,000 \text{ s}) \approx (0,58 \text{ d})$ with boundaries as shown in Fig. 3.

The initial and boundary conditions are given by

$$ p(\cdot, 0) = p^0(\cdot) = 0 \quad \text{in } \Omega \times \{0\}, $$

$$ u(\cdot, 0) = u^0(\cdot) = 0 \quad \text{in } \Omega \times \{0\}, $$

$$ \frac{K}{\nu} \nabla_x p \cdot n = 0 \quad \text{on } \partial\Omega \setminus \Gamma_{\text{right}} \times I, \text{ (No flow condition, homogeneous Neumann)} $$

$$ \sigma(u) \cdot n = -\bar{t}e_y \quad \text{on } \Gamma_{\text{top}} \times I, \text{ (Traction condition, inhomogeneous Neumann)} $$

$$ p = 0 \quad \text{on } \Gamma_{\text{right}} \times I, \text{ (Constant zero pressure, homogeneous Dirichlet)} $$

$$ \sigma(u) \cdot n = 0 \quad \text{on } \Gamma_{\text{right}} \times I, \text{ (Traction free, homogeneous Neumann)} $$

$$ u_y = 0 \text{ and } \partial_y u_x = 0 \quad \text{on } \Gamma_{\text{bottom}} \times I, \text{ (Confined conditions, mixed Dirichlet/Neumann)} $$
<figure><img src="image_5.png" /><figcaption>Fig. 3 Domain for Mandel's problem</figcaption></figure>

<table><caption>**Table 1** Parameters in Mandel's problem</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>M</td><td>1.75 × 10<sup>7</sup> Pa</td></tr><tr><td>c</td><td>1/M</td></tr><tr><td>α</td><td>1 Pa m</td></tr><tr><td>ν</td><td>1 × 10<sup>-3</sup> m<sup>2</sup>s<sup>-1</sup></td></tr><tr><td>K</td><td>1 × 10<sup>-13</sup> m<sup>2</sup></td></tr><tr><td>ρ</td><td>1 kg m<sup>-3</sup></td></tr><tr><td>l̄</td><td>1 × 10<sup>7</sup> Pa m</td></tr><tr><td>μ</td><td>1 × 10<sup>8</sup></td></tr><tr><td>λ</td><td>2/3 × 10<sup>8</sup></td></tr></tbody></table>

$u_x = 0$ and $\partial_x u_y = 0$ on $\Gamma_{\text{left}} \times I$. (Confined conditions, mixed Dirichlet/Neumann)

This results in the spatial function spaces (then modified in Section 2.1) being given by

$$V_u := \{u \in H^1(\Omega)^2 \mid u_y = 0 \text{ on } \Gamma_{\text{bottom}}, u_x = 0 \text{ on } \Gamma_{\text{left}}\},$$

$$V_p := \{p \in H^1(\Omega) \mid p = 0 \text{ on } \Gamma_{\text{right}}\}.$$

The parameters for Mandel's problem are summarized in Table 1.

In space, we use Taylor-Hood elements, i.e., quadratic finite elements for the displacements $u$ and linear finite elements for the pressure $p$. For this numerical test, a fixed spatial mesh with 80 spatial cells in the $x$-direction and 16 spatial cells in the $y$-direction is designed, which leads to an isotropic mesh with 10,626 DoFs for displacement and 1,377 DoFs for pressure. For the temporal discretization, the end time is $T = 5\,000\,000$ s $\approx$ 58 days with 5,000 temporal elements i.e., the chosen time step size $k = 1000$ s. For the quantity of interest, we choose the time-integrated pressure acting on the bottom boundary i.e.,

$$J(U) := \int_I \int_{\Gamma_{\text{bottom}}} p \, dx \, dt.$$

Further, for the MORe DWR enrichment, we choose $1 - 10^{-7}$ and $1 - 10^{-11}$ for the primal displacement and pressure energy tolerances (11). Accordingly, the dual energy tolerances are both set to $1 - 10^{-9}$. As argued in Sect. 3.4, it is crucial to have a dual space large enough for estimating accurately the error. Thus, the reduced dual basis is enriched during the first MORe DWR loops from the snapshots of the full-order dual solution on the first time steps $[t_0, t_{\ell-1}]$ i.e., for the last time points of the dual problem which is
<figure><img src="image_6.png" /><figcaption>Fig. 3 Domain for Mandel's problem</figcaption></figure>

<table><caption>**Table 1** Parameters in Mandel's problem</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>M</td><td>1.75 × 10<sup>7</sup> Pa</td></tr><tr><td>c</td><td>1/M</td></tr><tr><td>α</td><td>1 Pa m</td></tr><tr><td>ν</td><td>1 × 10<sup>-3</sup> m<sup>2</sup>s<sup>-1</sup></td></tr><tr><td>K</td><td>1 × 10<sup>-13</sup> m<sup>2</sup></td></tr><tr><td>ρ</td><td>1 kg m<sup>-3</sup></td></tr><tr><td>l̄</td><td>1 × 10<sup>7</sup> Pa m</td></tr><tr><td>μ</td><td>1 × 10<sup>8</sup></td></tr><tr><td>λ</td><td>2/3 × 10<sup>8</sup></td></tr></tbody></table>

$u_x = 0$ and $\partial_x u_y = 0$ on $\Gamma_{\text{left}} \times I$. (Confined conditions, mixed Dirichlet/Neumann)

This results in the spatial function spaces (then modified in Section 2.1) being given by

$$V_u := \{u \in H^1(\Omega)^2 \mid u_y = 0 \text{ on } \Gamma_{\text{bottom}}, u_x = 0 \text{ on } \Gamma_{\text{left}}\},$$

$$V_p := \{p \in H^1(\Omega) \mid p = 0 \text{ on } \Gamma_{\text{right}}\}.$$

The parameters for Mandel's problem are summarized in Table 1.

In space, we use Taylor-Hood elements, i.e., quadratic finite elements for the displacements $u$ and linear finite elements for the pressure $p$. For this numerical test, a fixed spatial mesh with 80 spatial cells in the $x$-direction and 16 spatial cells in the $y$-direction is designed, which leads to an isotropic mesh with 10,626 DoFs for displacement and 1,377 DoFs for pressure. For the temporal discretization, the end time is $T = 5\,000\,000$ s $\approx$ 58 days with 5,000 temporal elements i.e., the chosen time step size $k = 1000$ s. For the quantity of interest, we choose the time-integrated pressure acting on the bottom boundary i.e.,

$$J(U) := \int_I \int_{\Gamma_{\text{bottom}}} p \, dx \, dt.$$

Further, for the MORe DWR enrichment, we choose $1 - 10^{-7}$ and $1 - 10^{-11}$ for the primal displacement and pressure energy tolerances (11). Accordingly, the dual energy tolerances are both set to $1 - 10^{-9}$. As argued in Sect. 3.4, it is crucial to have a dual space large enough for estimating accurately the error. Thus, the reduced dual basis is enriched during the first MORe DWR loops from the snapshots of the full-order dual solution on the first time steps $[t_0, t_{\ell-1}]$ i.e., for the last time points of the dual problem which is
<figure><img src="image_6.png" /><figcaption>(a) Evolution of the true and estimated relative errors throughout the MORe DWR iterations.</figcaption></figure>

<caption>(b) Evolution of the POD basis sizes throughout the MORe DWR iterations.</caption>

<caption>**Fig. 4** Adaptivity in the MORe DWR method of Mandel's problem</caption>

<table><caption>**Table 2** Performance of MORe DWR method for the Mandel problem depending on the tolerance in the goal functional</caption><thead><tr><th>TOL<sup>rel</sup> [%]</th><th>e<sup>rel</sup> [%]</th><th>Speedup</th><th>FOM solves</th><th>ROM size</th><th>l<sub>eff</sub></th><th>l<sub>ind</sub></th></tr></thead><tbody><tr><td>0.1</td><td>0.123</td><td>8.5</td><td>97</td><td>5 / 33 + 37 / 26</td><td>1.399</td><td>1.377</td></tr><tr><td>1</td><td>0.821</td><td>9.5</td><td>83</td><td>5 / 26 + 34 / 25</td><td>1.089</td><td>1.092</td></tr><tr><td>2</td><td>2.03</td><td>10.1</td><td>77</td><td>5 / 24 + 33 / 25</td><td>1.038</td><td>1.042</td></tr><tr><td>5</td><td>4.89</td><td>11.2</td><td>69</td><td>5 / 21 + 31 / 23</td><td>1.052</td><td>1.054</td></tr><tr><td>10</td><td>10.7</td><td>12.7</td><td>59</td><td>4 / 16 + 27 / 20</td><td>1.113</td><td>1.110</td></tr><tr><td>20</td><td>19.6</td><td>14.3</td><td>51</td><td>4 / 13 + 24 / 18</td><td>1.143</td><td>1.142</td></tr></tbody></table>

formulated backward in time. Here, it is enriched for the first five MORe DWR iterations based on the snapshots corresponding with $[t_0, t_4]$ in the primal temporal discretization.

First, the adaptive enrichments of the method are shown in Fig. 4. The true relative error between the reference full-order solution and the reduced-order solution is given by

$$e^{rel} := \frac{|J(U^{FOM}) - J(U^{ROM})|}{|J(U^{FOM})|}.$$

Here, we observe a clear decrease with the iterations of the MORe DWR algorithm; see Fig. 4a. In contrast to the original MORe DWR paper [32], the relative error is now much closer to or even slightly exceeds the error tolerance. This can be explained by the fact that we do not enforce the error tolerance on each temporal element, but only on the whole integrated time domain. The estimated relative error is very close to the true error after iteration 2. Thereafter, only slight deviations can be noticed locally. On the contrary, a large discrepancy can be observed between the true and estimated errors for the first iterations, where the estimate yields significant underestimations due to the inaccuracy of the dual reduced space. The target tolerance of 0.1% is reached after 36 iterations. For larger error tolerances, the algorithm would terminate sooner e.g., for 1% relative error, we need 29 iterations, for 5% error we need 22 iterations and for 20% error only 13 iterations are required (see Table 2). The enrichment of the POD basis is illustrated in Fig. 4b. The primal reduced spaces increased progressively and regularly with the iterations. We observe that the primal pressure solution requires significantly more modes than the primal displacement. This implies that the choice of separating the bases for displacement and pressure is beneficial. Interestingly, the dual displacement solution (37 modes) requires more POD modes than the primal displacement solution (5 modes),
<figure><img src="image_7.png" /><figcaption>(a) Evolution of the true and estimated relative errors throughout the MORe DWR iterations.</figcaption></figure>

<caption>(b) Evolution of the POD basis sizes throughout the MORe DWR iterations.</caption>

<caption>**Fig. 4** Adaptivity in the MORe DWR method of Mandel's problem</caption>

<table><caption>**Table 2** Performance of MORe DWR method for the Mandel problem depending on the tolerance in the goal functional</caption><thead><tr><th>TOL<sup>rel</sup> [%]</th><th>e<sup>rel</sup> [%]</th><th>Speedup</th><th>FOM solves</th><th>ROM size</th><th>l<sub>eff</sub></th><th>l<sub>ind</sub></th></tr></thead><tbody><tr><td>0.1</td><td>0.123</td><td>8.5</td><td>97</td><td>5 / 33 + 37 / 26</td><td>1.399</td><td>1.377</td></tr><tr><td>1</td><td>0.821</td><td>9.5</td><td>83</td><td>5 / 26 + 34 / 25</td><td>1.089</td><td>1.092</td></tr><tr><td>2</td><td>2.03</td><td>10.1</td><td>77</td><td>5 / 24 + 33 / 25</td><td>1.038</td><td>1.042</td></tr><tr><td>5</td><td>4.89</td><td>11.2</td><td>69</td><td>5 / 21 + 31 / 23</td><td>1.052</td><td>1.054</td></tr><tr><td>10</td><td>10.7</td><td>12.7</td><td>59</td><td>4 / 16 + 27 / 20</td><td>1.113</td><td>1.110</td></tr><tr><td>20</td><td>19.6</td><td>14.3</td><td>51</td><td>4 / 13 + 24 / 18</td><td>1.143</td><td>1.142</td></tr></tbody></table>

formulated backward in time. Here, it is enriched for the first five MORe DWR iterations based on the snapshots corresponding with $[t_0, t_4]$ in the primal temporal discretization.

First, the adaptive enrichments of the method are shown in Fig. 4. The true relative error between the reference full-order solution and the reduced-order solution is given by

$$e^{rel} := \frac{|J(U^{FOM}) - J(U^{ROM})|}{|J(U^{FOM})|}.$$

Here, we observe a clear decrease with the iterations of the MORe DWR algorithm; see Fig. 4a. In contrast to the original MORe DWR paper [32], the relative error is now much closer to or even slightly exceeds the error tolerance. This can be explained by the fact that we do not enforce the error tolerance on each temporal element, but only on the whole integrated time domain. The estimated relative error is very close to the true error after iteration 2. Thereafter, only slight deviations can be noticed locally. On the contrary, a large discrepancy can be observed between the true and estimated errors for the first iterations, where the estimate yields significant underestimations due to the inaccuracy of the dual reduced space. The target tolerance of 0.1% is reached after 36 iterations. For larger error tolerances, the algorithm would terminate sooner e.g., for 1% relative error, we need 29 iterations, for 5% error we need 22 iterations and for 20% error only 13 iterations are required (see Table 2). The enrichment of the POD basis is illustrated in Fig. 4b. The primal reduced spaces increased progressively and regularly with the iterations. We observe that the primal pressure solution requires significantly more modes than the primal displacement. This implies that the choice of separating the bases for displacement and pressure is beneficial. Interestingly, the dual displacement solution (37 modes) requires more POD modes than the primal displacement solution (5 modes),
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Goal functionals for Mandel's problem and different relative error tolerances</figcaption></figure>

whereas the dual pressure solution (26 modes) requires less POD modes than the primal pressure solution (33 modes). We also see that the sizes of the dual bases increase faster than in the primal case during the first 5 iterations as we enforce an additional dual basis enrichment. But, its slope decreases after these initial enrichments.

The time trajectories of the goal functional from the full-order solution $U_h$ and the reduced-order solution $U_N$ are compared for each temporal element in Fig. 5. For these comparisons, we choose relative error tolerances between 0.1% and 20%. Both quantity of interest trajectories are almost indistinguishable for tolerances between 0.1% and 2%. Some differences appear between the full-order and the reduced-order results for larger relative error tolerances. They might be seen as a too large error. However, for certain real-world applications an error of e.g., 10% might still be acceptable for an appealing computational cost. We can see that choosing the relative error tolerance, the MORe DWR approach allows to control the quality of the approximation of the true quantity of interest, even within a reduced-order solution framework.

Table 2 gives an overview of the results obtained for the Mandel problem using the MORe DWR algorithm with various error tolerances $TOL^{rel}$ between 0.1% and 20%. For each tolerance value, the relative true error $e^{rel}$ is given as well as the computational speedup, the total number of FOM solves, the POD basis sizes for the primal displacement/pressure and the dual displacement/pressure, the effectivity index $I_{eff}$ from (9) and the indicator index $I_{ind}$ defined in (10). Here, the number of FOM solves sums up all primal and dual solves and the basis sizes are shown in the pattern: primal displacement / primal pressure + dual displacement / dual pressure.
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Goal functionals for Mandel's problem and different relative error tolerances</figcaption></figure>

whereas the dual pressure solution (26 modes) requires less POD modes than the primal pressure solution (33 modes). We also see that the sizes of the dual bases increase faster than in the primal case during the first 5 iterations as we enforce an additional dual basis enrichment. But, its slope decreases after these initial enrichments.

The time trajectories of the goal functional from the full-order solution $U_h$ and the reduced-order solution $U_N$ are compared for each temporal element in Fig. 5. For these comparisons, we choose relative error tolerances between 0.1% and 20%. Both quantity of interest trajectories are almost indistinguishable for tolerances between 0.1% and 2%. Some differences appear between the full-order and the reduced-order results for larger relative error tolerances. They might be seen as a too large error. However, for certain real-world applications an error of e.g., 10% might still be acceptable for an appealing computational cost. We can see that choosing the relative error tolerance, the MORe DWR approach allows to control the quality of the approximation of the true quantity of interest, even within a reduced-order solution framework.

Table 2 gives an overview of the results obtained for the Mandel problem using the MORe DWR algorithm with various error tolerances $TOL^{rel}$ between 0.1% and 20%. For each tolerance value, the relative true error $e^{rel}$ is given as well as the computational speedup, the total number of FOM solves, the POD basis sizes for the primal displacement/pressure and the dual displacement/pressure, the effectivity index $I_{eff}$ from (9) and the indicator index $I_{ind}$ defined in (10). Here, the number of FOM solves sums up all primal and dual solves and the basis sizes are shown in the pattern: primal displacement / primal pressure + dual displacement / dual pressure.
<figure><img src="image_8.png" /><figcaption>Fig. 6 Domain for 3D footing problem</figcaption></figure>

The relative error as well as the speedup increase with a rise in tolerance. The relative error is rather close to the error tolerance, but sometimes slightly exceeds it. Indeed, here error estimates are employed, not error bounds. This can be explained by the fact that we do not enforce the error tolerance on each temporal element, but only on the whole integrated time domain. In contrast, enforcing the error tolerance on each temporal element leads to much lower error than the error tolerance [32]. The speedup is explained by the decreasing amount of FOM solves and smaller POD bases for both the primal and dual problems, as well as displacement and pressure w.r.t. the given tolerance. The effectivity index and the indicator index are close to 1 as expected for a linear PDE and goal functional. This validates the usage of dual-weighted residual error estimates to control the reduced-order-modeling error because the estimated error accurately approximates the true error.

## Footing problem in 3D

In this second numerical example, a three-dimensional footing problem inspired by [34] is studied. Let $\Omega := (-32 \text{ m}, 32 \text{ m}) \times (-32 \text{ m}, 32 \text{ m}) \times (0 \text{ m}, 64 \text{ m})$ and $I := (0 \text{ s}, 5 000 000 \text{ s})$ with boundaries as shown in Fig. 6.

The initial and boundary conditions are given by

$$p(0) = p^0 = 0 \quad \text{in } \Omega \times \{0\},$$

$$u(0) = u^0 = 0 \quad \text{in } \Omega \times \{0\},$$

$$\frac{K}{\nu} \nabla_x p \cdot n = 0 \quad \text{on } \partial\Omega \setminus \Gamma_{\text{bottom}} \times I, \quad (\text{No flow condition, homogeneous Neumann})$$

$$\sigma(u) \cdot n = -\bar{t}e_z \quad \text{on } \Gamma_{\text{compression}} \times I, \quad (\text{Traction condition, inhomogeneous Neumann})$$

$$\sigma(u) \cdot n = 0 \quad \text{on } \Gamma_{\text{top}} \setminus \Gamma_{\text{compression}} \times I, \quad (\text{Traction-free, homogeneous Neumann})$$

$$p = 0 \quad \text{on } \Gamma_{\text{bottom}} \times I, \quad (\text{Constant zero pressure, homogeneous Dirichlet})$$

$$u = 0 \quad \text{on } \Gamma_{\text{bottom}} \times I, \quad (\text{Fixed displacements, homogeneous Dirichlet})$$

$$\sigma(u) \cdot n = 0 \quad \text{on } \Gamma_{\text{wall}} \times I, \quad (\text{Traction-free, homogeneous Neumann})$$

This results in the spatial function spaces (then modified in Sect. 2.1) being given by

$$V_u := \{u \in H^1(\Omega)^3 \mid u = 0 \text{ on } \Gamma_{\text{bottom}}\},$$
<figure><img src="image_9.png" /><figcaption>Fig. 6 Domain for 3D footing problem</figcaption></figure>

The relative error as well as the speedup increase with a rise in tolerance. The relative error is rather close to the error tolerance, but sometimes slightly exceeds it. Indeed, here error estimates are employed, not error bounds. This can be explained by the fact that we do not enforce the error tolerance on each temporal element, but only on the whole integrated time domain. In contrast, enforcing the error tolerance on each temporal element leads to much lower error than the error tolerance [32]. The speedup is explained by the decreasing amount of FOM solves and smaller POD bases for both the primal and dual problems, as well as displacement and pressure w.r.t. the given tolerance. The effectivity index and the indicator index are close to 1 as expected for a linear PDE and goal functional. This validates the usage of dual-weighted residual error estimates to control the reduced-order-modeling error because the estimated error accurately approximates the true error.

## Footing problem in 3D

In this second numerical example, a three-dimensional footing problem inspired by [34] is studied. Let $\Omega := (-32 \text{ m}, 32 \text{ m}) \times (-32 \text{ m}, 32 \text{ m}) \times (0 \text{ m}, 64 \text{ m})$ and $I := (0 \text{ s}, 5000000 \text{ s})$ with boundaries as shown in Fig. 6.

The initial and boundary conditions are given by

$$p(0) = p^0 = 0 \quad \text{in } \Omega \times \{0\},$$

$$u(0) = u^0 = 0 \quad \text{in } \Omega \times \{0\},$$

$$\frac{K}{\nu} \nabla_x p \cdot n = 0 \quad \text{on } \partial\Omega \setminus \Gamma_{\text{bottom}} \times I, \quad (\text{No flow condition, homogeneous Neumann})$$

$$\sigma(u) \cdot n = -\bar{t}e_z \quad \text{on } \Gamma_{\text{compression}} \times I, \quad (\text{Traction condition, inhomogeneous Neumann})$$

$$\sigma(u) \cdot n = 0 \quad \text{on } \Gamma_{\text{top}} \setminus \Gamma_{\text{compression}} \times I, \quad (\text{Traction-free, homogeneous Neumann})$$

$$p = 0 \quad \text{on } \Gamma_{\text{bottom}} \times I, \quad (\text{Constant zero pressure, homogeneous Dirichlet})$$

$$u = 0 \quad \text{on } \Gamma_{\text{bottom}} \times I, \quad (\text{Fixed displacements, homogeneous Dirichlet})$$

$$\sigma(u) \cdot n = 0 \quad \text{on } \Gamma_{\text{wall}} \times I, \quad (\text{Traction-free, homogeneous Neumann})$$

This results in the spatial function spaces (then modified in Sect. 2.1) being given by

$$V_u := \{u \in H^1(\Omega)^3 \mid u = 0 \text{ on } \Gamma_{\text{bottom}}\},$$
$$V_p := \{p \in H^1(\Omega) | p = 0 \text{ on } \Gamma_{\text{bottom}}\}.$$

The material parameters as listed in Table 1 are used again.

As before, Taylor-Hood elements, namely quadratic finite elements in space for the displacement $u$ and linear finite elements for the pressure $p$, are employed. The spatial mesh is fixed, it comprises 16 spatial cells in each direction i.e., an isotropic mesh with 107, 811 DoFs for displacement and 4, 913 DoFs for pressure. The temporal domain from 0 to the end time $T = 5\,000\,000$ s $\approx 58$ d is discretized with 5, 000 temporal elements, which leads to a time step size $k = 1000$ s. For the goal functional, we choose the time-integrated pressure acting at the compression boundary i.e.,

$$J(U) := \int_I \int_{\Gamma_{\text{compression}}} p \, dx \, dt.$$

A linear system of equations with 112, 724 unknowns is obtained, which has to be solved for each temporal element. The previous 2D Mandel problem was solved with a direct solver. Due to the much larger size of the equation system, we now have to resort to an iterative solution scheme for maintaining a low computational cost. The equation system is not symmetric and the memory for the calculations is not limited; thus we choose the generalized minimal residual method (GMRES). Specifically, SciPy [74] implementation of GMRES is utilized with a tolerance for convergence of $5 \cdot 10^{-8}$. For the preconditioning of the system, we have tested a diagonal Jacobi preconditioner as well as the smoothed aggregation algebraic multigrid (SA-AMG) method utilizing the implementation provided by [13]. Both approaches have been benchmarked on the footing problem for the first 500 temporal elements. The SA-AMG ansatz needs a mean iteration count of 6.33 for convergence per temporal element solve, while the Jacobi preconditioner requires 283.81 iterations to reach the same accuracy. Although the number of iterations is significantly lower for the SA-AMG method, the mean wall time per temporal element solve is 93.37 s for the SA-AMG approach, but only 3.89 s for the Jacobi preconditioner. The Jacobi preconditioner is significantly faster than the SA-AMG method in our case, most likely due to easier parallelization capabilities on multicore machines. Thus, the Jacobi preconditioner is employed for this application case.

Similar to Mandel's problem, we choose $1 - 10^{-7}$ and $1 - 10^{-11}$ for the primal displacement and pressure energy tolerance (11). Again, the dual energy tolerances are both set to $1 - 10^{-9}$. As argued in Sect. 3.4, the dual spaces are enriched with additional dual full-order solutions. In this case, we perform this additional enrichment for the first 8 MORe DWR iterations and include the dual solutions for the last five time steps of the dual problem i.e., the time steps $[t_0, t_4]$ of the primal problem discretization.

The adaptive enrichment of the method is illustrated in Fig. 7. The relative error between the full-order and reduced-order solutions, as well as the error estimator is shown in Fig. 7a. We observe a similar behavior for the first iterations as in Mandel's problem: a discrepancy between the estimate and the error, with a severe underestimation of the error. Thereafter, both quantities align for the rest of the iterations. In addition, we see a sharp decline after 18 iterations that is followed by a non-monotonic behavior, highlighting the greedy enrichment process of the MORe DWR method.

The evolution of the size of the primal/dual POD bases for displacement/pressure is shown in Fig. 7b. We see again a steep increase in the dual POD sizes in the first iterations
$$V_p := \{p \in H^1(\Omega) | p = 0 \text{ on } \Gamma_{\text{bottom}}\}.$$

The material parameters as listed in Table 1 are used again.

As before, Taylor-Hood elements, namely quadratic finite elements in space for the displacement $u$ and linear finite elements for the pressure $p$, are employed. The spatial mesh is fixed, it comprises 16 spatial cells in each direction i.e., an isotropic mesh with 107, 811 DoFs for displacement and 4, 913 DoFs for pressure. The temporal domain from 0 to the end time $T = 5\,000\,000$ s $\approx 58$ d is discretized with 5, 000 temporal elements, which leads to a time step size $k = 1000$ s. For the goal functional, we choose the time-integrated pressure acting at the compression boundary i.e.,

$$J(U) := \int_I \int_{\Gamma_{\text{compression}}} p \, dx \, dt.$$

A linear system of equations with 112, 724 unknowns is obtained, which has to be solved for each temporal element. The previous 2D Mandel problem was solved with a direct solver. Due to the much larger size of the equation system, we now have to resort to an iterative solution scheme for maintaining a low computational cost. The equation system is not symmetric and the memory for the calculations is not limited; thus we choose the generalized minimal residual method (GMRES). Specifically, SciPy [74] implementation of GMRES is utilized with a tolerance for convergence of $5 \cdot 10^{-8}$. For the preconditioning of the system, we have tested a diagonal Jacobi preconditioner as well as the smoothed aggregation algebraic multigrid (SA-AMG) method utilizing the implementation provided by [13]. Both approaches have been benchmarked on the footing problem for the first 500 temporal elements. The SA-AMG ansatz needs a mean iteration count of 6.33 for convergence per temporal element solve, while the Jacobi preconditioner requires 283.81 iterations to reach the same accuracy. Although the number of iterations is significantly lower for the SA-AMG method, the mean wall time per temporal element solve is 93.37 s for the SA-AMG approach, but only 3.89 s for the Jacobi preconditioner. The Jacobi preconditioner is significantly faster than the SA-AMG method in our case, most likely due to easier parallelization capabilities on multicore machines. Thus, the Jacobi preconditioner is employed for this application case.

Similar to Mandel's problem, we choose $1 - 10^{-7}$ and $1 - 10^{-11}$ for the primal displacement and pressure energy tolerance (11). Again, the dual energy tolerances are both set to $1 - 10^{-9}$. As argued in Sect. 3.4, the dual spaces are enriched with additional dual full-order solutions. In this case, we perform this additional enrichment for the first 8 MORe DWR iterations and include the dual solutions for the last five time steps of the dual problem i.e., the time steps $[t_0, t_4]$ of the primal problem discretization.

The adaptive enrichment of the method is illustrated in Fig. 7. The relative error between the full-order and reduced-order solutions, as well as the error estimator is shown in Fig. 7a. We observe a similar behavior for the first iterations as in Mandel's problem: a discrepancy between the estimate and the error, with a severe underestimation of the error. Thereafter, both quantities align for the rest of the iterations. In addition, we see a sharp decline after 18 iterations that is followed by a non-monotonic behavior, highlighting the greedy enrichment process of the MORe DWR method.

The evolution of the size of the primal/dual POD bases for displacement/pressure is shown in Fig. 7b. We see again a steep increase in the dual POD sizes in the first iterations
<figure><img src="image_9.png" /><figcaption>(a) Evolution of the true and estimated relative errors throughout the MORe DWR iterations.</figcaption></figure>

<caption>(b) Evolution of the POD basis sizes throughout the MORe DWR iterations.</caption>

<figure><img src="image_10.png" /><figcaption>(a) Evolution of the true and estimated relative errors throughout the MORe DWR iterations.</figcaption></figure>

<caption>(b) Evolution of the POD basis sizes throughout the MORe DWR iterations.</caption>

<caption>**Fig. 8** Adaptivity in the MORe DWR method of the 3D footing problem *without additional dual basis enrichments*</caption>

that then flattens out and is finally surpassed by the primal pressure POD size. The size of the primal displacement POD size is almost constant, with a maximum of 4.

Next, in Fig. 8 we compare the latter results with the MORe DWR method *neglecting* the additional dual reduced space enrichment proposed in Sect. 3.4. A minimum number of 20 iterations is enforced to ensure the consistency of the results because it is expected from the previous study that the error is severely underestimated during the first iterations without extra enrichment (Fig. 7a). Figure 8a compares the relative error between the full-order and reduced-order solutions with its estimate over the course of MORe DWR iterations. In contrast to the results in Fig. 7a, we observe that the first phase of discrepancy between the estimate and the error, with the typical severe underestimation of the error, lasts much longer with more variation. The true error and its estimate align after 35 iterations. Thereafter, the error estimate is as accurate as when using the extra dual enrichment. Thus, the extra enrichment of the dual basis appears necessary for an accurate error estimate from the early stage of the algorithm. The ROM solution that complies to a tolerance of 0.5% is obtained with a total of 104 FOM evaluations (52 primal plus 52 dual evaluations) without the dual enrichment, whereas only 80 FOM evaluations (20 primal plus 20 dual plus 40 extra dual) were required with the enriched dual basis. Thus, although we added dual FOM evaluations in the first iterations of the extra-enriched approach, the trade-off seems beneficial, as this approach reduces the total number of iterations and thus the computational costs.
<figure><img src="image_10.png" /><figcaption>(a) Evolution of the true and estimated relative errors throughout the MORe DWR iterations.</figcaption></figure>

<caption>(b) Evolution of the POD basis sizes throughout the MORe DWR iterations.</caption>

<figure><img src="image_11.png" /><figcaption>(a) Evolution of the true and estimated relative errors throughout the MORe DWR iterations.</figcaption></figure>

<caption>(b) Evolution of the POD basis sizes throughout the MORe DWR iterations.</caption>

<caption>**Fig. 8** Adaptivity in the MORe DWR method of the 3D footing problem *without additional dual basis enrichments*</caption>

that then flattens out and is finally surpassed by the primal pressure POD size. The size of the primal displacement POD size is almost constant, with a maximum of 4.

Next, in Fig. 8 we compare the latter results with the MORe DWR method *neglecting* the additional dual reduced space enrichment proposed in Sect. 3.4. A minimum number of 20 iterations is enforced to ensure the consistency of the results because it is expected from the previous study that the error is severely underestimated during the first iterations without extra enrichment (Fig. 7a). Figure 8a compares the relative error between the full-order and reduced-order solutions with its estimate over the course of MORe DWR iterations. In contrast to the results in Fig. 7a, we observe that the first phase of discrepancy between the estimate and the error, with the typical severe underestimation of the error, lasts much longer with more variation. The true error and its estimate align after 35 iterations. Thereafter, the error estimate is as accurate as when using the extra dual enrichment. Thus, the extra enrichment of the dual basis appears necessary for an accurate error estimate from the early stage of the algorithm. The ROM solution that complies to a tolerance of 0.5% is obtained with a total of 104 FOM evaluations (52 primal plus 52 dual evaluations) without the dual enrichment, whereas only 80 FOM evaluations (20 primal plus 20 dual plus 40 extra dual) were required with the enriched dual basis. Thus, although we added dual FOM evaluations in the first iterations of the extra-enriched approach, the trade-off seems beneficial, as this approach reduces the total number of iterations and thus the computational costs.
<figure><img src="image_11.png" /><figcaption>**Fig. 9** Goal functionals for the 3D footing problem, depending on the tolerance in the goal functional</figcaption></figure>

In addition, the development of the POD bases in case of no additional dual enrichment is shown in Fig. 8b. In contrast to Fig. 7b, except for the primal displacement, all the POD bases start growing at nearly the same rate. The dual pressure basis size stagnates after 42 iterations, and the primal displacement has, after 7 iterations, a constant size of 4. Comparing the POD basis sizes for both approaches for a relative tolerance of 0.5%, the final dual POD basis sizes are nearly the same. In contrast, in the case of the extra enrichment, the primal pressure POD size of 19 is considerably smaller than the 47 modes needed when neglecting this enrichment.

The time trajectory of the goal functional from the reduced-order solution is compared with the trajectory given by the full-order solution for error tolerances between 0.1% and 20% in Fig. 9. The trajectories of quantities of interest are nearly identical for a tolerance of 0.1%. Differences emerge between the full-order and reduced-order goal functionals for larger tolerances. In essence, the reduced-order solution reaches the accuracy chosen by the user.

The time trajectories of the goal functional from the reduced-order solution is compared with the trajectory given by the full-order solution for error tolerances between 0.1% and 20% in Fig. 9. The trajectories of both quantities of interest are nearly identical for a tolerance of 0.1%. Differences emerge between the full-order and reduced-order goal functionals for larger tolerances. In essence, the reduced-order solution reaches the accuracy chosen by the user.

In Table 3, we give an overview of the simulation results of the three-dimensional footing problem for different error tolerances between 0.1% and 20%. Tolerances are chosen based on their potential applicability in real-world scenarios. Lower tolerances could necessitate too many POD modes for the reduced-order model, rendering it less
<figure><img src="image_12.png" /><figcaption>**Fig. 9** Goal functionals for the 3D footing problem, depending on the tolerance in the goal functional</figcaption></figure>

In addition, the development of the POD bases in case of no additional dual enrichment is shown in Fig. 8b. In contrast to Fig. 7b, except for the primal displacement, all the POD bases start growing at nearly the same rate. The dual pressure basis size stagnates after 42 iterations, and the primal displacement has, after 7 iterations, a constant size of 4. Comparing the POD basis sizes for both approaches for a relative tolerance of 0.5%, the final dual POD basis sizes are nearly the same. In contrast, in the case of the extra enrichment, the primal pressure POD size of 19 is considerably smaller than the 47 modes needed when neglecting this enrichment.

The time trajectory of the goal functional from the reduced-order solution is compared with the trajectory given by the full-order solution for error tolerances between 0.1% and 20% in Fig. 9. The trajectories of quantities of interest are nearly identical for a tolerance of 0.1%. Differences emerge between the full-order and reduced-order goal functionals for larger tolerances. In essence, the reduced-order solution reaches the accuracy chosen by the user.

The time trajectories of the goal functional from the reduced-order solution is compared with the trajectory given by the full-order solution for error tolerances between 0.1% and 20% in Fig. 9. The trajectories of both quantities of interest are nearly identical for a tolerance of 0.1%. Differences emerge between the full-order and reduced-order goal functionals for larger tolerances. In essence, the reduced-order solution reaches the accuracy chosen by the user.

In Table 3, we give an overview of the simulation results of the three-dimensional footing problem for different error tolerances between 0.1% and 20%. Tolerances are chosen based on their potential applicability in real-world scenarios. Lower tolerances could necessitate too many POD modes for the reduced-order model, rendering it less
<table><caption>**Table 3** Performance of MORe DWR method for the 3D footing problem, depending on the tolerance in the goal functional</caption><thead><tr><th>TOL<sup>rel</sup> [%]</th><th>e<sup>rel</sup> [%]</th><th>Speedup</th><th>FOM solves</th><th>ROM size</th><th>l<sub>eff</sub></th><th>l<sub>ind</sub></th></tr></thead><tbody><tr><td>0.1</td><td>0.0971</td><td>8.6</td><td>220</td><td>4 / 55 + 53 / 28</td><td>0.999</td><td>1.207</td></tr><tr><td>0.5</td><td>0.5333</td><td>21.2</td><td>80</td><td>4 / 19 + 38 / 27</td><td>1.068</td><td>3.441</td></tr><tr><td>1</td><td>0.579</td><td>22.4</td><td>78</td><td>4 / 18 + 38 / 26</td><td>1.084</td><td>3.378</td></tr><tr><td>2</td><td>0.579</td><td>21.7</td><td>78</td><td>4 / 18 + 38 / 26</td><td>1.084</td><td>3.378</td></tr><tr><td>5</td><td>0.579</td><td>22.2</td><td>78</td><td>4 / 18 + 38 / 26</td><td>1.084</td><td>3.378</td></tr><tr><td>10</td><td>8.49</td><td>22.4</td><td>76</td><td>4 / 17 + 38 / 26</td><td>1.008</td><td>1.099</td></tr><tr><td>20</td><td>19.9</td><td>26.2</td><td>66</td><td>3 / 13 + 33 / 24</td><td>1.005</td><td>1.031</td></tr></tbody></table>

efficient. The tolerances are met besides in the case of 0.5%, where the error is slightly underestimated. We observe that with a rise in the tolerance, the relative error as well as the speedup increase. We note that the algorithm has the same behavior for the tolerances 1%, 2% and 5%. Indeed, these tolerances are reached within the same number of MORe DWR iterations, yielding all other quantities to be the same. This behavior can be explained by the sharp decline in the error after 19 MORe DWR iterations, c.f. Figure 7a. The evaluation of the dual-weighted residual error estimates by means of the effectivity and indicator indices shows that the estimate integrated over the temporal domain is almost perfect, as the effectivity index is near 1.0 in all cases. However, the indicator indices vary from 1.00 up to 3.44 which translates to inaccurate temporal error localizations in the case of high indicator indices.

## Conclusion

We developed an adaptive incremental reduced-order (MORe DWR) model for single-phase flow problems in porous media, namely the Biot system. To this end, we further developed and extended our prior work [32]. Therein, the following ingredients are combined. First, a space-time formulation and Galerkin space-time finite element discretization for the coupled problem. Second, a ROM method based on an incremental POD. Third, a goal-oriented a posteriori error estimator using the dual-weighted residual method for adaptively refining the ROM bases. The efficiency of the methodology has been demonstrated, by using the incremental ROM for the two-dimensional Mandel benchmark and a three-dimensional footing problem. The algorithm is stopped according to the total relative error considering the whole temporal domain. The speed-ups offered by this approach compared with the reference full-order approach vary from 8.5x in two dimensions for the smaller relative tolerances, to 26.2x for the three-dimensional case with large relative tolerances. The accuracy of the error estimator has been improved by additional enrichment of the dual reduced bases. The true error and the error estimator almost coincide with effectivity indices close to one, which makes the error estimates reliable. It is important to note that the reduced-order model starts with no other prior information than a single primal snapshot and a dual FOM snapshot. Hence, the full-order model needs to be queried in the online stage to gather snapshots for basis enrichment. Consequently, our approach could be further sped up by efficient full-order model solvers. Therefore, for the three-dimensional test, we used GMRES as our iterative solver and compared a simple diagonal Jacobi preconditioner and an algebraic multigrid preconditioner. An aspect for future work is the combination of temporal adaptivity and reduced basis adaptivity. The
<table><caption>**Table 3** Performance of MORe DWR method for the 3D footing problem, depending on the tolerance in the goal functional</caption><thead><tr><th>TOL<sup>rel</sup> [%]</th><th>e<sup>rel</sup> [%]</th><th>Speedup</th><th>FOM solves</th><th>ROM size</th><th>l<sub>eff</sub></th><th>l<sub>ind</sub></th></tr></thead><tbody><tr><td>0.1</td><td>0.0971</td><td>8.6</td><td>220</td><td>4 / 55 + 53 / 28</td><td>0.999</td><td>1.207</td></tr><tr><td>0.5</td><td>0.5333</td><td>21.2</td><td>80</td><td>4 / 19 + 38 / 27</td><td>1.068</td><td>3.441</td></tr><tr><td>1</td><td>0.579</td><td>22.4</td><td>78</td><td>4 / 18 + 38 / 26</td><td>1.084</td><td>3.378</td></tr><tr><td>2</td><td>0.579</td><td>21.7</td><td>78</td><td>4 / 18 + 38 / 26</td><td>1.084</td><td>3.378</td></tr><tr><td>5</td><td>0.579</td><td>22.2</td><td>78</td><td>4 / 18 + 38 / 26</td><td>1.084</td><td>3.378</td></tr><tr><td>10</td><td>8.49</td><td>22.4</td><td>76</td><td>4 / 17 + 38 / 26</td><td>1.008</td><td>1.099</td></tr><tr><td>20</td><td>19.9</td><td>26.2</td><td>66</td><td>3 / 13 + 33 / 24</td><td>1.005</td><td>1.031</td></tr></tbody></table>

efficient. The tolerances are met besides in the case of 0.5%, where the error is slightly underestimated. We observe that with a rise in the tolerance, the relative error as well as the speedup increase. We note that the algorithm has the same behavior for the tolerances 1%, 2% and 5%. Indeed, these tolerances are reached within the same number of MORe DWR iterations, yielding all other quantities to be the same. This behavior can be explained by the sharp decline in the error after 19 MORe DWR iterations, c.f. Figure 7a. The evaluation of the dual-weighted residual error estimates by means of the effectivity and indicator indices shows that the estimate integrated over the temporal domain is almost perfect, as the effectivity index is near 1.0 in all cases. However, the indicator indices vary from 1.00 up to 3.44 which translates to inaccurate temporal error localizations in the case of high indicator indices.

## Conclusion

We developed an adaptive incremental reduced-order (MORe DWR) model for single-phase flow problems in porous media, namely the Biot system. To this end, we further developed and extended our prior work [32]. Therein, the following ingredients are combined. First, a space-time formulation and Galerkin space-time finite element discretization for the coupled problem. Second, a ROM method based on an incremental POD. Third, a goal-oriented a posteriori error estimator using the dual-weighted residual method for adaptively refining the ROM bases. The efficiency of the methodology has been demonstrated, by using the incremental ROM for the two-dimensional Mandel benchmark and a three-dimensional footing problem. The algorithm is stopped according to the total relative error considering the whole temporal domain. The speed-ups offered by this approach compared with the reference full-order approach vary from 8.5x in two dimensions for the smaller relative tolerances, to 26.2x for the three-dimensional case with large relative tolerances. The accuracy of the error estimator has been improved by additional enrichment of the dual reduced bases. The true error and the error estimator almost coincide with effectivity indices close to one, which makes the error estimates reliable. It is important to note that the reduced-order model starts with no other prior information than a single primal snapshot and a dual FOM snapshot. Hence, the full-order model needs to be queried in the online stage to gather snapshots for basis enrichment. Consequently, our approach could be further sped up by efficient full-order model solvers. Therefore, for the three-dimensional test, we used GMRES as our iterative solver and compared a simple diagonal Jacobi preconditioner and an algebraic multigrid preconditioner. An aspect for future work is the combination of temporal adaptivity and reduced basis adaptivity. The
MORE DWR algorithm already requires the solution of the primal and dual problems. This additional computational effort can be reused to refine the temporal mesh adaptively using the dual-weighted residual method and could reduce the number of time steps, thus further speeding up the computations.

### Acknowledgements

The authors acknowledge the funding of the German Research Foundation (DFG) within the framework of the International Research Training Group on Computational Mechanics Techniques in High Dimensions GRK 2657 under Grant Number 433082294. The support of the French-German University through the French-German Doctoral college "Sophisticated Numerical and Testing Approaches" (CDFA-DFDK 19-04) is also acknowledged. Moreover, T. Wick gratefully acknowledges the J. Tinsley Oden Faculty Fellowship Research Program at UT Austin for the visit in July/August 2023.

### Author contributions

H. Fischer: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing—original draft. Writing—review & editing, Visualization. J. Roth: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing—original draft, Writing—review & editing, Visualization. L. Chamoin: Conceptualization, Formal analysis, Writing—review & editing, Supervision, Funding acquisition. A. Fau: Conceptualization, Formal analysis, Writing—review & editing, Supervision, Funding acquisition. M. Wheeler: Conceptualization, Formal analysis, Discussions, Writing—review & editing. T. Wick: Conceptualization, Formal analysis, Resources, Writing—review & editing, Supervision, Funding acquisition.

### Funding Information

Open Access funding enabled and organized by Projekt DEAL.

### Data Availability

The incremental POD in Python can be found at https://github.com/Hendrik240298/Incremental_POD. All other algorithms and data will be made available on reasonable request.

### Declarations

#### Competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Received: 15 November 2023 Accepted: 19 March 2024

Published online: 18 April 2024

### References

1. Abousleiman YN, Cheng A, Cui L, Detournay E, Roegiers JC. Mandel's problem revisited. Geotechnique. 1996;46:187-95.

2. Alameddin S, Fau A, Néron D, Ladevèze P, Nackenhorst U. Toward optimality of proper generalised decomposition bases. Math Comput Appl. 2019;24(1).

3. Almani T, Kumar K, Dogru A, Singh G, Wheeler M. Convergence analysis of multirate fixed-stress split iterative schemes for coupling flow with geomechanics. Comput Methods Appl Mech Eng. 2016;311:180-207.

4. Alnæs M, Blechta J, Hake J, Johansson A, Kehlet B, Logg A, Richardson C, Ring J, Rognes M E, Wells G N. The FEniCS Project Version 1.5. Arch Numer Softw. 2015;3.

5. Anselmann M, Bause M, Margenberg N, Shamko P. An energy-efficient gmres-multigrid solver for space-time finite element computation of dynamic poro- and thermoelasticity; 2023.

6. Babuska I, Rheinboldt WC. A-posteriori error estimates for the finite element method. Int J Numer Meth Eng. 1978;12(10):1597-615.

7. Baiges J, Codina R, Idelsohn S. Explicit reduced-order models for the stabilized finite element approximation of the incompressible Navier-Stokes equations. Int J Numer Methods Fluids. 2013;72(12):1219-43.

8. Bangerth W, Rannacher R. Adaptive finite element methods for differential equations. Birkhäuser Verlag, 2003.

9. Bastian P, Helmig R. Efficient fully-coupled solution techniques for two-phase flow in porous media: parallel multigrid solution and large scale computations. Adv Water Resour. 1999;23(3):199-216.

10. Bause M, Radu F, Köcher U. Space-time finite element approximation of the biot poroelasticity system with iterative coupling. Comput Methods Appl Mech Energy. 2017;320:745-68.

11. Becker R, Rannacher R. A feed-back approach to error control in finite element methods: basic analysis and examples. East-West J Numer Math. 1996;4:237-64.

12. Becker R, Rannacher R. An optimal control approach to a posteriori error estimation in finite element methods. Acta Numer. 2001;10:1-102.

13. Bell N, Olson LN, Schroder J, Southworth B. PyAMG: algebraic multigrid solvers in python. J Open Source Softw. 2023;8(87):5495.

14. Benner P, Cohen A, Ohlberger M, Willcox K. Model reduction and approximation: theory and algorithms. SIAM Philadelphia, 2015.

15. Benner P, Schilders W, Grivet-Talocia S, Quarteroni A, Rozza G, Miguel Silveira L. Model Order Reduction: Volume 2: Snapshot-Based Methods and Algorithms. De Gruyter, 2020.
MORE DWR algorithm already requires the solution of the primal and dual problems. This additional computational effort can be reused to refine the temporal mesh adaptively using the dual-weighted residual method and could reduce the number of time steps, thus further speeding up the computations.

**Acknowledgements**

The authors acknowledge the funding of the German Research Foundation (DFG) within the framework of the International Research Training Group on Computational Mechanics Techniques in High Dimensions GRK 2657 under Grant Number 433082294. The support of the French-German University through the French-German Doctoral college "Sophisticated Numerical and Testing Approaches" (CDFA-DFDK 19-04) is also acknowledged. Moreover, T. Wick gratefully acknowledges the J. Tinsley Oden Faculty Fellowship Research Program at UT Austin for the visit in July/August 2023.

**Author contributions**

H. Fischer: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing—original draft. Writing—review & editing, Visualization. J. Roth: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing—original draft, Writing—review & editing, Visualization. L. Chamoin: Conceptualization, Formal analysis, Writing—review & editing, Supervision, Funding acquisition. A. Fau: Conceptualization, Formal analysis, Writing—review & editing, Supervision, Funding acquisition. M. Wheeler: Conceptualization, Formal analysis, Discussions, Writing—review & editing. T. Wick: Conceptualization, Formal analysis, Resources, Writing—review & editing, Supervision, Funding acquisition.

**Funding Information**

Open Access funding enabled and organized by Projekt DEAL.

**Data Availability**

The incremental POD in Python can be found at https://github.com/Hendrik240298/Incremental_POD. All other algorithms and data will be made available on reasonable request.

**Declarations**

**Competing interests**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Received: 15 November 2023 Accepted: 19 March 2024

Published online: 18 April 2024

**References**

1. Abousleiman YN, Cheng A, Cui L, Detournay E, Roegiers JC. Mandel's problem revisited. Geotechnique. 1996;46:187-95.

2. Alameddin S, Fau A, Néron D, Ladevèze P, Nackenhorst U. Toward optimality of proper generalised decomposition bases. Math Comput Appl. 2019;24(1).

3. Almani T, Kumar K, Dogru A, Singh G, Wheeler M. Convergence analysis of multirate fixed-stress split iterative schemes for coupling flow with geomechanics. Comput Methods Appl Mech Eng. 2016;311:180-207.

4. Alnæs M, Blechta J, Hake J, Johansson A, Kehlet B, Logg A, Richardson C, Ring J, Rognes M E, Wells G N. The FEniCS Project Version 1.5. Arch Numer Softw. 2015;3.

5. Anselmann M, Bause M, Margenberg N, Shamko P. An energy-efficient gmres-multigrid solver for space-time finite element computation of dynamic poro- and thermoelasticity; 2023.

6. Babuska I, Rheinboldt WC. A-posteriori error estimates for the finite element method. Int J Numer Meth Eng. 1978;12(10):1597-615.

7. Baiges J, Codina R, Idelsohn S. Explicit reduced-order models for the stabilized finite element approximation of the incompressible Navier-Stokes equations. Int J Numer Methods Fluids. 2013;72(12):1219-43.

8. Bangerth W, Rannacher R. Adaptive finite element methods for differential equations. Birkhäuser Verlag, 2003.

9. Bastian P, Helmig R. Efficient fully-coupled solution techniques for two-phase flow in porous media: parallel multigrid solution and large scale computations. Adv Water Resour. 1999;23(3):199-216.

10. Bause M, Radu F, Köcher U. Space-time finite element approximation of the biot poroelasticity system with iterative coupling. Comput Methods Appl Mech Energy. 2017;320:745-68.

11. Becker R, Rannacher R. A feed-back approach to error control in finite element methods: basic analysis and examples. East-West J Numer Math. 1996;4:237-64.

12. Becker R, Rannacher R. An optimal control approach to a posteriori error estimation in finite element methods. Acta Numer. 2001;10:1-102.

13. Bell N, Olson LN, Schroder J, Southworth B. PyAMG: algebraic multigrid solvers in python. J Open Source Softw. 2023;8(87):5495.

14. Benner P, Cohen A, Ohlberger M, Willcox K. Model reduction and approximation: theory and algorithms. SIAM Philadelphia, 2015.

15. Benner P, Schilders W, Grivet-Talocia S, Quarteroni A, Rozza G, Miguel Silveira L. Model Order Reduction: Volume 2: Snapshot-Based Methods and Algorithms. De Gruyter, 2020.
16. Bhattacharyya M, Fau A, Desmorat R, Alameddin S, Néron D, Ladevèze P, Nackenhorst U. A kinetic two-scale damage model for high-cycle fatigue simulation using multi-temporal Latin framework. Eur J Mech A Solids. 2019;77: 103808.

17. Bhattacharyya M, Fau A, Nackenhorst U, Néron D, Ladevèze P. A Latin-based model reduction approach for the simulation of cycling damage. Comput Mech. 2018;62:725.

18. Biot M. Consolidation settlement under a rectangular load distribution. J Appl Phys. 1941;12(5):426-30.

19. Biot M. General theory of three-dimensional consolidation. J Appl Phys. 1941;12(2):155-64.

20. Biot M. Theory of elasticity and consolidation for a porous anisotropic solid. J Appl Phys. 1955;25:182-5.

21. Biot M. Theory of finite deformations of porous solids. Indiana Univ Math J. 1971/72;21:597-620.

22. Both J, Borregales M, Nordbotten J, Kumar K, Radu F. Robust fixed stress splitting for biot's equations in heterogeneous media. Appl Math Lett. 2017;68:101-8.

23. Brand M. Incremental Singular Value Decomposition of Uncertain Data with Missing Values. In: European Conference on Computer Vision, pages 707-720. Springer, 2002.

24. Brand M. Fast low-rank modifications of the thin singular value decomposition. Linear Algebra Appl. 2006;415(1):20-30.

25. Caiazzo A, Iliescu T, John V, Schyschlowa S. A numerical investigation of velocity-pressure reduced order models for incompressible flows. J Comput Phys. 2014;259:598-616.

26. Cheng AH-D. A direct boundary element method for plane strain poroelasticity. Int J Numer Anal Meth Geomech. 1988;12:551-72.

27. Christensen EA, Brøns M, Sørensen JN. Evaluation of proper orthogonal decomposition-based decomposition techniques applied to parameter-dependent nonturbulent flows. SIAM J Sci Comput. 1999;21(4):1419-34.

28. Coussy O. Poromechanics. Wiley; 2004.

29. Cryer C. A comparison of the three-dimensional consolidation theories of Biot and Terzaghi. Q J Mech Appl Math. 1663;16:401-12.

30. Ekre F, Larsson F, Runesson K, Jänicke R. A posteriori error estimation for numerical model reduction in computational homogenization of porous media. Int J Numer Methods Eng. 2020;121(23):5350-80.

31. Endtmayer B, Langer U, Wick T. Two-side a posteriori error estimates for the dual-weighted residual method. SIAM J Sci Comput. 2020;42(1):A371-94.

32. Fischer H, Roth J, Wick T, Chamoin L, Fau A. MORe DWR: space-time goal-oriented error control for incremental POD-based ROM for time-averaged goal functionals. J Comp Phys. 2024;504: 112863.

33. Gai X. A coupled geomechanics and reservoir flow model on parallel computers. Doctoral Thesis, University of Texas at Austin, 2004.

34. Gaspar FJ, Gracia JL, Lisbona FJ, Oosterlee CW. Distributive smoothers in multigrid for problems with dominating grad-div operators. Numer Linear Algebra Appl. 2008;15(8):661-83.

35. Girault V, Pencheva G, Wheeler M, Wildey T. Domain decomposition for Poroelasticity and Elasticity with DG jumps and mortars. Math Models Methods Appl Sci. 2011;21:10.

36. Gräßle C, Hinze M. POD reduced-order modeling for evolution equations utilizing arbitrary finite element discretizations. Adv Comput Math. 2018;44(6):1941-78.

37. Grepl MA, Patera AT. A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations. ESAIM: M2AN. 2005;39(1):157-81.

38. Gubisch M, Volkwein S. Chapter 1: proper orthogonal decomposition for linear-quadratic optimal control: Theory and algorithms. In: Model Reduction and Approximation, pages 3-63. Society for Industrial and Applied Mathematics, 2017.

39. Gunzburger MD, Peterson JS, Shadid JN. Reduced-order modeling of time-dependent PDEs with multiple parameters in the boundary data. Comput Methods Appl Mech Eng. 2007;196(4-6):1030-47.

40. Guzman H A F. Domain decomposition methods in geomechanics. Doctoral Thesis, University of Texas at Austin, 2012.

41. Harris CR, Millman KJ, van der Walt SJ, Gommers R, Virtanen P, Cournapeau D, Wieser E, Taylor J, Berg S, Smith NJ, Kern R, Picus M, Hoyer S, van Kerkwijk MH, Brett M, Haldane A, del Río JF, Wiebe M, Peterson P, Gérard-Marchant P, Sheppard K, Reddy T, Weckesser W, Abbasi H, Gohlke C, Oliphant TE. Array programming with NumPy. Nature. 2020;585(7825):357-62.

42. Hong Q, Kraus J, Lymbery M, Philo F. Conservative discretizations and parameter-robust preconditioners for biot and multiple-network flux-based poroelasticity models. Numer Linear Algebra Appl. 2019;26(4): e2242.

43. Kergrene K, Chamoin L, Laforest M, Prudhomme S. On a goal-oriented version of the proper generalized decomposition method. J Sci Comput. 2019;81(1):92-111.

44. Kim J, Tchelepi H, Juanes R. Stability and convergence of sequential methods for coupled flow and geomechanics: fixed-stress and fixed-strain splits. Comp Methods Appl Mech Energy. 2011;200(13-16):1591-606.

45. Kunisch K, Volkwein S. Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics. SIAM J Numer Anal. 2002;40(2):492-515.

46. Kühl N, Fischer H, Hinze M, Rung T. An incremental singular value decomposition approach for large-scale spatially parallel & distributed but temporally serial data-applied to technical flows. Comput Phys Commun. 2024;296: 109022.

47. Lassila T, Manzoni A, Quarteroni A, Rozza G. Model Order Reduction in Fluid Dynamics: Challenges and Perspectives. In: Reduced Order Methods for Modeling and Computational Reduction, volume 9, pages 235-273. Springer International Publishing, 2014.

48. Lee S, Lee Y-J, Wheeler MF. A locally conservative enriched Galerkin approximation and efficient solver for elliptic and parabolic problems. SIAM J Sci Comput. 2016;38(3):A1404-29.

49. Lewis RW, Schrefler B. The finite element method in the static and dynamic deformation and consolidation of porous media, 2nd Edition. Wiley. 1999.

50. Liu R. Discontinuous Galerkin Finite Element Solution for Poromechanics. PhD thesis, The University of Texas at Austin, 2004.

51. Mandel J. Consolidation des sols (Étude mathématique). Géotechnique. 1953;3(7):287-99.
16. Bhattacharyya M, Fau A, Desmorat R, Alameddin S, Néron D, Ladevèze P, Nackenhorst U. A kinetic two-scale damage model for high-cycle fatigue simulation using multi-temporal Latin framework. Eur J Mech A Solids. 2019;77: 103808.

17. Bhattacharyya M, Fau A, Nackenhorst U, Néron D, Ladevèze P. A Latin-based model reduction approach for the simulation of cycling damage. Comput Mech. 2018;62:725.

18. Biot M. Consolidation settlement under a rectangular load distribution. J Appl Phys. 1941;12(5):426-30.

19. Biot M. General theory of three-dimensional consolidation. J Appl Phys. 1941;12(2):155-64.

20. Biot M. Theory of elasticity and consolidation for a porous anisotropic solid. J Appl Phys. 1955;25:182-5.

21. Biot M. Theory of finite deformations of porous solids. Indiana Univ Math J. 1971/72;21:597-620.

22. Both J, Borregales M, Nordbotten J, Kumar K, Radu F. Robust fixed stress splitting for biot's equations in heterogeneous media. Appl Math Lett. 2017;68:101-8.

23. Brand M. Incremental Singular Value Decomposition of Uncertain Data with Missing Values. In: European Conference on Computer Vision, pages 707-720. Springer, 2002.

24. Brand M. Fast low-rank modifications of the thin singular value decomposition. Linear Algebra Appl. 2006;415(1):20-30.

25. Caiazzo A, Iliescu T, John V, Schyschlowa S. A numerical investigation of velocity-pressure reduced order models for incompressible flows. J Comput Phys. 2014;259:598-616.

26. Cheng AH-D. A direct boundary element method for plane strain poroelasticity. Int J Numer Anal Meth Geomech. 1988;12:551-72.

27. Christensen EA, Brøns M, Sørensen JN. Evaluation of proper orthogonal decomposition-based decomposition techniques applied to parameter-dependent nonturbulent flows. SIAM J Sci Comput. 1999;21(4):1419-34.

28. Coussy O. Poromechanics. Wiley; 2004.

29. Cryer C. A comparison of the three-dimensional consolidation theories of Biot and Terzaghi. Q J Mech Appl Math. 1663;16:401-12.

30. Ekre F, Larsson F, Runesson K, Jänicke R. A posteriori error estimation for numerical model reduction in computational homogenization of porous media. Int J Numer Methods Eng. 2020;121(23):5350-80.

31. Endtmayer B, Langer U, Wick T. Two-side a posteriori error estimates for the dual-weighted residual method. SIAM J Sci Comput. 2020;42(1):A371-94.

32. Fischer H, Roth J, Wick T, Chamoin L, Fau A. MORe DWR: space-time goal-oriented error control for incremental POD-based ROM for time-averaged goal functionals. J Comp Phys. 2024;504: 112863.

33. Gai X. A coupled geomechanics and reservoir flow model on parallel computers. Doctoral Thesis, University of Texas at Austin, 2004.

34. Gaspar FJ, Gracia JL, Lisbona FJ, Oosterlee CW. Distributive smoothers in multigrid for problems with dominating grad-div operators. Numer Linear Algebra Appl. 2008;15(8):661-83.

35. Girault V, Pencheva G, Wheeler M, Wildey T. Domain decomposition for Poroelasticity and Elasticity with DG jumps and mortars. Math Models Methods Appl Sci. 2011;21:10.

36. Gräßle C, Hinze M. POD reduced-order modeling for evolution equations utilizing arbitrary finite element discretizations. Adv Comput Math. 2018;44(6):1941-78.

37. Grepl MA, Patera AT. A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations. ESAIM: M2AN. 2005;39(1):157-81.

38. Gubisch M, Volkwein S. Chapter 1: proper orthogonal decomposition for linear-quadratic optimal control: Theory and algorithms. In: Model Reduction and Approximation, pages 3-63. Society for Industrial and Applied Mathematics, 2017.

39. Gunzburger MD, Peterson JS, Shadid JN. Reduced-order modeling of time-dependent PDEs with multiple parameters in the boundary data. Comput Methods Appl Mech Eng. 2007;196(4-6):1030-47.

40. Guzman H A F. Domain decomposition methods in geomechanics. Doctoral Thesis, University of Texas at Austin, 2012.

41. Harris CR, Millman KJ, van der Walt SJ, Gommers R, Virtanen P, Cournapeau D, Wieser E, Taylor J, Berg S, Smith NJ, Kern R, Picus M, Hoyer S, van Kerkwijk MH, Brett M, Haldane A, del Río JF, Wiebe M, Peterson P, Gérard-Marchant P, Sheppard K, Reddy T, Weckesser W, Abbasi H, Gohlke C, Oliphant TE. Array programming with NumPy. Nature. 2020;585(7825):357-62.

42. Hong Q, Kraus J, Lymbery M, Philo F. Conservative discretizations and parameter-robust preconditioners for biot and multiple-network flux-based poroelasticity models. Numer Linear Algebra Appl. 2019;26(4): e2242.

43. Kergrene K, Chamoin L, Laforest M, Prudhomme S. On a goal-oriented version of the proper generalized decomposition method. J Sci Comput. 2019;81(1):92-111.

44. Kim J, Tchelepi H, Juanes R. Stability and convergence of sequential methods for coupled flow and geomechanics: fixed-stress and fixed-strain splits. Comp Methods Appl Mech Energy. 2011;200(13-16):1591-606.

45. Kunisch K, Volkwein S. Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics. SIAM J Numer Anal. 2002;40(2):492-515.

46. Kühl N, Fischer H, Hinze M, Rung T. An incremental singular value decomposition approach for large-scale spatially parallel & distributed but temporally serial data-applied to technical flows. Comput Phys Commun. 2024;296: 109022.

47. Lassila T, Manzoni A, Quarteroni A, Rozza G. Model Order Reduction in Fluid Dynamics: Challenges and Perspectives. In: Reduced Order Methods for Modeling and Computational Reduction, volume 9, pages 235-273. Springer International Publishing, 2014.

48. Lee S, Lee Y-J, Wheeler MF. A locally conservative enriched Galerkin approximation and efficient solver for elliptic and parabolic problems. SIAM J Sci Comput. 2016;38(3):A1404-29.

49. Lewis RW, Schrefler B. The finite element method in the static and dynamic deformation and consolidation of porous media, 2nd Edition. Wiley. 1999.

50. Liu R. Discontinuous Galerkin Finite Element Solution for Poromechanics. PhD thesis, The University of Texas at Austin, 2004.

51. Mandel J. Consolidation des sols (Étude mathématique). Géotechnique. 1953;3(7):287-99.
52. Meyer M, Matthies HG. Efficient model reduction in non-linear dynamics using the Karhunen-Loève expansion and dual-weighted-residual methods. Comput Mech. 2003;31(1):179-91.

53. Mikelić A, Wheeler MF. Convergence of iterative coupling for coupled flow and geomechanics. Comput Geosci. 2012;17(3):455-62.

54. Molenaar J. Multigrid methods for fully implicit oil reservoir simulation. Technical report, TWI: Delft University of Technology; 1995.

55. Murad MA, Thomée V, Loula AFD. Asymptotic behavior of semidiscrete finite-element approximations of biot's consolidation problem. SIAM J Numer Anal. 1996;33(3):1065-83.

56. Peherstorfer B, Willcox K. Online adaptive model reduction for nonlinear systems via low-rank updates. SIAM J Sci Comput. 2015;37(4):A2123-50.

57. Peherstorfer B, Willcox K. Data-driven operator inference for nonintrusive projection-based model reduction. Comput Methods Appl Mech Eng. 2016;306:196-215.

58. Phillips P. Finite element methods for linear poroelasticity: theoretical and computational results. PhD thesis, The University of Texas at Austin, 2005.

59. Phillips PJ, Wheeler MF. A coupling of mixed and discontinuous galerkin finite-element methods for poroelasticity. Comput Geosci. 2008;12(4):417-35.

60. Phillips PJ, Wheeler MF. Overcoming the problem of locking in linear elasticity and poroelasticity: an heuristic approach. Comput Geosci. 2009;13(1):5-12.

61. Ravindran SS. A reduced-order approach for optimal control of fluids using proper orthogonal decomposition. Int J Numer Methods Fluids. 2000;34(5):425-48.

62. Richter T, Wick T. Variational localizations of the dual weighted residual estimator. J Comput Appl Math. 2015;279:192-208.

63. Roth J, Soszyńska M, Richter T, Wick T. A monolithic space-time temporal multirate finite element framework for interface and volume coupled problems. J Comput Appl Math. 2024;446:115831.

64. Roth J, Thiele JP, Köcher U, Wick T. Tensor-product space-time goal-oriented error control and adaptivity with partition-of-unity dual-weighted residuals for nonstationary flow problems. Comput Methods Appl Math. 2024;24(1):185-214.

65. Rozza G, Huynh DBP, Patera AT. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. Arch Computat Methods Eng. 2008;15(3):229.

66. Schmich M. Adaptive finite element methods for computing nonstationary incompressible flows. Doctoral Thesis, Heidelberg University, 2009. ISSN: 0001-0200.

67. Showalter R. Diffusion in poro-elastic media. J Math Anal Appl. 2000;251(1):310-40.

68. Sirovich L. Turbulence and the dynamics of coherent structures. I. Coherent structures. Q Appl Math. 1987;45(3):561-71.

69. Sun S, Liu J. A locally conservative finite element method based on piecewise constant enrichment of the continuous galerkin method. SIAM J Sci Comput. 2009;31(4):2528-48.

70. Thiele J P, Wick T. Numerical modeling and open-source implementation of variational partition-of-unity localizations of space-time dual-weighted residual estimators for parabolic problems. J Sci Comput. 2024;99(25).

71. Tolstoy I. Acoustic, elasticity, and thermodynamics of porous media. New York: Twenty-one papers by M.A. Biot. Acoustical Society of America; 1992.

72. van Duijn CJ, Mikelić A. Mathematical proof of the mandel-cryer effect in poroelasticity. Multiscale Model Simul. 2021;19(1):550-67.

73. van Duijn CJ, Mikelić A, Wick T. Mandel's problem as a benchmark for two-dimensional nonlinear poroelasticity. Appl Anal. 2022;101(12):4267-93.

74. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, Burovski E, Peterson P, Weckesser W, Bright J, van der Walt SJ, Brett M, Wilson J, Millman KJ, Mayorov N, Nelson ARJ, Jones E, Kern R, Larson E, Carey CJ, Polat I, Feng Y, Moore EW, VanderPlas J, Laxalde D, Perktold J, Cimrman R, Henriksen I, Quintero EA, Harris CR, Archibald AM, Ribeiro AH, Pedregosa F, van Mulbregt P. SciPy 1.0 Contributors. SciPy 1.0 Fundamental Algorithms for Scientific Computing in Python. Nat Methods. 2020;17:261-72.

75. Wheeler MF, Wick T, Lee S. IPACS: integrated phase-field advanced crack propagation simulator. An adaptive, parallel, physics-based-discretization phase-field framework for fracture propagation in porous media. Comput Methods Appl Mech Energy. 2020;367: 113124.

76. White JA, Borja RL. Block-preconditioned Newton-Krylov solvers for fully coupled flow and geomechanics. Comput Geosci. 2011;15(4):647-59.

77. White JA, Castelletto N, Tchelepi HA. Block-partitioned solvers for coupled poromechanics: a unified framework. Comput Methods Appl Mech Eng. 2016;303:55-74.

78. Wick T. Multiphysics phase-field fracture: modeling, adaptive discretizations, and solvers. Berlin, Boston: De Gruyter; 2020.

79. Wick T. Space-time methods: formulations, discretization, solution, goal-oriented error control and adaptivity. https://thomaswick.org/links/Wi23_st_book_preprint_Aug_8_2023.pdf, to appear in Springer, Compact Textbooks in Mathematics, 2023.

80. Willcox K, Peraire J. Balanced model reduction via the proper orthogonal decomposition. AIAA J. 2002;40(11):2323-30.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
52. Meyer M, Matthies HG. Efficient model reduction in non-linear dynamics using the Karhunen-Loève expansion and dual-weighted-residual methods. Comput Mech. 2003;31(1):179-91.

53. Mikelić A, Wheeler MF. Convergence of iterative coupling for coupled flow and geomechanics. Comput Geosci. 2012;17(3):455-62.

54. Molenaar J. Multigrid methods for fully implicit oil reservoir simulation. Technical report, TWI: Delft University of Technology; 1995.

55. Murad MA, Thomée V, Loula AFD. Asymptotic behavior of semidiscrete finite-element approximations of biot's consolidation problem. SIAM J Numer Anal. 1996;33(3):1065-83.

56. Peherstorfer B, Willcox K. Online adaptive model reduction for nonlinear systems via low-rank updates. SIAM J Sci Comput. 2015;37(4):A2123-50.

57. Peherstorfer B, Willcox K. Data-driven operator inference for nonintrusive projection-based model reduction. Comput Methods Appl Mech Eng. 2016;306:196-215.

58. Phillips P. Finite element methods for linear poroelasticity: theoretical and computational results. PhD thesis, The University of Texas at Austin, 2005.

59. Phillips PJ, Wheeler MF. A coupling of mixed and discontinuous galerkin finite-element methods for poroelasticity. Comput Geosci. 2008;12(4):417-35.

60. Phillips PJ, Wheeler MF. Overcoming the problem of locking in linear elasticity and poroelasticity: an heuristic approach. Comput Geosci. 2009;13(1):5-12.

61. Ravindran SS. A reduced-order approach for optimal control of fluids using proper orthogonal decomposition. Int J Numer Methods Fluids. 2000;34(5):425-48.

62. Richter T, Wick T. Variational localizations of the dual weighted residual estimator. J Comput Appl Math. 2015;279:192-208.

63. Roth J, Soszyńska M, Richter T, Wick T. A monolithic space-time temporal multirate finite element framework for interface and volume coupled problems. J Comput Appl Math. 2024;446:115831.

64. Roth J, Thiele JP, Köcher U, Wick T. Tensor-product space-time goal-oriented error control and adaptivity with partition-of-unity dual-weighted residuals for nonstationary flow problems. Comput Methods Appl Math. 2024;24(1):185-214.

65. Rozza G, Huynh DBP, Patera AT. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. Arch Computat Methods Eng. 2008;15(3):229.

66. Schmich M. Adaptive finite element methods for computing nonstationary incompressible flows. Doctoral Thesis, Heidelberg University, 2009. ISSN: 0001-0200.

67. Showalter R. Diffusion in poro-elastic media. J Math Anal Appl. 2000;251(1):310-40.

68. Sirovich L. Turbulence and the dynamics of coherent structures. I. Coherent structures. Q Appl Math. 1987;45(3):561-71.

69. Sun S, Liu J. A locally conservative finite element method based on piecewise constant enrichment of the continuous galerkin method. SIAM J Sci Comput. 2009;31(4):2528-48.

70. Thiele J P, Wick T. Numerical modeling and open-source implementation of variational partition-of-unity localizations of space-time dual-weighted residual estimators for parabolic problems. J Sci Comput. 2024;99(25).

71. Tolstoy I. Acoustic, elasticity, and thermodynamics of porous media. New York: Twenty-one papers by M.A. Biot. Acoustical Society of America; 1992.

72. van Duijn CJ, Mikelić A. Mathematical proof of the mandel-cryer effect in poroelasticity. Multiscale Model Simul. 2021;19(1):550-67.

73. van Duijn CJ, Mikelić A, Wick T. Mandel's problem as a benchmark for two-dimensional nonlinear poroelasticity. Appl Anal. 2022;101(12):4267-93.

74. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, Burovski E, Peterson P, Weckesser W, Bright J, van der Walt SJ, Brett M, Wilson J, Millman KJ, Mayorov N, Nelson ARJ, Jones E, Kern R, Larson E, Carey CJ, Polat I, Feng Y, Moore EW, VanderPlas J, Laxalde D, Perktold J, Cimrman R, Henriksen I, Quintero EA, Harris CR, Archibald AM, Ribeiro AH, Pedregosa F, van Mulbregt P. SciPy 1.0 Contributors. SciPy 1.0 Fundamental Algorithms for Scientific Computing in Python. Nat Methods. 2020;17:261-72.

75. Wheeler MF, Wick T, Lee S. IPACS: integrated phase-field advanced crack propagation simulator. An adaptive, parallel, physics-based-discretization phase-field framework for fracture propagation in porous media. Comput Methods Appl Mech Energy. 2020;367: 113124.

76. White JA, Borja RL. Block-preconditioned Newton-Krylov solvers for fully coupled flow and geomechanics. Comput Geosci. 2011;15(4):647-59.

77. White JA, Castelletto N, Tchelepi HA. Block-partitioned solvers for coupled poromechanics: a unified framework. Comput Methods Appl Mech Eng. 2016;303:55-74.

78. Wick T. Multiphysics phase-field fracture: modeling, adaptive discretizations, and solvers. Berlin, Boston: De Gruyter; 2020.

79. Wick T. Space-time methods: formulations, discretization, solution, goal-oriented error control and adaptivity. https://thomaswick.org/links/Wi23_st_book_preprint_Aug_8_2023.pdf, to appear in Springer, Compact Textbooks in Mathematics, 2023.

80. Willcox K, Peraire J. Balanced model reduction via the proper orthogonal decomposition. AIAA J. 2002;40(11):2323-30.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.