Projection-based reduced order modeling of an iterative scheme for linear thermo-poroelasticity

This paper explores an iterative coupling approach to solve linear thermo-poroelasticity problems, with its application as a high-fidelity discretization utilizing finite elements during the training of projection-based reduced order models. One of the main challenges in addressing coupled multi-physics problems is the complexity and computational expenses involved. In this study, we introduce a decoupled iterative solution approach, integrated with reduced order modeling, aimed at augmenting the efficiency of the computational algorithm. The iterative coupling technique we employ builds upon the established fixed-stress splitting scheme that has been extensively investigated for Biot's poroelasticity. By leveraging solutions derived from this coupled iterative scheme, the reduced order model employs an additional Galerkin projection onto a reduced basis space formed by a small number of modes obtained through proper orthogonal decomposition. The effectiveness of the proposed algorithm is demonstrated through numerical experiments, showcasing its computational prowess.


Introduction
Thermo-hydro-mechanical (THM) processes refer to the coupled interactions between temperature, fluid flow, and mechanical deformation that occur in a wide range of natural and engineered systems.These processes are relevant in various fields, including environmental science, civil engineering, and material science, to name a few.Therefore, the ability to understand and predict tightly coupled THM processes in natural and engineering systems has significant impacts, for example, on the environment, public safety, and the economy.
One widely employed mathematical model for describing THM processes is based on Biot's non-isothermal consolidation theory [1], known as the thermo-poroelasticity model.This model extends Biot's poroelasticity model [2], which characterizes the interaction between a deformable porous medium and fluid flow within it under isothermal conditions.The governing system of partial differential equations (PDEs) of the thermo-poroelasticity model consists of a heat transfer equation, mass conservation equation, and momentum conservation equation.These equations are fully coupled through linear and nonlinear coupling terms.For simplicity, this paper focuses on a fully-coupled, linearized thermo-poroelasticity model, where the heat transfer is assumed to be diffusion-dominated and the coupling strength between the THM processes is rather weak.
Due to the complex nature of the mathematical model for coupled multi-physics thermo-poroelasticity, it is highly challenging to develop an accurate and computationally efficient numerical method.There are various approaches to handling the coupled nature of the problems.First, one can solve the coupled system using a monolithic numerical method, where all subproblems (PDE equations) are solved simultaneously in each time step.In general, monolithic schemes are unconditionally stable but computationally expensive and may require some special linear or nonlinear solvers, often making them less practical.Another approach is to use a sequential method, where the subproblems are separated and solved sequentially in each time step.This approach is much more economical compared to the monolithic method.However, a sequential method may not converge to the right solution, or effective stabilization is often required to make such methods unconditionally stable [3].To overcome this limitation, one can resort to an iterative approach, which is a staggered yet tightly coupled scheme.In iterative coupling schemes, the subproblems are solved in a staggered way while ensuring full convergence of the solution in each time step.
Iterative coupling schemes for poroelasticity have been extensively studied in recent years [4,5,6,7,8,9,10,11].In particular, one of the most successful approaches is based on the fixed-stress spitting method [12,13,14,15,16].In this fixed-stress scheme, one solves the flow problem while fixing the volumetric stress (mean stress), subsequently solving the mechanics problem using the updated pressure solution.The convergence and stability of this fixed-stress splitting method have been studied in the context of various numerical methods, including CG, DG, and mixed finite element methods.In contrast to poroelasticity, the subject of iterative methods for thermo-poroelasticity has not been heavily studied.Some recent papers on this topic include [17,3,18].
This paper first describes an iterative scheme for solving the thermal-poroelasticity model based on the classical continuous Galerkin finite element method, which is considered to be the high-fidelity solver or offline solver in the context of the reduced order modeling approach.We especially utilize the fixed stress iterative method and derive the stability and convergence properties that are close to those of the monolithic approach.
Next, we utilize this fixed stress iterative scheme to train a projection-based reduced order model (ROM) as outlined in [19,20].ROMs generally consist of two distinct phases: the offline (or training) stage, and the online (or evaluation) phase.In the offline stage, a high-fidelity scheme is employed to generate a dataset containing solutions for the specific problem.In our context, we utilize the fixed stress iterative method, and this dataset is compressed using proper orthogonal decomposition [21], yielding a reduced basis space of significantly lower dimension.During the online phase, the ROM incorporates a Galerkin method that leverages the reduced basis space, in contrast to the original high-fidelity space.If the reduced basis space possesses dimensions substantially smaller than its high-fidelity counterpart, and if the problem operators can be suitably pre-assembled, then the ROM evaluation is anticipated to be significantly more computationally efficient than assessing the high-fidelity scheme.The application of projection-based ROMs to THM systems is still limited, see in particular [22,23].Furthermore, the use of iterative schemes as high-fidelity solvers in ROMs [24,25,26,27,28,29,30,31] is not as common as monolithic ones, and none of the aforementioned references uses a scheme with the specific properties of the fixed-stress iterative scheme.Therefore, by applying a fixed-stress iterative ROM to a THM problem, this work advances the state-ofthe-art and could pave the way for applications of the proposed ROM in geothermal flows.
The rest of the paper is organized as follows.We describe the governing equations for the thermo-poroelasticity model in Section 2. We then introduce the high-fidelity discretization based on the fixed-stress iterative scheme in Section 3. In there, we also prove that the solution of the iterative scheme converges to the monolithic scheme, assuming some conditions for certain physical parameters.In Section 4, we introduce the projection-based ROM, which builds upon the fixed-stress iterative scheme, and show that its solution converges to the one of the ROM trained on the monolithic scheme.Finally, we present several numerical experiments in Section 5 to validate and present the capabilities of the proposed algorithm.

