Monolithic and splitting based solution schemes for fully coupled quasi-static thermo-poroelasticity with nonlinear convective transport

This paper concerns splitting-based iterative procedures for the coupled nonlinear thermo-poroelasticity model problem. The thermo-poroelastic model problem we consider is formulated as a three-field system of PDE's, consisting of an energy balance equation, a mass balance equation and a momentum balance equation, where the primary variables are temperature, fluid pressure, and elastic displacement. Due to the presence of a nonlinear convective transport term in the energy balance equation, it is convenient to have access to both the pressure and temperature gradients. Hence, we introduce these as two additional variables and extend the original three-field model to a five-field model. For the numerical solution of this five-field formulation, we compare three approaches that differ by how we treat the coupling/decoupling between the flow and/from heat and/from mechanics; these approaches have in common a simultaneous application of the fixed-stress splitting scheme on both the non-linearity and the coupling structure of the problem. More precisely, the derived procedures transform a nonlinear and fully coupled problem into a set of simpler subproblems to be solved sequentially in an iterative fashion. We provide a convergence proof for the derived algorithms, and validate our results through several numerical examples.


Problem statement
The field of poroelasticity aims to describe the interaction between viscous fluid flow and elastic solid deformation within a porous material, and it was pioneered through the works of K. Terzhagi [1] and M. A. Biot [2,3]. In the fully-✩ This work forms part of Norwegian Research Council project 250223. saturated, quasi-static regime, the mathematical modeling of such processes constitutes a coupled two-field linear model where the primary variables are the fluid pressure and the elastic displacement of the solid. This is known as the quasi-static Biot model.
In many important applications, such as geothermal energy extraction, nuclear waste disposal and carbon storage, temperature also plays a vital role and must therefore be included in the aforementioned model. Thus, we consider here a thermo-poroelastic system which can be seen as a generalization of the Biot system to the non-isothermal case; i.e., the coupled processes are heat, flow, and geomechanics. Since it is the cornerstone of many complex models, we focus on the following nonlinear and coupled quasi-static thermo-poroelastic equations as described in [4][5][6]: Find the temperature T , the pressure p, and the displacement u such that ∂ t ψ(p, u, T ) + c f (K∇p) · ∇T − ∇ · (Θ ∇T ) = z, in Ω × (0, t f ), (1.1a) −∇ · θ(u) + α∇p + β∇T = f, in Ω × (0, t f ), (1.1b) ∂ t ϕ(p, T , u) − ∇ · (K∇p) = g, in Ω × (0, t f ), (1.1c) T = 0, u = 0, p = 0, on ∂Ω × (0, t f ), (1.1d) T (·, 0) = T 0 , u(·, 0) = u 0 , p(·, 0) = p 0 , in Ω. (1.1e) In the above model, Ω is a bounded (connected and open) domain in R d , d = 2 or 3, and t f > 0 is the final time. The function z is the heat source, g is the mass source, and f is the body force. The functionals ψ and ϕ denote the heat content and fluid content, respectively; i.e., ψ(p, u, T ) := a 0 T − b 0 p + β∇ · u, and ϕ(p, u, T ) := c 0 p − b 0 T + α∇ · u, where c 0 is the constrained-specific storage coefficient, a 0 is the effective volumetric heat capacity divided by reference temperature, b 0 is the thermal dilation coefficient, α is the Biot-Willis constant, and β is the thermal stress coefficient. The parameter c f is the volumetric heat capacity of the fluid divided by reference temperature, K = (K ij ) d i,j=1 is the permeability divided by fluid viscosity, and Θ = (Θ ij ) d i,j=1 is the effective thermal conductivity divided by reference temperature. The function θ denotes the effective stress tensor, i.e., θ(u) := 2µε(u) + λ∇ · uI, where ε(u) := (∇u + ∇u ⊤ )/2 the symmetric part of ∇u, and I is the identity tensor. Finally, T 0 is the initial temperature, u 0 is the initial displacement and p 0 is the initial pressure. For the present purposes we consider (1.1a)-(1.1e) to be given in dimensionless form, i.e., coefficients and variables are without units.
Note that the above model introduces a nonlinearity in a coupling term, which is the convective transport term in the energy balance equation (1.1a). The presence of this nonlinear coupling term strongly complicates the problem compared to the isothermal case (i.e., to the linear Biot model). Note that if b 0 = β = 0, the flow and mechanics decouples from the heat, and Biot's model is recovered. For the derivation of the constitutive equations of thermo-poroelasticity we refer to the works [6][7][8], and particularly to [4][5][6] where the above model was derived within the framework of the two-scale asymptotic expansion method (see, e.g., [9] for a review of this technique).

