Well-posedness of the fully coupled quasi-static thermo-poro-elastic equations with nonlinear convective transport

This paper is concerned with the analysis of the quasi-static thermo-poroelastic model. This model is nonlinear and includes thermal effects compared to the classical quasi-static poroelastic model (also known as Biot's model). It consists of a momentum balance equation, a mass balance equation, and an energy balance equation, fully coupled and nonlinear due to a convective transport term in the energy balance equation. The aim of this article is to investigate, in the context of mixed formulations, the existence and uniqueness of a weak solution to this model problem. The primary variables in these formulations are the fluid pressure, temperature and elastic displacement as well as the Darcy flux, heat flux and total stress. The well-posedness of a linearized formulation is addressed first through the use of a Galerkin method and suitable a priori estimates. This is used next to study the well-posedness of an iterative solution procedure for the full nonlinear problem. A convergence proof for this algorithm is then inferred by a contraction of successive difference functions of the iterates using suitable norms.


Introduction
The field of poroelasticity is concerned with describing the interaction between viscous fluid flow and elastic solid deformation within a porous material, and goes back to the works of K. Terzhagi [32] and M.A. Biot [6,7]. Porous materials are by definition solid materials comprising a great number of interconnected pores, typically at the order of micrometers, where the interconnectivity of the pores is sufficient to allow for fluid flow through the material. For this reason, porous materials are usually modeled at the ✩ This work forms part of Norwegian Research Council project 250223.
continuum scale, such that the complex micro-structure needs not be explicitly accounted for in the modeling, but rather implicitly through so-called effective parameters such as e.g. porosity and permeability. Porous materials are primarily associated with objects such as rocks and clays, but biological tissue, foams and paper products also fall within this category. Consequently, the field of poroelasticity is of great importance in a range of different engineering disciplines, such as petroleum engineering, agricultural science and biomedicine, among others. A number of comprehensive text books related to the field exists; see e.g. [13,14,36].
Mathematical modeling of fluid saturated deformable porous media on the continuum scale relies on the theory of linear elasticity, adapted to porous materials by using the so-called total stress tensor instead of the Cauchy stress in the momentum balance equation. In particular, the total stress tensor is a linear combination of the Cauchy stress for the empty elastic skeleton and the isotropic stress coming from the fluid, i.e. the pore pressure. Within the quasi-static framework inertial terms are ignored, thus giving a purely elliptic equation for the momentum balance. A second equation of parabolic type accounts for the mass balance as fluid is displaced by the deformation of the solid, and relates change in porosity to volumetric fluid flow, i.e. the Darcy flux. This is essentially Biot's poroelastic model for quasi-static deformation (see e.g. [6,13]). There is an extensive literature on this model problem and on its numerical approximation.
To mention a few, the well-posedness based on the canonical two-field formulation with displacement and pressure as variables was carried out in [30], while three and four-field formulations have also been analyzed (taking Darcy flux and total stress as independent variables), and can be found in several studies, e.g. [1,26,37]. A key feature of this model, one which greatly facilitates the analysis, is the symmetric coupling between the equations.
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 modeling. Using the method of formal two-scale expansions (see e.g. [12,18] for a detailed review of this method), a thermo-poroelastic model was derived in [10], which accounts for fluid pressure, elastic displacement, and temperature distribution within a fine-grained, fully saturated poroelastic material within the framework of quasi-static deformation. This model is similar to other thermo-poroelastic models which exists in the literature; see e.g. [13,16,22,31,34], although there are also some notable differences among these works, in particular from the modeling point of view; i.e. allowable flow rates and deformation, choice of coordinate frames etc. (see [10,34] for a comparison of existing thermo-poroelastic models). However, from the point of view of analysis the important factor is the coupling structure between the equations, and the model we analyze exhibits a fully coupled structure.
The aim of the present work is to establish the well-posedness of the nonlinear thermo-poroelastic model as described in [10], where we also provide a priori energy estimates and regularity properties of the solutions. We restrict our attention to an isotropic material such that the elastic coefficients are given by the Lamé parameters, and the Biot coefficient and thermal stress coefficient are given by scalar quantities. Some algebraic constraints on these coefficients must be imposed in order to obtain our results. Although the literature on the analysis of poroelastic models is quite extensive, there is not much literature on the analysis of thermo-poroelastic models; in [34] a corresponding energy functional for the thermo-poroelastic model was derived. This functional was then shown to be monotonically decreasing in time for a small enough characteristic temperature difference.
We undertake our analysis with a future mixed finite-element implementation in mind, and therefore double the number of variables from three to six, and investigate the existence and uniqueness of a weak solution corresponding to this fully coupled six-field model. The primary variables in this model are; fluid pressure, temperature, elastic displacement, Darcy flux, heat flux, and total stress. This makes the problem suitable for combinations of well-known stable finite-elements, such as Raviart-Thomas(-Nédélec) [25,29] and Arnold-Winther [2,3]. From an implementation point of view there are several advantages of a mixed formulation over the canonical three-field formulation; the discretization respects mass and energy conser-vation, produces continuous normal fluxes regardless of mesh quality, and in general a mixed formulation is advantageous for domain decomposition techniques. We restrict our attention to two spatial dimensions, as this will be the most relevant case for the subsequent work, although the results we present can be extended to higher dimensions in a straightforward manner. In particular, the definition of the isotropic compliance tensor must reflect the choice of spatial dimension.
The main difficulty we face in the following analysis is the nonlinear coupling between the equations, i.e. the nonlinear convective transport term in the energy balance equation, which takes the form ∇T · w, where w is the Darcy flux, and T is the temperature distribution. The first part of the paper is therefore concerned with analyzing a linearized version of the model, where we write the convective transport term as η · w, for some given η ∈ L ∞ (the remaining coupling terms are retained). Once we have obtained the existence and uniqueness of a weak solution to this problem, we introduce an iterative algorithm where we approximate the convective transport term as ∇T m−1 · w m , where m ≥ 1 is the iteration index. Due to the results we obtained for the linearized problem, and by a natural assumption that the temperature gradient admits L ∞ -regularity in space, we construct a well defined sequence of iterates as m → ∞. This we show to converge in adequate norms to the solution of the original nonlinear problem, thus establishing the existence and uniqueness of its weak solution. The convergence proof relies on the Banach Fixed Point Theorem, which we use to obtain local solutions in time. Here, the time interval is supposed to be small to ensure a contraction of the successive difference functions of the iterates. Then, using piecewise continuation in time, we extend these local solutions to global solutions for any finite final time. The idea is that such an iterative scheme can also be applied numerically to a discretized formulation, and in this sense our analysis sets the stage for subsequent numerical experiments. We mention also some of the literature on iterative schemes in poroelasticity; in [5,8,20,24] there can be found several iterative procedures for solving Biot's equations, and in [23,27,28] iterative methods for solving Richards' equation were analyzed.
We summarize the main contribution of the article as follows: under a natural hypothesis on the regularity of the convective term, we give a proof of existence and uniqueness of a weak solution to the fully coupled six-field thermo-poroelastic problem within the quasi-static framework.
The article is organized as follows: Section 2 recalls the physical model and the assumptions on the data, introduces the relevant function spaces and introduces the mixed weak formulations. In section 3 we define a linear version of the original mixed variational problem, and proceed to analyze this in the following way; we construct approximate solutions using a Galerkin method, the existence of which is established by the theory of DAEs (Differential Algebraic Equations). Suitable a priori estimates are then derived which enables us to pass to the limit, thanks to the weak compactness of the spaces. Section 4 is devoted to analyzing an iterative solution procedure for the original nonlinear problem and to establish the convergence of the algorithm in suitable norms. In Appendix A we propose an alternative to the hypothesis on the temperature gradient, i.e. we show how the required regularity can be obtained by sufficient regularity of the data. For easy reference of the notation used in this article we provide some tables in Appendix B.

Presentation of the problem
, be an open and bounded domain, where we denote the boundary by Γ := ∂Ω, which is assumed to be Lipschitz continuous. Let a time interval J = (0, T f ) be given with T f > 0 and define Q := Ω × (0, T f ] to be the space-time domain. The thermo-poroelastic model problem we consider, as it is exposed in [10], is as follows: given a heat source h, a body force f , and a mass source g, find (T, u, p) such that where a 0 is the effective thermal capacity, b 0 is the thermal dilation coefficient, β is the thermal stress coefficient, K = (K ij ) d i,j=1 is the permeability divided by fluid viscosity, Θ = (Θ ij ) d i,j=1 is the effective thermal conductivity, μ and λ are the Lamé parameters, α is the Biot-Willis constant and c 0 is the specific storage coefficient. The primary variables are the temperature distribution T , displacement u and fluid pressure p. To close the system, we prescribe homogeneous Dirichlet conditions on the boundary, i.e., T = 0, u = 0, and p = 0, on Γ × J, (2.1d) and we assume the following initial conditions for some known functions T 0 , u 0 and p 0 . In practice, we may use nonhomogeneous Dirichlet and Neumann boundary conditions for which the analysis remains valid. Note also that if β = b 0 = 0, the above system decouples from the energy equation, and the well-known quasi-static Biot equations are recovered (see e.g. [1] where both the two-and four-field formulations are presented).

Assumptions on the data
Before transcribing the mixed variational formulation of the problem (2.1), we make precise the assumptions on the data (further generalizations are possible, bringing more technicalities):

Mixed variational formulation
We now give the mixed variational formulation of the problem (2.1), for which we need to introduce the total stress tensor; σ(u, p, T ) := 2με(u) + λ∇ · uI − αpI − βT I, where I is the identity tensor and ε(u) is the linearized strain tensor given by ε(u) := (∇u + ∇ T u)/2, the Darcy flux w := −K∇p, and the heat flux r := −Θ∇T . For simplicity, we now restrict our attention to the case d = 2, in which case the fourth order compliance tensor, A, is given by as seen in [37] (see also [20] for the general formula). Note that A is bounded and symmetric positive definite uniformly with respect to x ∈ Ω, and defines an L 2 -equivalent norm, i.e.
where τ 2 A = Ω Aτ : τ dx. Applying A to the total stress tensor, it is inferred that and by taking the trace on both sides, we get the following relationship We also introduce the following notation The above definitions yields an equivalent mixed form to (2.1): We now set The following mixed variational formulation of the problem (2.1) can be obtained by multiplying by adequate test functions and then integrating by parts: and such that the initial conditions (2.1e) holds true in the weak sense, i.e.
Remark 2.1. Note that a different variational formulation of the problem (2.7) is possible, using a weakly symmetric space for the stress tensor. This formulation will then involve a new variable acting as a Lagrange multiplier which is enforcing the symmetry of the stress (see e.g. [2,4,20]). For simplicity of presentation we shall keep the formulation (2.8) throughout. The analysis presented next can nevertheless also be extended to the previously mentioned formulation using the same techniques, as done in [1] for the four-field Biot equations.
Remark 2.2. The nonlinear coupling in the above problem makes the analysis difficult. The next section is therefore devoted to analyzing a linearized problem, the results from which will be helpful when analyzing the full nonlinear problem in the last section. We mention also that other nonlinearities can be added, e.g. nonlinear compressibility or nonlinear Lamé parameters.

Analysis of the linear problem
In this section we introduce a linear version of the problem (2.8). Precisely, we replace the convective transport term (Θ −1 r · w, S) in the energy balance equation (2.8a), by −(η · w, S), for some given η ∈ L ∞ (Ω). We denote by γ := η ∞ . We introduce the resulting linear problem which reads: find and such that initial conditions (2.8g) holds true. The remaining part of this section is devoted to proving the well-posedness of this system. In what follows, we assume the following hypothesis on the effective thermal capacity a 0 , the thermal dilation coefficient b 0 , the specific storage coefficient c 0 and the Lamé parameters μ, λ; These constraints are typically needed in order to ensure a gradient flow structure. Similar constraints were used to analyze the Biot equations in mixed form in [1]. We also refer the reader to [21] for a more detailed discussion about the scaling of Biot's (isothermal) equations. However, compared to these works, our constraints involve also the thermal coefficients. We omit any further discussion on the justification for these constraints, other than they are necessary to prove the results we present. The well-posedness of problem (3.1) is then given in the following result.
The proof will follow from a series of partial results to be done in the sequel. The analysis uses a Galerkin method together with the theory of differential algebraic equations (DAEs), as well as weak compactness arguments (cf. [1,37,26,15]).

Construction of approximate solutions
In order to employ Galerkin's method we introduce a finite dimensional approximation of the problem (3.1). We need to introduce the following finite dimensional subspaces. Let (i, j, k, l, m, n) ∈ N 6 be fixed and strictly positive, and let T i := span{S ∈ T : = 1, · · · , i}, R j := span{y ∈ R : = 1, · · · , j}, P k := span{q ∈ P : = 1, · · · , k}, W l := span{z ∈ W : = 1, · · · , l}, S m := span{τ ∈ S : = 1, · · · , m} and U n := span{v ∈ U : = 1, · · · , n}, where the functions S , y , q , z , τ and v , for ∈ N, constitute Hilbert bases for the spaces T , R, P, W, S and U, respectively. Let now ( be the solution to the following problem: We introduce the coefficient vectors of the solutions: let Thus, we impose the initial conditions by We also define the following linear operators: Finally, we define the vectors: (L 1 ) := (f , v ), for 1 ≤ ≤ n, (L 2 ) := (g, q ), for 1 ≤ ≤ k and (L 3 ) := (h, S ), for 1 ≤ ≤ i. We rewrite using the above notation the problem (3.5) as a system of ODEs After rearranging, these ODE equations can be written in the form of a DAE system where X(t) := (P k (t), Σ m (t), T i (t), W l (t), U n (t), R j (t)) T , L(t) := (L 2 (t), 0, L 3 (t), 0, L 1 (t), 0) T and and From the theory of DAEs, equation (3.7) together with initial conditions (3.5g) has a solution if the matrix pencil, sΦ + Ψ, is nonsingular for some s = 0 (see [9]). Note that we can write sΦ + Ψ as a block 2 × 2 matrix as follows Let B = S m × P k × T i and C = U n × W l × R j , such that the bilinear form associated with sΦ + Ψ can be decomposed into the bilinear forms associated with each block, i.e. φ A :

A priori estimates
In this section, we derive a priori estimates for the unknowns which will allow us to pass to the limit in problem (3.5) by weak compactness arguments [38,15]. Throughout this section we denote by C > 0 a generic positive constant which may change value from one line to the next, but it will always be independent of the relevant parameters, i.e. of the tuple (i, j, k, l, m, n). We summarize these estimates in the following theorem.

Theorem 3.3 (A priori estimates). Under the Assumption 1, there exists a constant C > 0, independent of
Proof. By Thomas' Lemma [33] there exist σ ∈ H 1 (J; S m ) such that −∇ ·σ(·, t) = u n (·, t) on Ω for t ∈ J, and with σ(t) ≤ C u n (t) . Thus, we set τ =σ(t) in (3.5e) and obtain Next, we take τ = σ m in (3.5e) and v = u n in (3.5f), and add the resulting equations together to obtain Applying the C-S and Young inequalities together with the above estimate (3.14b) yields Choosing suitable values for the epsilons, i.e. 1 = 1 3α , 2 = 1 3β and 3 = 1 6C(μ + λ) , we obtain It then follows immediately that Take now σ ∈ L 2 (J; S m ) such that −∇ ·σ(·, t) = ∂ t u n (·, t) on Ω, for t ∈ J, and with σ(t) ≤ C ∂ t u n (t) . Then, by differentiating equation (3.5e) with respect to time, and setting τ =σ, we get in the same way as before We continue by differentiating equations (3.5e) and (3.5f) with respect to time, and take ∂ t σ m and ∂ t u n as test functions, respectively, and get analogously and Next, we take ∂ t σ m , p k , w l , T i and r j as a test functions in (3.5e), (3.5c), (3.5d), (3.5a) and (3.5b), respectively. We differentiate then (3.5f) with respect to time, and take u n as a test function. Adding together the resulting equations yields Using the properties of K and Θ, in addition to the C-S and Young inequalities yields (3.20b) Choosing = γ k m , integrating from 0 to t and substituting the inequalities (3.14b) and (3.15c), we deduce Since from (3.5g) we have we obtain the first estimate (i) using Grönwall's inequality, i.e.
For the second estimate, we differentiate (3.5e), (3.5f), (3.5d) and (3.5b) with respect to time and use ∂ t σ m , ∂ t u n , w l and r j as test functions, respectively. In (3.5c) and (3.5a), we use ∂ t p k and ∂ t T i as test functions, respectively. Summing the resulting equations yields By applying the C-S and Young inequalities, and substituting the estimates (3.17) and (3.18), we deduce Choosing suitable values for the epsilons, i.e. 1 = αβ C(μ + λ) , and 4 = αβ Simplifying the above expression, integrating over (0, t) and using the initial conditions yields (3.23d) It remains to provide estimates for w l (0) 2 and r j (0) 2 . To this end, take w l as a test function in equation (3.5d), and set t = 0. This gives which holds true for any k, l ≥ 1. Use now the properties of K to bound the left-hand side, tend k → ∞ and then integrate by parts in the right-hand side to obtain Thus, we have Similarly, using (3.5b), we obtain It remains to obtain the estimate (iv), for which we need just to bound the divergences. Since ∇ ·r j (t) ∈ L 2 (Ω) for t ∈ J, we can write ∇ · r j (t) = ∞ =1 ξ (t)S , for some functions ξ (t) ∈ R. Now, we multiply equation (3.5a) with ξ , sum over = 1, .., i and use the C-S and Young inequalities to obtain Using (3.18), integrating in time and using (3.25) we get It remains to tend i → ∞ to obtain Similarly, we obtain the following from equations (3.5c) and (3.5f) Combining the estimates (3.27c)-(3.27d) with (i) and (iii), we get the estimate (iv). This ends the proof. 2 The following estimates prove that the solution has improved regularity given some additional regularity on the data. We state the result as a lemma: Lemma 3.4 (Estimates for improved regularity). Assume that f ∈ H 2 (J; L 2 (Ω)) and g, h ∈ H 1 (J; L 2 (Ω)). Then there exists a constant C > 0 independent of (i, j, k, l, m, n) such that (i) Proof. We begin by differentiating equations (3.5e), (3.5c), (3.5d), (3.5a) and (3.5b) with respect to time, and take ∂ tt σ m , ∂ t p k , ∂ t w l , ∂ t T i and ∂ t r j as a test functions respectively. Then, we differentiate (3.5f) twice with respect to time, and take ∂ t u n as a test function. Summing the resulting equations yields Using the properties of K and Θ, in addition to the C-S and Young inequalities, we get By integrating over (0, t), using the initial conditions and substituting the inequalities (3.18) and (3.19), it is inferred that We proceed to bound ∂ t p k (0) and ∂ t T i (0) . To this end, we discard the terms under the time differential on the left-hand side of (3.23c) and set t = 0 to obtain We use (3.24c) to bound the initial value of the Darcy flux, i.e., Now we substitute this in (3.28c), using also (i) from Theorem 3.3 and apply Grönwall's Lemma to obtain Summing the above estimate with estimate (iii) from Theorem 3.3 produces the estimate (ii). Going back to the estimate (3.27a), we now substitute in the right-hand side with (3.30) and (3.31), let i → ∞ to obtain From equations (3.5c) and (3.5f) we obtain using the same technique and Summing the estimates (3.32)-(3.34) and combining with (ii) and (iii) from Theorem 3.3 produces the estimate (iii). This ends the proof. 2

End of the proof of Theorem 3.1:
The proof of the first part of Theorem 3.1 follows the steps below: , and {r j } ∞ 0 is bounded in L 2 (J; H(div, Ω)) ∩ L ∞ (J; L 2 (Ω)). By the weak compactness properties of the spaces there exist subsequences (denoted the same way as before) and functions σ ∈ L ∞ (J; H s (div, Ω)) ∩ H 1 (J; L 2 (Ω)), u ∈ H 1 (J; L 2 (Ω)), p ∈ H 1 (J; L 2 (Ω)), w ∈ L 2 (J; H(div, Ω)) ∩ L ∞ (J; L 2 (Ω)), T ∈ H 1 (J; L 2 (Ω)), and r ∈ L 2 (J; H(div, Ω)) ∩ L ∞ (J; L 2 (Ω)), such that In order to pass to the limit in problem (3.5), we fix a tuple (i, j, k, l, m, n) ≥ 1 and take (S, y, q, z, τ , v) ∈ C 1 (J; T i ×R j ×P k ×W l ×S m ×U n ) as test functions, and then integrate equations (3.5a)-(3.5f) with respect to time to obtain Passing to the limit yields Finally, by the density of the test function space, C 1 (J; T i × R j × P k × W l × S m × U n ) in L 2 (J; T × R × P × W × S × U) as (i, j, k, l, m, n) → ∞, the equations (3.1) hold true for a.e. t ∈ J. It remains now to show that the initial conditions are satisfied, i.e. T (0) = T 0 , u(0) = u 0 and p(0) = p 0 , in the weak sense. To this end, take q ∈ C 1 (J; P k ) such that q(T f ) = 0 as a test function in (3.35c) and integrate the first term by parts in time On the other hand, from (3.36c) we obtain Since q(0) was arbitrary, and since p n (0) → p 0 in L 2 (Ω), we get that p(0) = p 0 . We obtain in the same way that u(0) = u 0 , and T (0) = T 0 .
• To finish the proof we show the uniqueness of a weak solution to problem (3.1). To this end, we assume that (T 1 (t), r 1 (t), p 1 (t), w 1 (t), σ 1 (t), u 1 (t)) and (T 2 (t), r 2 (t), p 2 (t), w 2 (t), σ 2 (t), u 2 (t)) are two solution tuples in T × R × P × W × S × U, and let (e T (t), e r (t), e p (t), e w (t), e σ (t), e u (t)) be the corresponding difference. This then satisfies the following variational problem: find (e T (t), e r (t), e p (t), e w (t), e σ (t), e u (t)) ∈ T × R × P × W × S × U such that for a.e. t ∈ J there holds Integrating the above equation from 0 to t and using the properties of K and Θ, in addition to the C-S and Young inequalities yields which after application of the Grönwall inequality yields Then, using Thomas' Lemma [33] we take τ =σ(·, t) ∈ S in (3.39e), such that for t ∈ J, −∇ ·σ(t) = e u (t) in Ω, with σ(t) ≤ C e u (t) for some constant C > 0. Thus, we obtain =⇒ e u ≤ C( e σ + e p + e T ). (3.47) This implies that e T (t) = e r (t) = e p (t) = e w (t) = e σ (t) = e u (t) = 0, in Ω, for a.e. t ∈ J, implying the uniqueness of a weak solution to problem (3.1). Finally, thanks to Lemma 3.4, we can finish the proof of the second part of Theorem 3.1 using similar arguments. 2

Analysis of the non-linear problem
We now consider the analysis of the mixed variational formulation for the original nonlinear problem (2.8). The analysis uses the results derived previously for the linear case, in addition to the Banach Fixed Point Theorem (see e.g. [11]) in order to obtain a local solution to (2.8) in time. We then proceed to extend this local solution by small increments until a global solution is obtained for any finite final time (see e.g. [19,35] where similar techniques are used). Precisely, an iterative solution procedure is introduced based on linearizing the heat flux term in (2.8a), which is shown to be well-defined, and which converges to the weak solution of the nonlinear problem in adequate norms. Note that we now must require the iterates to be continuous in time, hence we shall invoke Lemma 3.4. The iterative linearization algorithm we consider is then as follows: let m ≥ 1, and at the iteration m, we solve for (T m , r m , p m , w m , σ m , u m ) ∈ T × R × P × W × S × U such that for t ∈ J there holds , ∀q ∈ P, (4.1c) together with initial conditions, (2.8g), and where the algorithm is initialized by given initial guess r 0 . We consider the following hypothesis on the heat flux: Hypothesis 1 (The heat flux). We suppose that for all m ≥ 1, the heat flux is such that r m (t) ∈ L ∞ (Ω), for t ∈ J.
The above hypothesis is a natural one, and it is necessary for the solution to the iterative procedure (4.1) to be well-defined for each m ≥ 1. This hypothesis is satisfied with sufficiently regular data and domain boundary. We provide some formal arguments in Appendix A on the specific requirements such that the solution to the problem (3.1) yields r ∈ C([0, T f ], L ∞ (Ω)) (or alternatively w, r ∈ C([0, T f ]; L 4 (Ω))), thus making the above hypothesis superfluous. We delegate this discussion to the Appendix in order to avoid overly strict assumptions on the data. Remark 4.1. Note that if we had approximated the convective term in equation (4.1a) instead as (w m−1 · Θ −1 r m , S), Hypothesis 1 would be on the regularity of the Darcy flux w, and the above algorithm would be initialized by some w 0 . However the analysis presented next remains true and follows exactly the same lines.
Based on the development of the previous sections, we now state the main result of this article.  Using Thomas' lemma [33], we take τ =σ(·, t) in equation (4.4e) such that e m u (·, t) = ∇ ·σ(·, t) with σ(t) ≤ C e m u (t) for t ∈ J, and combine with (4.5b) to obtain Therefrom, we proceed as on done in the first time interval [0, t 1 ] to show the convergence of the successive approximations (T m , r m , p m , w m , σ m , u m )| [t k−1 ,t k ] , k ∈ N, to (T, r, p, w, σ, u)| [t k−1 ,t k ] . This solution is similarly extended to any time t ≥ t k given by Finally, since the series Some remarks on the above proof are in order.

Remark 4.2.
We could also define a fully explicit iterative scheme where both the Darcy and heat fluxes in the convective term are given at the previous iteration. If such an explicit scheme was chosen we would have the advantage of a symmetric linearized problem, as the convective terms in the iterative procedure can be viewed as part of the source term on the right hand side. Remark 4.3. Assume that f is in H 1 (J; L 2 (Ω)), g, h in L 2 (J; L 2 (Ω)), p 0 , T 0 in H 1 0 (Ω), and u 0 in (L 2 (Ω)) d . Suppose that instead of Hypothesis 1, we have r m , w m in ∈ H 1 (0, T ; L ∞ (Ω)). Then, we can reproduce the proof of Theorem 4.1 to prove the convergence of the scheme given by (4.1) to a weak solution of the nonlinear problem (2.8).

Conclusions
In this article we have given mixed formulations for the fully coupled quasi-static thermo-poroelastic model. The model is nonlinear, with the nonlinearity appearing on a coupling term. This makes the analysis challenging. A linearization of the model was therefore employed as an intermediate step in analyzing the full nonlinear model. For the linear case, the well-posedness is established using the theory of DAEs, and energy estimates together with a Galerkin method. This result together with derived energy estimates are combined with the Banach Fixed Point Theorem to obtain local solutions in time of the nonlinear problem. Due to the continuity in time of the convergent (local) solutions, we can infer a (global) convergence proof of an iterative procedure approximating the weak solution to the original nonlinear problem. Work underway addresses discretization of this model problem using an appropriate mixed finite element method as well as a priori and a posteriori error analysis.