Governing Equations
Let Ω be a bounded, connected, and Lipschitz domain in R d , d = 2, 3, with the boundary ∂Ω and let I = (0, T ] with T > 0.Then, let u : Ω → R d be the vector-valued displacement of the solid, p : Ω → R the scalar-valued fluid pressure, and θ : Ω → R the scalar-valued temperature.Then, governing equations for thermo-poroelasticity are derived by coupling momentum balance for mechanics based on linear elasticity, mass balance for the pressure, and energy balance for the temperature as follows: In the momentum balance equation (1a), σ(u) is the standard stress tensor from linear elasticity.It satisfies the constitutive equation σ(u) := 2µϵ(u) + λ(∇ • u)I, where ϵ(u) := 1 2 [∇u + (∇u) T ] is the strain tensor, I is the d × d identity tensor, and µ, λ are the Lamé constants.The Lamé constants are assumed to be in the range µ ∈ [µ 0 , µ 1 ] and λ ∈ [0, ∞) for some 0 < µ 0 < µ 1 < ∞.Also, f is the body force, α is the Biot-Willis constant, K dr := (dλ + 2µ)/d is the drained isothermal bulk modulus, and α T is the volumetric skeleton thermal dilation coefficient.The total stress tensor is given by σ(u, p, θ) = σ(u) − αpI − 3α T K dr θI.
The second equation (1b) is the mass balance equation for the fluid, assuming the Darcy law for the volumetric fluid flux: q = −K∇p.We ignore the gravity effect and set the fluid viscosity to be one here for a simple presentation of the numerical method.However, including the gravity term and the fluid viscosity in the numerical formulation is straightforward.Here, K ∈ R d×d is the permeability tensor, which is symmetric and uniformly positive-definite and satisfies the following assumption: there exist positive constants k min , k max such that for any x ∈ Ω, In addition, c 0 = 1/M , where M is Biot's modulus, 3α m is the thermal dilation coefficient, and g is the volumetric fluid source/sink term.Finally, the energy balance equation, or the heat transfer equation (1c), is obtained by assuming local thermal equilibrium between solid and fluid in pores.Therefore, this energy balance equation is expressed in terms of a single temperature variable θ with the effective total heat conductivity C d .Additionally, η is the volumetric heat source/sink term, and θ 0 is a reference temperature, which is assumed to be nonzero.In addition, the bulk thermal conductivity tensor D is symmetric and uniformly positive-definite and assumed to satisfy the following: there exist positive constants d min , d max such that for any x ∈ Ω, Note that we adopt a linearized heat transfer equation in this paper, contrary to using a nonlinear equation in other work (e.g.[3]).This linearization of the model by introducing the reference temperature and dropping the nonlinear convection term can be justified if we assume small magnitudes of α T and α m and that diffusion is a dominant process for heat transfer.Similar linearized models were used in [32,33].To complete the system of governing equations (1), we have to provide initial conditions and boundary conditions.As our main focus here is on the iterative coupling scheme, we will consider the pure homogeneous Dirichlet boundary conditions for all three variables, whereas the initial conditions are given as

High-fidelity Discretization
In this section, we present the discretization schemes employed during the reduced order model training in section 4, which will be referred to as high-fidelity schemes therein.We discuss first the monolithic scheme (section 3.1), then introduce the proposed fixed-stress iterative scheme (section 3.2).In section 3.3, we prove a convergence result for the proposed fixed-stress iterative scheme by utilizing continuous Galerkin finite element methods.