Remark 1.1 (Boundary Condition).
We present the problem (1.1a)-(1.1e) with homogeneous Dirichlet boundary conditions only to keep the following presentation as concise as possible. Extending to non-homogeneous or Neumann boundary conditions is straightforward. All results presented in the sequel are valid also for Neumann boundary conditions.

Weak solution and well-posedness of the continuous problem
The common structure of mathematical models that are based on (systems of) scalar conservation laws of the form (1.1a) where nonlinear gradient terms appear suggests introducing the heat flux, r := −Θ∇T , or the Darcy flux, w := −K∇p, as an additional variable. Thus, either the term , e.g., [10,11]. Precisely, it is well known that such terms, dealing non-linearly with the coupled convection, can be quite difficult to approximate correctly in their actual forms. This altogether leads to challenging numerical issues. Furthermore, the choice to introduce the heat flux or the Darcy flux as a new variable depends strongly on which process (flow or heat) that dominates and may result in a different treatment of the convective term. Here, to avoid some of these complexities, we adopt from [12] the mixed form for both the heat and flow subproblems (1.1a) and (1.1b), keeping in mind that Mixed Finite Element (also Finite Volume) literature has developed techniques to handle convective terms [13,14]. Throughout the paper we consider the following assumptions to hold true: The source terms are such that z, g ∈ L 2 (0, t f ; L 2 (Ω)) and f ∈ H 1 (0, t f ; L 2 (Ω)). Furthermore, z, g and f are piecewise constant in time with respect to the temporal mesh of Section 2.
The above variational problem was analyzed in [12]. There, it was shown that under the assumption that the heat flux (or Darcy flux) is such that r(t) ∈ (L ∞ (Ω)) d , for t ∈ (0, t f ), the problem (1.2) has a unique weak solution. Moreover, it was shown that with additional regularity on the data, i.e., , the fluxes are bounded functions. We also note that in [12] constraints on the parameters similar to (A4) were needed for the well-posedness of the above variational problem.