Monolithic high-fidelity (M-HF) method
In this section, we present a weak formulation and a fully discrete continuous Galerkin method for the model problem (1).The standard notation for the L 2 -and Sobolev spaces and their associated inner products and norms will be used here.In particular, (•, •) and ∥ • ∥ denote the L 2 inner product and L 2 -norm, respectively, and (•, •) 1 and ∥ • ∥ 1 denote the (full) H 1 -inner product and (full) H 1 -norm, respectively.
To derive a weak formulation, we multiply (1a), (1b), and (1c , and s ∈ H 1 0 (Ω), respectively, and integrate by parts, resulting in the following weak problem: At every t ∈ (0, T ], find (u, p, θ) for any (v, w, s) To present a fully-discrete method, first let ∆t > 0 be a time step size and t n = n∆t, where n = 0, . . ., N with t N = T .For temporal discretization, we consider the backward Euler time-stepping scheme for simplicity.For spatial discretization, let us consider a shape-regular triangulation T h of Ω.Then, we let V h , W h , Θ h be finite-dimensional subspaces of [H 1 0 (Ω)] d , H 1 0 (Ω), H 1 0 (Ω), respectively, on the mesh T h .Then, the fully-discrete monolithic CG method for (5) reads as follows: Given for any (v, w, s) The well-posedness and convergence of the CG method for the same model as (1) but with an extra nonlinear coupling term −K∇p n • ∇θ n+1 in the temperature equation were studied in [34].In their paper, they employed the standard continuous piecewise P k finite element for all V h , W h , Θ h for k ≥ 1, where P k is a space of all polynomials of degree at most k.This space, V h × W h × Θ h = P k × P k × P k , does not satisfy the inf-sup stability condition.Despite the lack of the inf-sup stability condition, however, they proved that the resulting method is still well-posed and convergent with some constraints on the physical parameters.These constraints are the same as the ones for the convergence of our FS-HF iterative scheme, to be established in (16).

Fixed-stress high-fidelity (FS-HF) iterative scheme
Solving the monolithic system resulting from the M-HF method ( 6) is very expensive computationally.Our goal is to develop an iteratively coupled method whose solutions converge to the solution of the M-HF method.To this end, we will employ a version of the fixed-stress splitting scheme.Similarly to [3,18], the main idea behind the fixed-stress split is to fix (or lag) the total mean stress while solving the flow and heat problems.Here, the total mean stress, σ, is defined as σ = tr( σ)/d, where tr(•) is the trace operator on tensors.By noting the relation between volumetric stress and strain, we can see that In light of this relationship, let σ n+1 h be the approximate total mean stress at time t n+1 , defined by Now, to present our fixed-stress iterative scheme, we will denote by ϕ n+1,i the i-th iterate of ϕ n+1 for any variable ϕ = u h , p h , θ h .To derive a fixed-stress iterative scheme, we want the condition while solving the flow and heat equations.Note that the flow and heat equations can be solved simultaneously or sequentially, leading to different versions of the fixed-stress scheme.In this section, we consider a method where we solve these two problems sequentially, solving the flow equation first, then the heat problem, before solving the mechanics problem.Since we will solve the flow and heat equations sequentially, we will use a modified version of (8a) when solving them.Specifically, we will assume that the temperature and pressure variables, besides the total mean stress, are fixed while solving the flow and heat problems, respectively.This means that we use when solving the flow problem and when solving the heat problem.
Given an initialization u n+1,0 h , p n+1,0 h and θ n+1,0 h , the following is our FS-HF iterative scheme at each time step n+1 to generate infinite sequences {u n+1,i+1 h Note that L > 0 to appear in the algorithm is a stabilization coefficient.If we drive a scheme using (8b) and (8a), then L = 1.However, one can consider tuning this coefficient to achieve faster convergence, see for instance [35].
Step 1.Given (u n+1,i h Steps 1, 2, and 3 are repeated by increasing i to i+1 until an appropriate stopping criterion is satisfied.The specific choice of the stopping criterion depending on a given problem, as well as the procedure to initialize u n+1,0 h , p n+1,0 h , and θ n+1,0 h , will be discussed in Section 5.

Convergence analysis of the FS-HF scheme
In this section, we will prove that the solution of the FS-HF scheme, (9), converges to that of the M-HF scheme, (6).To facilitate the analysis, we introduce some notations for the errors in the iterates.Let e n+1,i h , respectively.To derive a system of error equations, subtract (6a), (6b), and (6c) from (9c), (9a), and (9b), respectively, and do some algebraic manipulations.We note that the resulting system involves quantities only at the current time step n + 1.Therefore, in the following, we will use superscripts indicating only the iteration numbers.Then, the system of error equations reads as follows: for any (v, w, s) For the subsequent analysis, we provide two preliminary results here.
Lemma 3.1.The following coercivity condition is satisfied for any v ∈ V h .
Proof.It is easy to check that ∥∇ • v∥ 2 ≤ d∥ε(v)∥ 2 using Young's inequality.Then, (11) follows from the definition of K dr .□ Lemma 3.2.The following inequality holds: Proof.First, subtract the mechanics equation (10a) at the iteration i from the iteration i + 1: Then, take v = e i+1 u −e i u in (13) and use the coercivity condition (11) on the left-hand side and the Cauchy-Schwarz inequality on the right-hand side to obtain which leads to Using (15) back in (14), we now get from which we can get the desired result (12) after multiplying out the terms and using Young's inequality. and For a fixed n + 1, assume that the FS-HF solution at time step n + 1 is initialized to the M-HF solution at time step n, i.e.
Then, the FS-HF solution defined in (9) converges to the solution of the M-HF method (6): Proof.We first recall that holds for any e i , e i+1 in L 2 (Ω).Now, take v = e i+1 u , w = ∆t e i+1 p , s = ∆t e i+1 θ /θ 0 in (10), add the three equations, and apply (18) to obtain The last two terms on the left-hand side of the above equation can be replaced by u − e i u ) in (10a).Then, applying Young's inequality to the last two terms on the left-hand side and to the two terms on the right-hand side, we obtain where δ > 0 is a positive constant.After collecting terms with the superscript i + 1 on the left-hand side and the terms with i on the right-hand side, we bound the lefthand side below.Specifcally, we use the assumptions in ( 16) for the coefficients of ∥e i+1 p ∥ 2 and ∥e i+1 θ ∥ 2 , and the uniform ellipticity ( 2) and (3), and Poincaré inequality for the flow and heat flux terms, and (12) for the terms involving e i+1 u − e i u .Then, we can arrive at the following inequality: where C Ω is a constant from the Poincaré inequality.Now, due to the assumptions in ( 17), both 1 − 1 2δ and L 2 − δ are nonnegative, hence we have Then, we have then, (20) can be rewritten as Therefore, (e i p , e i θ ) → 0 as i → ∞, which also implies that ∥e i p ∥ → 0 and ∥e i θ ∥ → 0. In order to prove the convergence of ∥e in the same manner as in the proof of Lemma 3.2, and recall Korn's inequality, stating that there is a constant C korn > 0 such that Using these two results, we have As the quantity on the right-hand side goes to 0 as i → ∞, ∥e i u ∥ 1 also goes to 0 as i → ∞. □ Remark 3.1.Similar constraints to the ones in (16) were needed to prove the wellposedness and convergence of numerical methods for thermo-poroelasticity in [36,17,37].On the other hand, Kim [3] proved unconditional stability for an extended fixed-stress split for a nonlinear thermo-poroelasticity model.

Projection-Based Reduced Order Models
In this section, we discuss projection-based reduced order models (ROMs) based on compression of the time history of the solutions provided by one of the highfidelity methods (either M-HF or FS-HF method) from Section 3.Such compression comprises the offline stage (or training) of the ROM and will result in a reduced basis of small dimension r.The reduced basis will be obtained by means of the proper orthogonal decomposition (POD) algorithm, see Section 4.1.The online stage (or evaluation) of the ROM will then consist of a Galerkin method applied on either the fully-coupled scheme (Section 4.2) or the fixed-stress iterative scheme (Section 4.3).Convergence of the fixed-stress iterative ROM scheme to the solution of the fullycoupled ROM scheme will be proven in section 4.4.

Proper orthogonal decomposition
Let ϕ = u h , p h , θ h denote one of the components of the solution field, the displacement, pressure, and temperature, respectively.Then, {ϕ n } N n=0 denote the sequence obtained by collecting such a component at different times, t 0 , . . ., t N .The proper orthogonal decomposition [21,20,19] is a data analysis method that aims to obtain an equivalent representation of the sequence {ϕ n } N n=0 as a linear combination of at most N + 1 orthogonal modes φ 0 , . . ., φ N and define information content indices ν 0 , . . ., ν N associated to each corresponding mode.Assuming the latter to be sorted as ν 0 ≥ ν 1 ≥ . . .≥ ν N , POD is typically used as a data compression tool to obtain an approximate representation in the linear space Φ r = span(φ 0 , . . ., φ r−1 ), truncating the combination to the first r modes for some r ≤ N and neglecting the modes φ r , . . ., φ N associated to indices ν r , . . ., ν N with small information content.
In the context of ROMs, the primary interest of the application of POD is in the generation of the space Φ r , to be called the reduced basis space.To this end, we will focus on the so-called method of snapshots, even though alternative presentations are possible, for instance by means of the singular value decomposition or the principal component analysis [21,20,19].First of all, the method of snapshots requires solving the following eigenvalue problem Since the matrix ) is symmetric positive definite by construction, its eigenvalues ν 0 , . . ., ν N are real and positive and, without loss of generality, can be assumed to be sorted in decreasing order.Denoting by v n ∈ R N +1 the eigenvector associated with the eigenvalue ν n , the n-th POD mode is then computed as Indeed, the obtained POD modes are (•, •) 1 -orthonormal.To see this, first note that v T n v m = δ nm , where δ nm is the Kronecker delta function.Then, we obtain Finally, we notice that POD modes φ 0 , . . ., φ N satisfy homogeneous Dirichlet boundary conditions if the input sequence ϕ 0 , . . ., ϕ N is zero on ∂Ω.This easily follows by taking the trace of (24) on ∂Ω.Therefore, any element in Φ r is by construction equal to zero on ∂Ω.

Monolithic reduced order model (M-ROM)
In this section, we describe a projection-based ROM, which builds upon the fully-coupled monolithic (M-HF) scheme introduced in Section 3.1 as a high-fidelity discretization.During the offline stage, the high-fidelity discretization is queried to obtain three sequences, {u n h } N n=0 , {p n h } N n=0 , and {θ n h } N n=0 , representing the time history of the displacement, pressure, and temperature fields, respectively.Upon selecting a reduced basis size r ≤ N and applying POD to each sequence {u n h } N n=0 , {p n h } N n=0 , and {θ n h } N n=0 as discussed in Section 4.1, we obtain the following reduced basis spaces respectively.
During the online stage, the ROM is a Galerkin method for problem (5) on the reduced basis space V r × W r × Θ r .Therefore, the fully-coupled monolithic ROM (M-ROM) scheme for (5) reads as follows: for any (v, w, s) ∈ V r × W r × Θ r .Even though ( 6) and ( 26) are both obtained by means of a Galerkin method, the fundamental difference is that the reduced basis space V r × W r × Θ r employed in ( 26) is of small dimension 3r, owing to the compression carried out by POD.
In order to highlight further differences between ( 6) and ( 26), let u n+1 r ∈ R r be the vector whose components are the M-ROM degrees of freedom for u n+1 r ∈ V r , i.e.
Similarly, let p n+1 r ∈ R r and θ n+1 r ∈ R r be the vector whose components are the M-ROM degrees of freedom for p n+1 r ∈ W r and θ n+1 r ∈ Θ r , respectively.The system (26) can thus be written in the following block-matrix form where the expressions of matrices and vectors appearing in the block formulation are summarized in Table 1.Indeed, matrices in Table 1 can be assembled once and for all at the end of the offline stage since such computations involve integration over the finite element mesh T h .Similarly, vectors in Table 1 can be pre-assembled at the end of the offline stage for every n = 1, . . ., N and stored.During the online stage, iterating in time through (27) can finally be carried out at a vastly decreased computational cost since i) the assembly of the block system only requires loading r×r matrices and vectors of dimension r in Table 1 from storage without necessitating any operation involving the finite element T h , and ii) the resulting linear system is of small size 3r × 3r.

Fixed-stress reduced order model (FS-ROM)
In this section, we further describe a projection-based ROM built upon the fixedstress high-fidelity (FS-HF) iterative scheme described in Section 3.2.During the offline stage, the FS-HF scheme is queried to obtain the time evolution {u n,∞ h } N n=0 , {p n,∞ h } N n=0 and {θ n,∞ h } N n=0 of displacement, pressure and temperature fields, respectively, where the superscript n, ∞ denotes the converged solutions at time step n.We notice that in practice the number of iterations will not be ∞ upon defining suitable stopping criteria in Section 5.1.Those sequences are then compressed by means of POD to obtain the reduced basis spaces V F S r , W F S r and Θ F S r upon proceeding as in Section 4.2.We notice that the obtained reduced basis spaces V F S r , W F S r and Θ F S r are in principle different from V r , W r and Θ r obtained in Section 4.2 since the former is obtained by applying POD to solutions of the FS-HF scheme, while the latter is the result of a compression of solutions computed by the M-HF scheme.For the sake of a simpler notation, however, we will drop the suffix F S from the reduced basis spaces V F S r , W F S r , and Θ F S r , still understanding that reduced basis spaces employed in this section are generated from the FS-HF scheme, rather than the M-HF one.
Steps 1-ROM, 2-ROM, 3-ROM are repeated by increasing i to i + 1 until appropriate stopping criteria are satisfied.The specific choice of the stopping criteria, as well as the procedure to initialize u n+1,0 r , p n+1,0 r and θ n+1,0 r , will be discussed in Section 5.1.The FS-ROM scheme can be equivalently reformulated in matrix form as follows.It generates the sequences {u n+1,i+1 r Navier-Stokes problems do not need enrichment.Based on this experience, and the fact that (i) in [39], ROM solutions derived from spaces which lack inf-sup stability are very visibly incorrect, with relative errors up to 100%, (ii) as in [42], our goal is to discuss a scheme with iterative coupling, and (iii) as in [41], the iterative scheme contains some stabilization, in the next numerical results we will not perform any enrichment.This will only empirically confirm that the resulting ROM is inf-sup stable, while further theoretical justifications are of relevant mathematical interest, but beyond the scope of this work.

Convergence analysis of the FS-ROM
In this section, we retrace the convergence analysis carried out in Section 3.3 in order to adapt it to the FS-ROM introduced in Section 4.3.
Theorem 4.1.Assume that ( 16) and (17) hold.Furthermore, assume that the M-ROM online stage employs a Galerkin projection onto the reduced basis spaces associated to FS-ROM.For a fixed n + 1, assume that the FS-ROM solution at time step n + 1 is initialized to the M-ROM solution at time step n, i.e.
Then, the FS-ROM solution defined in (28) converges to the solution of the M-ROM method (26): Upon subtracting (26a), (26b), and (26c) from (28c), (28a), and (28b) we obtain again (10), with the only exception that it must be now intended for any (v, w, s) ∈ V r × W r × Θ r .Also, Lemma 3.1 implies that the coercivity condition (11) also holds for any v ∈ V r , since V r ⊂ V h .Even with the modified definitions of e i u , e i p and e i θ defined above, the proof of Lemma 3.2 and the proof of the convergence in Theorem 3.1 follow in the same manner as the original proofs.□ Remark 4.2.We notice that we cannot provide bounds to with the same technique used in the proofs of Lemma 3.2 and Theorem 3.1.Indeed, if we define our errors as then we can obtain the same error equation, (10), as before but for any (v, w, s) ∈ V r × W r × Θ r .However, since (e i u , e i p , e i θ ) ∈ V h × W h × Θ r , we cannot use these errors as test functions.As a result, the proofs of Lemma 3.2 and Theorem 3.1 would not follow through.Instead, we will empirically quantify (30) by means of the numerical experiments in the next section.

Numerical Experiments
In this section, we provide numerical experiments to validate and demonstrate the computational efficiency of the proposed algorithm.All the following computations are conducted by utilizing FEniCSx1 for HF computations, and RBniCSx2 for ROM computations.FEniCSx and RBniCSx rely on PETSc3 for sparse and dense linear algebra, respectively.RBniCSx further queries LAPACK4 for the solution of the dense eigenproblem (23).

Example 1. Validation of ROM on test cases with smooth analytical solution
In this example, we validate our numerical methods against four different test cases with the given analytical solution, namely where (x, y) ∈ Ω = (0, 1) 2 and t ∈ I = (0, T ] for some T > 0 to be specified in each test case.We note that the factor xy(1 − x)(1 − y) in each solution variable guarantees that the solution satisfies homogeneous Dirichlet boundary conditions, in agreement with the assumptions in the previous sections.The right-hand side functions f , g, and η in (1) are obtained from the expressions in (31).Furthermore, the physical parameters are chosen to be , and D = 10 −5 .
To compute the initial conditions at n = 0, the M-HF and FS-HF schemes employ the linear interpolation of (31), while M-ROM and FS-ROM employ the L 2projection of (31) on the reduced basis spaces.The stabilization parameter L = 1 is employed in the FS-HF and FS-ROM schemes.Initialization of the iterations at time n + 1 is based on the assignment of the converged solution at time n.Iterations for FS-HF will continue until the following conditions are satisfied where ϵ = 10 −10 .Similarly, the FS-ROM iterations will continue until where ∥ • ∥ R r denotes the Euclidean norm; we note that such conditions still represent a stopping criteria on the relative ∥ • ∥ 1 -norm of the increment because, due to (25), the matrix representing the (•, •) 1 inner product on V r , W r and Θ r is the identity.Furthermore, we set the maximum iteration number as 20 throughout this example.

Example 1A. Validation of HF and ROM solvers by a mesh refinement analysis
As the first test case in Example 1, we validate the four distinct solvers discussed in the previous sections: M-HF (Section 3.1), FS-HF (Section 3.2), M-ROM (Section 4.2), and FS-ROM (Section 4.3).In the domain Ω of the unit square and for a final time T = 1, we employ an initial spatial discretization of h = 0.25, linear finite elements, and a temporal discretization of ∆t = 0.0025.Subsequently, with each cycle, we reduce h by half and divide ∆t by four.The convergence rate is then computed by employing corresponding norms to calculate the errors (maximum over (0, T ]) against the provided analytical solutions in (31).
To validate the ROM algorithm, for each cycle, we initially run the HF solver to generate sequences Then, we employ the POD algorithm as discussed in Section 4.1 to compress these sequences, focusing on retaining the first five modes i.e., r ≤ 5.The ROM's performance is evaluated across five distinct scenarios, specifically for values of r equal to 1, 2, 3, 4, and 5.
Results of the convergence analyses are plotted in Figure 1 for the monolithic schemes (M-HF and M-ROM) and in Figure 2 for the fixed-stress iterative schemes  (FS-HF and FS-ROM).Specifically, the left column illustrates results in the ∥ • ∥ 1norm, whereas the panels on the right show results in the ∥•∥-norm.The convergence analysis for the displacement u, pressure p, and temperature θ are arranged in three rows within the figure.The first row pertains to the displacement u, the second row to the pressure p, and the third row to the temperature θ.For the M-HF scheme we observe a convergence rate of 2 for errors in the ∥ • ∥-norm and of 1 for errors in the ∥ • ∥ 1 -norm.Furthermore, since the solution of FS-HF converges to that of M-HF, as presented in Theorem 3.1, FS-HF is expected to converge at the same rate as M-HF.Indeed, numerical results show that FS-HF achieves the same convergence rates.Furthermore, with the exception of r = 1, the ROM scheme also achieves identical convergence rates to the HF scheme used in its training.Figure 3 highlights one of the main advantages of the ROM schemes, which is the reduction of computational cost.This advantage becomes evident through the significant reduction in the number of degrees of freedom achieved by the ROM schemes in comparison to the HF schemes.By accompanying this reduction with the precomputation of matrices and vectors in Table 1 we ensure that the leading cost in the evaluation of the ROM scales only with the ROM dimension r, and is not affected by the HF mesh size.Consequently, ROMs exhibit a remarkable speedup of two orders of magnitude when compared to the underlying HF discretization for both monolithic and fixed-stress cases.
We further remark that, while decreasing the mesh size h demands a higher CPU time for solving with the HF method, the results in Figure 3 show that the ROM  CPU time is not affected by mesh refinements.Indeed, the evaluation of the ROM only requires solutions of r × r linear systems for the FS-ROM or 3r × 3r linear systems for the M-ROM regardless of the value of h.Moreover, increasing r from 1 to 5 does not seem to increase the required CPU time either, since the solution of linear systems of dimension at most 15 × 15 is computationally inexpensive.
On the other hand, a comparison between Figure 3a and Figure 3b highlights that the M-HF solver places lower computational demands than the FS-HF solver.Indeed, the ratio between the CPU time of the M-HF solver and the FS-HF one ranges from 0.13 (largest value of h) to 0.32 (smallest value of h).We anticipate that the M-HF solver's computational cost will eventually surpass that of the FS-HF solver in an asymptotic sense, but we refrain from further reducing the value of h due to the computational costs.A similar comparison illustrates that the M-ROM solver is an order of magnitude less expensive than the FS-ROM solver.Indeed, the ratio between the M-ROM CPU time and the FS-ROM one is around 0.15 for every value of h and r.This is again due to the fact that, in the current realization, solving linear systems of dimension up to 15 × 15 requires a CPU time that is almost independent of the dimension of the linear system.The M-ROM solver needs to solve only one small linear system, while the FS-ROM requires solving 3 linear systems per iteration.
Finally, Figure 4 demonstrates that the average number of iterations necessary by FS-ROM for any given value of r falls within the range of 5 to 6.Moreover, except for the case when r = 1, the FS-ROM exhibits an equivalent number of iterations as the FS-HF solver.

Example 1B. Comparison of HF and ROM solvers in time
In the second test case in Example 1, we focus on the fixed stress iteration scheme and compare the errors between HF and ROM solvers.Here, we use a mesh size of h = 1/16, linear finite elements, and a time step size of ∆t = 0.001.We set the final time to be T = 1, which is larger than the one used in Example 1A.We then train both M-ROM and FS-ROM schemes, compute at most 10 POD modes, and evaluate the ROMs on the same time interval (0, T ] with the same time step size ∆t.The larger number of POD modes selected in this test case compared to Example 1A is due to the fact that the final time is now ten times larger than the one employed in Example 1A.
First, Figure 5 illustrates the relative errors with respect to the analytical solutions in (31) for the fixed-stress schemes.In particular, the left column presents the results in the ∥ • ∥ 1 -norm, whereas the right column shows the results in the ∥ • ∥norm.The first row pertains to the displacement u, the second row to the pressure p, and the third row to the temperature θ.We observe that the ROM for r = 1 is inaccurate, with relative errors up to 100% over the time interval (0, T ], regardless of the solution component u, p or θ.For r ≥ 3, the ROM solutions are as accurate as the HF results, with only minor discrepancies appearing at the beginning of the time interval (0, T ] in Figures 5d and 5f.Furthermore, the error rises as time progresses.
Next, Figure 6 demonstrates the performance of the FS-ROM scheme by comparing its solutions to those of the FS-HF scheme.Figures 6a-6b show that u r converges very fast to u h in the ∥ • ∥ 1 and ∥ • ∥-norm as r increases.In particular, we observe an error reduction of five orders of magnitude when increasing the value of r from r = 1 to r = 5.Similarly, Figures 6c-6d show a decrease of more than one order of magnitude for p within the first five POD modes, and Figures 6e-6f a decrease of almost two orders of magnitude for θ.
Furthermore, Figure 7 reports the comparisons of the total CPU time for both HF and ROM solvers.Analogous to Example 1A, the ROM approach ensures a speedup of over two orders of magnitude, regardless of the specific value of r. Figure 7 also depicts the CPU time necessary for precomputation of the ROM operators summarized in Table 1.This observation leads us to conclude that the precomputation of the ROM operators at the end of the offline stage is essential to achieve the desired speedup.This becomes particularly crucial for larger values of r as attempting to carry out ROM assembly during an online ROM solve would result in an online evaluation that is nearly as computationally intensive as the HF solve itself.
Lastly, we present the fixed-stress iteration counts for the FS-ROM in Figure 8a.Similar to Example 1A, the FS-ROM exhibits an identical number of iterations when compared to the FS-HF solver.In addition, Figure 8b           numbers for each matrix in Step i-ROM, Step ii-ROM, and Step iii-ROM.We observe that the condition numbers increase mildly while r increases.Furthermore, Figure 8c presents the normalized eigenvalues of the POD compression, where ν 0 is the largest eigenvalue; the fast decay of the eigenvalues justifies limiting the ROM dimension to r = 10 at most.
5.1.3.Example 1C.ROM on a larger time interval than the HF solver Employing ROM in Examples 1A and 1B serves well for validation purposes; however, its practical utility remains limited.This is because the ROM, when trained on the time interval (0, 1] and evaluated on the same time interval, merely replicates pre-existing solutions within the same time span.In the context of the third test case, we delve into the ROM's extrapolation potential.Here, we train the model on the time span (0, 0.1] and assess its performance on a more extensive interval, (0, 1].
First, Figure 9 illustrates the relative errors with respect to the analytical solutions in (31) for the FS-HF and FS-ROM schemes in the time interval (0, 1].As in Figure 6, the left column presents the results in the ∥ • ∥ 1 -norm, whereas the right column shows the results in the ∥ • ∥-norm.The first row pertains to the displacement u, the second row to the pressure p, and the third row to the temperature θ.We note that the errors for the FS-HF solver are now limited to the interval (0, 0.1] but the errors displayed in the rest of the time interval [0.1, 1] is only reported for comparison purposes from the previous Example 1B.
Comparing the results in Figure 9 to those obtained in Figure 5 for Example 1B, one can notice that a ROM with r = 5 is not as accurate in Example 1B due to the errors from extrapolation.Still, increasing the reduced basis size to r = 7 allows us to obtain a solution that is as accurate as the solution of the FS-HF scheme.
Although this test case demonstrates a successful extrapolation in accuracy, such an achievement introduces additional challenges for the ROM. Figure 10a indicates that the average number of fixed-stress iterations from the FS-ROM is greater than that of the FS-HF solver for r ≥ 7.In addition, for r ≥ 9, the FS-ROM scheme surpasses the maximum allowable limit of 20 iterations at each time step.It is important to emphasize that this behavior does not contravene Theorem 4.1; rather, it stems from limitations in arithmetic precision.To support this statement, we plot the condition numbers of the matrices appearing on the left-hand side of Step i-ROM, Step ii-ROM, and Step iii-ROM (from the algorithm (29)) in Figure 10b.We notice that the condition numbers of the matrices from Step iii-ROM is above 10 9 for r ≥ 7 and the condition numbers of all matrices from Steps i-ROM, ii-ROM, and iii-ROM are above 10 9 for r ≥ 9. From Figure 10c we understand that this behavior is due to having added to the reduced basis some trailing POD modes associated with          the evaluation of the ROM still relies on time stepping through the interval (0, T ], adopting a larger time step presents a way to further enhance the speed of online evaluation as it entails a reduction in the number of time instances.Here, we set ∆t = 0.001 for FS-HF as in Example 1B, but FS-ROM with ∆t = 0.01, i.e. ten times larger than the one in Example 1B.We consider T = 1. Figure 11 displays the relative errors of the displacement, pressure, and temperature (in ∥ • ∥ 1 and ∥ • ∥ norms) with respect to the analytical solutions in (31) for the FS-HF and FS-ROM schemes.Comparing Figures 5 and 11, we conclude that employing the FS-ROM scheme with a larger time step introduces a larger time discretization error in the approximation of the pressure p and the temperature θ.However, in Figure 11, the errors in p r and θ r still remain in the same order of magnitude as those for p h and θ h , thus offering an approximation that could be accurate enough in several practical scenarios.Furthermore, the displacement u, which exhibits minimal variation over time in this setup, can be precisely approximated even by a ROM that employs a larger time step for time integration.In fact, there are negligible distinctions observed between u r and u h .
As expected, choosing a larger time step decreases the online computational cost.Comparing Figure 12 to Figure 7, we realize that the CPU time is decreased by a factor of ten since the number of time instances is divided by ten.Finally, choosing a larger time step does not have a detrimental effect on the average number of FS-ROM iterations, which are still equal to the ones of the FS-HF scheme, see Figure 13.

Example 2. A domain with parametric heterogeneities
The purpose of this example is to assess the capabilities of the proposed FS-ROM in a parametric setting with a more realistic application problem.The domain Ω = (0, 1) 2 is divided into two subdomains, Ω 1 and Ω 2 , with Ω 1 being characterized by higher permeability and conductivity than Ω 2 .
• Case i) The first testing set is composed of the same 25 parameter values used while training the ROM.These cases are represented by star-shaped markers ( ) in Figure 16.
• Case ii) The second testing set is obtained by taking a finer 7 × 7 uniform grid (i.e., 49 parameters) on P train and discarding any parameter values already present in the training set (e.g., a parameter ω = (−3, −1)), resulting in 40 parameter values.This set is used to evaluate the accuracy of the ROM on parameters that were not seen during training but still belong to the training range.The results are depicted in Figure 16 by diamond-shaped markers (♦).
• Case iii) Finally, the third testing set amounts to sampling a 7 × 7 uniform grid on P test and discarding parameter values lying in P train , for a total of 40 parameter configurations in P test \ P train .These values are employed to assess the accuracy of the ROM while extrapolating outside of the training range and are reported with plus-shaped markers (✚) in Figure 16.
The evaluation of M-ROM and FS-ROM follows through as described in Section 4, with a minor modification regarding the assembly of the ROM operators in Table 1.Here, we consider only [A θθ r ] βγ = (D∇φ θ γ , ∇φ θ β ) for the sake of exposition.Since now D = D(ω), the resulting A θθ r (ω) cannot be preassembled as reported in Table 1 without knowing the value of ω a priori.To retain the computational efficiency, the ROM makes use of the so-called affine parametric dependence or parameter separability [19,20]  are now independent of ω and can be preassembled at the end of the offline stage.
Figure 16 is divided into five panels, representing five different quantities of interest.The first three panels, Figures 16a-16c, depict the ∥ • ∥ 1 -norm relative errors for u r , p r , and θ r with respect to u h , p h and θ h , respectively.In each panel, the top row displays the results obtained from the M-ROM, while the bottom row illustrates those from the FS-ROM.Each row contains four images, corresponding to increasing values of r = 10, r = 30, r = 60 and r = 90.Markers are positioned at the provided ω = (ω 1 , ω 2 ) coordinates, signifying the errors for the three distinct cases: case i) denoted by , case ii) by ♦, and case iii) by ✚.
In the first column of Figures 16a, 16b and 16c, it is evident that both the M-ROM and HF-ROM with r = 10 exhibit inaccuracies.The relative errors are approximately 10 −2 within the training range (inside the inner rectangular region) and can escalate to as much as 1 outside of this range.The accuracy progressively increases when increasing r to 30, 60 and 90.In particular, when r = 90 we observe that both the M-ROM and HF-ROM become more accurate within the training range, including both seen (i.e., in the training set) and unseen data (i.e., not in the training set).We note that, at r = 90, the temperature (θ) exhibits the lowest error within the training range, reaching a maximum of 10 −6 for both the M-ROM and FS-ROM.This outcome aligns with the observation that the POD eigenvalues associated with temperature were the ones that exhibited the most rapid decrease.On the other hand, the displacement u and pressure p show differences of up to an order of magnitude between the M-ROM and FS-ROM for parameter values within the training range.Nevertheless, the resulting relative errors remain at most 10 −4 , a value that falls below the tolerance ϵ established by the FS-HF solver.
We also used both M-ROM and FS-ROM for extrapolation outside of the training range, which resulted in larger relative errors, especially in the bottom left of P test .The highest relative error is observed at ω = (−4, −2) and is of the order of 10 −2 for FS-ROM and 10 −1 for M-ROM for the pressure approximation, i.e. the variable that was characterized by the slowest POD decay.
Finally, Figures 17a and 17b discuss the number of iterations and the computational efficiency of the ROM.Specifically, Figure 17a illustrates the ratio of the total number of FS-ROM iterations to FS-HF iterations.Apart from the smallest values, r = 10 and r = 30, in the extrapolation region, FS-ROM converges in a comparable number of iterations to FS-HF across the entire parameter range.Also, Figure 17b depicts the speedup offered by the ROM, representing the ratio of CPU time required for a HF solve compared to a ROM solve.The magnitude of the speedup for the M-ROM diminishes from 10 4 (for r = 10) to 10 3 (for r = 90) as the reduced basis size r increases.On the other hand, FS-ROM exhibits a smaller speedup, approximately on the order of 10 2 , which remains relatively consistent across varying r values.In summary, any ROM yields query times that are at least 100 times faster than the corresponding HF scheme.