Goal and positioning of the paper
The simulation of thermo-poroelasticity problems is challenging due to the coexistence of different physics which necessitates a coupled set of equations. For these types of problems, there are typically three different approaches employed in modeling fluid flow coupled with reservoir geomechanics. They are known as the fully implicit, the explicit (loosely or weakly) coupling, and the splitting-iterative approaches. The main problem for the applicability of the fully implicit approach, which solves simultaneously the above three-processes (flow, heat and mechanics) problem, is that it results in a very large discrete system of equations to be solved at each time step. Moreover, it does not facilitate the (re-)use of existing codes dedicated to the various subproblems. On the other hand, the fully coupled approach has excellent stability properties [15,16]. An alternative is a weakly coupled approach, which results in a smaller discrete system and a lower computational cost compared to the fully implicit (monolithic) approach. On the other hand, accuracy may be sacrificed, and the sequential approach is only conditionally stable [17,18]. Herein, we adopt an iterative coupling approach, which provides a compromise between the implicit and explicit: At each iteration it has the cost of the sequential approach, yet it converges to the fully coupled implicit approach. We implement the idea of iterative coupling by resolving iteratively the two/three subsystems (depending on the choice of splitting procedure) and by exchanging the values of the shared state variables in an iterative fashion using a general framework of linearly stabilized schemes [19,20].
We argue that adopting an iterative method for the nonlinear and fully coupled three-processes problem, can be considered almost essential for efficient simulation, since the fully coupled approach leads to a prohibitively large system (particularly if MFE methods are adopted [15,[21][22][23]), incorporating different equations that are varied in type and with nonlinearities. The advantage of the iterative approaches considered in this paper is that, at each iteration, smaller, easier-to-solve systems are coupled iteratively through algorithms [22,24]. Another advantage that distinguishes our approaches is the possibility of reusing existing codes for different numerical schemes and coupling techniques specialized to each component of the problem (see e.g., [25,26]). For classical linear poroelasticity, the iterative coupling procedures mentioned above has been studied extensively [19,20,[27][28][29][30][31][32][33]. In particular, two such algorithms have received considerable attention: The "Undrained Split"(constant fluid mass during structure deformation) and the "Fixed Stress Split"(constant volumetric mean total stress during solution of flow problem). In [30], these were first shown to be unconditionally stable. In [20,32] contraction estimates and rates of convergence were derived.
The Undrained Split/Fixed Stress Split algorithms have been generalized in the context of the so-called L-schemes. In the context of coupled problems, these schemes involve adding an artificial stabilization term to one or more of the subproblems with a parameter L > 0. Here, the quantity held constant while solving of one of the subproblems needs not have any physical interpretation. In this sense, the L-scheme generalizes the Undrained Split/Fixed Stress Split algorithms and, due to the removal of physical constraints on the stabilization terms, allows for further optimization. The L-scheme can also be employed as a linearization procedure for nonlinear problems, with the parameter L > 0 mimicking the Jacobian from Newton iteration. To determine the parameter L > 0 for any given problem, derived convergence estimates are necessary. The L-scheme has been shown to perform robustly for Richards equation [34,35], for both linear and nonlinear coupled flow and geomechanics [19,36], for unsaturated/variably saturated porous media [37,38], for twophase flow [39], and for nonlinear diffusion problems [40]. In this paper, we utilize the L-scheme framework both as a decoupling strategy and as a linearization method.
Although the literature on iterative coupling procedures for (isothermal) poroelastic problems is quite extensive, thermo-poroelastic problems have not received the same amount of attention. Sequential iterative methods for linear thermo-poroelasticity were considered in [41]. Iterative splitting schemes for separate poroelasticity and thermoelasticity problems were considered in [42]. Compared to problems of (two-field) coupled flow and mechanics (which can be solved either sequentially or monolithically) we now have additional options in terms of partial decoupling, i.e., solving two of the subproblems together decoupled from the third. Combinatorially, this yields six variations of iterative procedures, ranging from monolithic to fully decoupled. In this work, we focus on the algorithmic developments necessary to handle the nonlinear coupling structure of the problem and propose and analyze all six iterative algorithms for nonlinear thermoporoelasticity. In particular, we employ variations of the L-scheme in all six algorithms, with artificial stabilization terms added to both the flow and heat subproblems. By proving a contraction of all schemes, we obtain explicit expressions for the linearization parameters that guarantee the stability and convergence of all schemes. The main advantage of the L-scheme is that it treats simultaneously the coupling and the non-linearity effects. Thus, no inner iterative approaches are required; see e.g., [43] where L-scheme type approaches are developed to treat iteratively a combined domain decomposition and nonlinearity problem. In most cases, the convergence is linear in the required energy norms. Furthermore, the necessary constraint on the time step is not severe.
The reason we propose six algorithms is the following: The coupling strength of the heat, flow and mechanics may vary depending on the physics at hand. Moreover, the practitioner may have access to existing software of various capabilities. Precisely, to develop robust and efficient solution procedures for the three-processes problem at hand, one should in principle take into account which process (the mechanics and/or flow and/or heat flow) dominate the full problem. In practice, one must also take into account implementation time and available frameworks. Thus, to be agnostic towards the dominating processes and other real-world constraints, we derive a complete framework for this model problem. The six variations of iterative coupling/decoupling algorithms for thermo-poroelasticity cover all possibilities of varying coupling strength between the three physical processes involved. Note that the developed algorithms are applicable on any numerical schemes used to obtain the solutions of the different processes [44,45]. For the convergence analysis, we derive energy-type estimates from which we infer the convergence of the iterate solutions as well as obtaining strict lower bounds on the stabilization parameters, and an upper bound on the time step. A "cut-off"operator M is introduced in the mixed setting in order to make the iterative schemes converge, but we emphasize that this does not affect the model in practice. Several numerical tests validate our proposed algorithms. In particular, we show that by using the derived stabilization estimates, the proposed algorithms perform robustly with respect to both mesh refinement and a wide range of different problem parameters.
The article is organized as follows: In Section 2 we present the fully discrete formulation of the thermo-poroelastic model, and in Section 3 we present all six iterative algorithms. In Section 4, convergence analysis based on contraction estimates is derived, from which the well-posedness of the discrete scheme is inferred in addition to the bounds on the stabilization parameters and time step. In Section 5 we provide several numerical experiments, and finally in Section 6 some concluding remarks.

Discrete setting
Let X h be a simplicial mesh of Ω, matching in the sense that for two distinct elements of X h their intersection is either an empty set or their common vertex or edge. Let h K denote the diameter of K ∈ X h and let h be the largest diameter of all such triangles, i.e., h := max K ∈X h h K . For the time partition, we let {t n : n = 0, 1, . . . , N} be the discrete time steps, where 0 := t 0 < t 1 < · · · < t N = t f , and let τ n = t n − t n−1 , n ≥ 1, be the difference between consecutive discrete times.
In other words, we have t n := ∑ n ℓ=1 τ ℓ , 1 ≤ n ≤ N, and therefrom t f = ∑ N n=1 τ n . For the discrete spaces, we let T h , R h , P h , W h and U h be suitable finite element spaces corresponding to the infinite dimensional spaces of Section 1.2, where we assume that For the time discretization we employ a backward Euler scheme. For the sake of simplicity, we take the source terms f, g and z to be piecewise constant in time. We denote by (T n h , r n h , p n h , w n h , u n h ) the discrete counterpart of the solution tuple to problem (1.2) at time t n .
Before giving the discrete version of the variational formulation (1.2a)-(1.2e), we need to introduce the so-called cut-off operator M as described in e.g., [10,11] as where M is a (large) positive constant. We note that the introduction of this operator in the following discrete variational formulation has little or no practical implications, but is necessary in order to facilitate the convergence analysis.
Obviously, if the exact fluxes are bounded, i.e., w n , r n ∈ (L ∞ (Ω)) d , and if we pick M large enough, we have M(w n )(x) = w n (x) and M(r n )(x) = r n (x). Thus, a precise value for the constant M is not necessary.
In the above scheme, we used (M(w n h ) · Θ −1 M(r n h ), S h ) for the approximation of the convective coupling term instead of the original (w n h ·Θ −1 r n h , S h ). The reason for this approximation will be clarified later. Eqs. (2.3a)-(2.3b) form the discrete mixed scheme of the heat subproblem, (2.3c)-(2.3d) form the discrete mixed scheme for the flow subproblem, and (2.3e) is the discrete form of the mechanics subproblem with the Galerkin finite element method. Together, these subproblems make up the nonlinear and fully coupled discrete version of the thermo-poroelastic problem. In the next section, their iterative solution procedure is detailed.
where two different cut-off operators, M and R are used (defined with different constants M and R, respectively). In that case, the underlying iterative methods of Section 3 as well as the convergence analysis of Section 4 remains true with minor modifications in the proofs. For simplicity, we let M = R (and thus M = R).

Remark 2.2 (Existence of M).
It was shown in [46] for a related poroelastic model that if the flux is bounded on the continuous level, then the discretized flux will inherit this property. Thus, with sufficient regularity of the domain, source and initial data, the existence of the constant M is guaranteed.

The L-type iterative schemes
We now present six iterative (splitting) algorithms for the discrete thermo-poroelastic problem (2.3). These algorithms involve either decoupling all the subproblems and solving each separately at every iteration (three-step algorithm), or decoupling only one subproblem from the other two which are then solved together (two-step algorithm), or solving a linearized problem monolithically at every iteration (one-step algorithm). We use the letters H (Heat), F (Flow), and M (Mechanics), to abbreviate the algorithms, e.g., a two-step algorithm where the heat and flow subproblems are solved together decoupled from the mechanics subproblem is referred to as (HF-M) and similarly for other combinations of coupling/decoupling of the subproblems. Throughout the rest of the article we will mostly refer to the discrete problems and therefore omit the h-subscript on the variables and test functions for cleaner notation. We shall also denote the time step simply by τ , keeping in mind it may depend on n.
At the time step n ≥ 1, let (T n−1 , r n−1 , p n−1 , w n−1 , u n−1 ) be given. We then approximate the solution at the actual time step n ∈ {1, . . . , N} using the sequence (T n,k , r n,k , p n,k , w n,k , u n,k ) for k ≥ 0, defined in an iterative fashion, and where the iterate (T n,0 , r n,0 , p n,0 , w n,0 , u n,0 ) is an initial guess (e.g., the solutions at the previous time step). All the algorithms involve adding the stabilization terms L T (T n,k − T n,k−1 , S) and L p (p n,k − p n,k−1 , q) to the left hand sides of Eqs. (2.3a) and (2.3c), respectively, where L T , L p > 0 are the stabilization parameters (to be chosen later). Furthermore, to make the notation easier, we introduce the parametrized fluid and heat content functionals: For a given L T , L p > 0, we define We are now able to present our six iterative algorithms:

The monolithic scheme (HFM)
At the each iteration k ≥ 1 of the L-type monolithic scheme, we solve the linearized thermo-poroelastic problem: Given (T n,k−1 , p n,k−1 , w n,k−1 ), find (T n,k , r n,k , p n,k , w n,k , u n,k ) such that This algorithm is continued until a fixed tolerance is reached. Clearly, in the above algorithm, the L-scheme acts only as a linearization procedure, where we approximate the convective transport term by M(w n,k−1 ) · Θ −1 M(r n,k ). One can also approximate this term by M(w n,k ) · Θ −1 M(r n,k−1 ), the analysis presented next remains true and follows exactly the same lines. The complexity in this algorithm is that it requires solving a large system generated by (3.2), which combines equations varied in type, and this at each iteration k ≥ 1. This encourages the development of efficient techniques for the resolution of these coupled systems.

The partially decoupled schemes
In the second set of iterative schemes, we only decouple the flow (F), mechanics (M) or heat (H) from the remaining two processes, which are being solved monolithically. Thus, we transform the monolithic solver (HFM) into a two-level iterative approach in which two simpler subproblems are solved sequentially. For the partially and fully decoupled schemes, we do not consider cyclical permutations of the order in which the subproblems are solved to yield different algorithms. The partially decoupled setting delivers the following three iterative approaches:

(HF-M): Coupled heat and flow
Decoupling the mechanics calculation from the coupled flow and heat flow calculation, the first two-level iterative scheme reads as follows: At the iteration k ≥ 1, do:

(HM-F): Coupled heat and mechanics
The second scheme in this subsection is obtained by decoupling the flow calculation from the remaining coupled thermo-elasticity calculation. This iterative scheme reads: At the iteration k ≥ 1, do:

(FM-H): Coupled flow and mechanics
The last two-level scheme is obtained by decoupling the poro-elasticity calculation (solved monolithically) from the heat flow. Note that a similar scheme was proposed in [47] for two-phase flow. This iterative scheme reads: At the iteration (3.5c) • Step 2: Given (p n,k , w n,k , u n,k , T n,k−1 ), find (T n,k , r n,k ) such that (a 0 + L T )(T n,k , S)

The fully decoupled schemes
In this set of iterative coupling schemes, we simply split the three processes, providing three subproblems to be solved sequentially. Fixing the mechanics calculation in the third level, two approaches are then derived in which either the problem of flow or heat transfer is solved first followed by solving the other system and then the mechanics using the already calculated information. These schemes lead to solving much simpler subsystems through the algorithm. In addition, they enable the reuse of existing codes for each component of the problem.

Convergence analysis
The starting point for our analysis is the existence and uniqueness of a solution to (2.3). To this aim, we will make use of the following Lemma (cf. [11]), stating the Lipschitz property of the cut-off operator M:  (4.2b) The proof of the next Theorem is based on showing that the scheme (3.2) is a contraction and then, by applying the Banach fixed-point theorem [48], deduce convergence of the scheme. In what follows, we will frequently use the following polarization and binomial identities: Hence, we can state the first of our main results: Then, the monolithic L-scheme HFM (Algorithm 3.1) defines a contraction satisfying (4.6) Therefrom, the limit is the unique solution of the problem (2.3).

Remark 4.1 (Bound on Time
Step). Note that a 0 − b 0 > 0 due to (A4) and that by the choice of M sufficiently large (thus, the right hand side of (4.5) is a positive number). If a priori bounds on the fluxes are available, and these are small enough such that M can be chosen to yield equality in (4.7), then there is no constraint on the time step.
We choose now S = e k T , y = τ e k r , q = e k p , z = τ e k w , and v = e k u as test functions in Eqs. (4.8a)-(4.8e), respectively. Then, summing the resulting equations and using the identity (4.3), together with applying the Cauchy-Schwarz and Young inequalities and some algebraic manipulations, we get, which leads to  (e k T , e k r , e k p , e k w , e k u ) := (T n,k − T n , r n,k − r n , p n,k − p n , w n,k − w n , u n,k − u n ). (4.13) The second of our main results is given through:

Theorem 4.3 (Convergence of the Partially Decoupled Schemes). Assuming (A1)-(A6) holds true, the stabilization parameters are such that
, (4.14) and the time step satisfies (4.5), then the partially decoupled L-scheme HF-M (Algorithm 3.2.1) is a contraction given by (4.16) Proof. We start by taking the difference of Eqs. (3.3a)-(3.3e) at iteration k with the corresponding equations solved by (T n , r n , p n , w n , u n ). This leads to the following set of difference equations: The aim now is to show a contraction of successive error functions, thereby implying convergence of the sequences (T n,k , r n,k , p n,k , w n,k , u n,k ) as k → ∞ for n ≥ 1. Taking as test functions q = e k p , z = τ e k w , S = e k T , y = τ e k r , and v = e k−1 u in (4.17a)-(4.17e), respectively, and adding the resulting equations together leads to Let now ξ ∈ (0, 1) and rewrite the above estimate as ) . (4.20) We now follow [19] and choose ξ = 2µ 2µ + dλ , which together with the Young inequality yields (4.21) Combining this with Eq. (4.18) leads to We thus need to impose some constraints on the stabilization parameters, i.e., L p ≥ .
With this, we can discard some positive terms on the left hand side of (4.22), and use the Cauchy-Schwarz and Young inequalities, together with the Lipschitz property of M to obtain Thus, if the time step τ satisfies (4.5), we can write (4.25) as for some ξ ∈ (0, 1). Following the same steps which led to (4.21), and choosing as before ξ = 2µ 2µ + dλ , we get by the Young inequality This shows a contraction of the residuals and therefore completes the proof. □ Before we state the last of our main results, we let the difference functions defined in (4.13) now be the difference between the solutions at the iteration k of problem (3.7) and the solutions to (2.3). The last of our main results then reads:

Numerical experiments
In the following, we present three numerical test cases using the algorithms from Section 3. The first is a constructed problem, posed on the unit square domain, with prescribed solutions for the temperature, pressure and displacements. Here, we consider five different parameter regimes, exhausting all possibilities of weak/strong coupling between the subproblems, and compare the number of iterations needed for convergence with decreasing mesh sizes for both stabilized and non-stabilized algorithms. Since analytical solutions are available, we present also discretization errors.
Next, we present two implementations of Mandel's problem [50], which is originally a benchmark problem in linear poroelasticity, extended here to nonlinear thermo-poroelasticity. For the original Mandel problem, analytical solutions for the pressure and displacement field are known. Due to the similarity of the thermo-poroelastic equations with the linear Biot's equations, and due to the lack of benchmark problems for thermo-poroelasticity, we choose to use this problem for our second and third numerical test cases. Even though the analytical solutions are no longer valid when including temperature, we have sufficiently weak temperature effects in the first implementation of Mandel's problem that the computed pressure and displacement field matches the (isothermal) analytical solutions. The second implementation of Mandel's problem includes a heat source, which has a significant effect on both the pressure and displacement. Regarding the spatial discretization, we choose the following finite element spaces: T h , P h := {ϕ ∈ L 2 (Ω) : ∀K ∈ X h , ϕ| K ∈ P 0 (K )}, where RT 0 (K ) denotes the lowest-order Raviart-Thomas finite-dimensional subspace associated with the element K ∈ X h , and P l (K ) is the space of polynomials on K ∈ X h of total degree less than or equal to l. Thus, the spaces (T h , R h ) and (P h , W h ) are the lowest order Raviart-Thomas mixed finite element spaces for the mixed flow and heat flow subproblems,  respectively. Note that both spaces satisfy the condition (2.1), see e.g., [51] for more details on (mixed) finite elements. The vector valued space U h is the first order Lagrange finite element space for the mechanics problem. We employ the following stopping criterion for the iterative algorithms, given in terms of the relative and absolute tolerances, aTOL and rTOL, i.e., where we set aTOL = rTOL = 1e − 6 for all the computations. For the solution of the linear subproblems, we make use of a direct sparse linear solver from the Python library SciPy [52], i.e., scipy.sparse.linalg.spsolve. The present approaches can also be combined with iterative solvers adapted to the various subproblems. All numerical tests are implemented in a finite element code written in Python, the complete source code is accessible at https://github.com/ matkbrun/FEM. Remark 5.1 (Stability). The discretization defined by (5.1a)-(5.1c) does not satisfy inf-sup stability unless U h is chosen as the space of piecewise quadratic polynomials [53]. We employ here a different stabilization strategy, which is iterative coupling using artificial stabilization parameters. Such strategies (e.g., the Fixed Stress Splitting algorithm) have proven very successful for stabilizing problems of poroelasticity [28].

Remark 5.2 (Parameter Robustness).
We believe that if the differences c 0 − b 0 and/or a 0 − b 0 are very small, the problem may become difficult to treat. In such cases it can be advantageous to formulate the mechanics problem like the Stokes equations, as done in [54]. However, such an investigation is outside the scope of the present paper. For more details on parameter robustness for the related problem of poroelasticity, see [53,55].

Test case 1: Example with manufactured solution
As a first test case, we let the domain be a regular triangularization of the unit square, i.e., Ω = [0, 1] × [0, 1] ⊂ R 2 , and prescribe the following smooth solutions for the temperature, pressure and displacement: For the analysis and comparison of our algorithms, we consider dimensionless equations, i.e., all parameters are set to 1.0e − 1, except for the three coupling coefficients {α, β, b 0 }, which we vary in order to weaken/strengthen the coupling between the three subproblems. In particular, we consider five different parameter regimes, PR1-PR5, specified in Table 1: We also set a 0 = c 0 = 2b 0 , thus satisfying (A4). We emphasize that the parameter regimes PR1-PR5 are not intended to have any physical meaning, they are only constructed in order to test the convergence properties of the proposed algorithms. Table 2 shows number of iterations needed for convergence using the six algorithms from Sections 3.1, 3.2 and 3.3 for a single time step with decreasing mesh sizes and stabilization parameters chosen according to equality in (4.14). We see that for the parameter regimes 1, 3 and 4 we have higher iterations numbers than for parameter regimes 2 and 5, for all six algorithms. This is because L T ∼ β 2 and L p ∼ α 2 , and larger stabilization results in higher iteration numbers.
Furthermore, as expected, the strongly coupled parameter regime (PR1) yields the highest iteration numbers, in particular for the algorithms HF-M, H-F-M and F-H-M. Apart from this, the algorithms are performing robustly both with respect to different coupling regimes and decreasing mesh sizes. For comparison we also provide in Table 3 the results without stabilization, i.e., L T = L p = 0.  We see here that the fully monolithic algorithm (HFM) has low iteration counts for all parameter regimes since this is only a linearization scheme, which does not require stabilization (cf. Theorem 4.2). For the two-level (Section 3.2) and three-level (Section 3.3) algorithms, which involves some splitting as well as linearization, we see that iteration counts for different parameter regimes correspond to the various coupling/decoupling of the subproblems present in the algorithms: Splitting of subproblems that are strongly coupled yields high iteration numbers compared to solving the strongly coupled subproblems together. This is in contrast to employing stabilization, which greatly improves the robustness of the algorithms with respect to variations in parameters. For the strongly coupled parameter regime (PR1), we even have no convergence for algorithm HF-M when no stabilization is applied.
Furthermore, in order to check the robustness of the proposed schemes with respect to the nonlinear coupling, we adjust the coefficient of the convective term, c f , in order to make this dominate the flow/heat coupling. Table 4 shows the number of iterations needed for convergence when c f = 1.0e1 for both the strongly coupled parameter regime (PR1) and the weakly coupled parameter regime (PR5). We also compare the results when no stabilization is applied. Note that we here only use a single mesh with h = 1/16.
For the weakly coupled parameter regime (PR5), there is no difference in iteration numbers between the stabilized and non-stabilized algorithms, even with a dominating nonlinearity. For the strongly coupled parameter regime (PR1), the stabilized algorithms have a significantly lower iteration count. This might be due to the fact that the nonlinearity appears as a coupling term.  Since analytical solutions are available for this problem, we provide also the discretization errors, denoted by (e h,T , e h,r , e h,p , e h,w , e h,u ), measured in the L 2 -norm. Due to almost no variation in discretization errors between the six algorithms and between the different parameter regimes (less than 5%), we provide in Table 5 the discretization errors using algorithm H-F-M applied on the weakly coupled parameter regime (PR5). We also include the convergence rates, defined by r T := e h j ,T /e h j+1 ,T , and similarly for the other variables.

Test case 2: Mandel's problem
We refer to [56] for a detailed description of Mandel's problem. Formulas for the analytical pressure and displacements can be found in [44]. We provide here only a brief description. Mandel's problem is posed on a rectangular domain representing a poroelastic slab of extent 2a in the horizontal direction, 2b in the vertical direction, and infinitely long in the third direction. The poroelastic slab is contained between two rigid plates, where at the initial time a downward force of magnitude 2F is applied to the top plate, with an equal but opposite force applied to the bottom plate. The top and bottom boundaries are treated as impermeable, while zero pressure (and temperature) is prescribed at the right and left boundaries. Due to the nature of Mandel's problem, the pressure, temperature and horizontal component of the displacement varies only in the horizontal direction, while the vertical component of the displacement varies only in the vertical direction. From symmetry considerations, it suffices to consider only the top right quarter rectangle, i.e., the We perform now all computations with dimensional equations and realistic choices of physical parameters. In particular, we take mechanics and flow parameters identical to [20], and heat parameters identical to [41]. However, in [41] the flow-heat coupling coefficient b 0 is taken to be identically zero, hence in order to preserve this coupling we instead choose a suitably small number (which satisfies (A4)). All parameters are listed in Table 6.
In terms of our previous notation, we now have 4) and the constraints (A4) should be understood as The magnitude of the compressive force is F = 2 × 10 8 Pa m, and the physical dimensions of the quarter rectangle is given by a = 100 m and b = 10 m, of which we make a regular triangularization. We impose the compressive force as a Dirichlet boundary condition on the top boundary (x 2 = b) for the vertical component of the displacement. We denote by n 1 and n 2 the number of subdivisions of the domain in the x 1 and x 2 directions, respectively. For the first implementation of Mandel's problem we prescribe homogeneous boundary conditions and zero source term and initial condition for the heat problem. Fig. 1 shows the solution profiles for the pressure, temperature and displacements for selected time steps, with the analytical (isothermal) solutions for the pressure and displacement included for comparison.
The computed solutions for pressure and displacement matches the analytical solutions, even though the analytical solutions are only valid for the linear isothermal problem. This is because the induced temperature effect in the system is   Mandel's problem: Solution profiles at t ∈ {100 s, 500 s, 1000 s}, computed using the monolithic scheme HFM, with z = 2 × 10 −4 W m −3 K −1 , and n 1 = n 2 = 40.

Table 7
Mandel's problem: Number of iterations with decreasing mesh sizes for Mandel's problem. Stabilization from theory.
Heat source  10  18  18  14  13  14  14  20  18  18  13  13  13  12  40  18  18  13 13 13 12 small enough that the heat decouples from the flow and mechanics. For the second implementation of Mandel's problem we prescribe a constant source term for the heat problem, i.e., z = 2 × 10 −4 W m −3 K −1 and zero initial condition. Fig. 2 shows the solution profiles for the pressure, temperature and displacements at selected time steps. The temperature source now interacts with the other processes and thus has an effect on the pressure and horizontal component of the displacement. Furthermore, the temperature change in the system is now increasing with increasing time. Table 7 shows the number of iterations for Mandel's problem using the derived algorithms.

Conclusions
Based on previous developments of iterative splitting schemes from linear poroelasticity, we have proposed six novel iterative procedures for nonlinear thermo-poroelasticity. These algorithms are using stabilization and linearization techniques similar to [19,34], which is known in the literature as the 'L-scheme'. The thermo-poroelastic problem we consider can be viewed as a coupling of three physical processes (or subproblems): Flow, geomechanics and heat. Solving this system either monolithically (all three subproblems simultaneously), partially decoupled (two subproblems simultaneously), or fully decoupled (each subproblem separately), yields six possible combinations of coupling/decoupling, which we have used to design our six algorithms. All of these involve a linearization of the convective term and added stabilization terms to both the flow and heat subproblems. In this sense, our use of the L-scheme is both as a stabilization for iterative splitting and a linearization of nonlinear problems.
For any given situation, the coupling strength between the three subproblems may vary. A-priori, the expectation is that solving together subproblems that are strongly coupled yields better efficiency and accuracy than splitting strategies. On the other hand, if the coupling between two or more subproblems is weak, a splitting procedure might be beneficial. For this reason, and due to the fact that splitting the three-way coupled multi-physics problem into smaller subproblems allows for combining existing codes that separately can handle any of the three processes involved (or two of them combined), six different algorithms are presented. These six algorithms cover all possibilities of strong/weak coupling between the three subproblems. Using the well-posedness of the continuous problem, we obtained lower bounds on the stabilization parameters and proved the convergence of our proposed algorithms under a constraint on the time step. In practice, however, we find that this bound is not tight; as long as the fluxes are not becoming unbounded (e.g., due to a singularity), a 'reasonable' time step can safely be chosen.
Our algorithms are tested in detail with several numerical examples. In particular, we find that all six algorithms are performing robustly with respect to both mesh refinement and different parameter regimes (i.e., strong/weak coupling between the subproblems and strong/weak nonlinear effects), using the stabilization revealed by our analysis. We also find that using no stabilization results in the algorithms being more sensitive to the parameter regimes, i.e., splitting subproblems that are strongly coupled yields high iteration numbers compared to solving these subproblems together. This phenomenon is also observed in the stabilized algorithms, but to a significantly lesser extent. In particular, our conclusion is that with no stabilization, each of the algorithms is suitable only for a certain parameter regime (i.e., one that corresponds to the coupling/decoupling structure present in the algorithm). This is in contrast to the stabilized algorithms, which can handle a much wider range of different parameter regimes.