Conclusions
This paper introduces a novel approach involving a fixed-stress iterative for the solution of linear thermo-poroelasticity problems in conjunction with reduced or-der modeling techniques.The approach is validated by means of several numerical examples, which also illustrate its computational capabilities.The benefits of this methodology are two-fold.Firstly, the utilization of fixed-stress iterations aids in the management of complex multi-physics coupling scenarios.Secondly, the incorporation of reduced order modeling significantly boosts computational efficiency.In future works we plan to demonstrate the versatility of the proposed techniques by employing them in conjunction with various numerical discretization techniques, including mixed finite elements, as well as discontinuous or enriched Galerkin finite element methods.

Figure 2 :
Figure 2: Example 1A: convergence of the errors for the fixed-stress iterative schemes (FS-HF and FS-ROM).

Figure 3 :
Figure 3: Example 1A: average CPU time per time step.

Figure 5 :
Figure 5: Example 1B: relative errors (δ = h, r) with respect to the analytical solution for the fixed-stress iterative schemes (FS-HF and FS-ROM).

Figure 6 :
Figure 6: Example 1B: errors between the FS-ROM and the FS-HF scheme.
Average number of FS-ROM iterations.

Figure
Figure Example 1C: relative errors for δ = h, r with respect to the analytical solution for the fixed-stress iterative schemes (FS-HF and FS-ROM).

Figure 11 :
Figure 11: Example 1D: relative errors (δ = h, r) with respect to the analytical solution for the fixed-stress iterative schemes (FS-HF and FS-ROM).