Numerical modelling of convection-driven cooling, deformation and fracturing of thermo-poroelastic media

Convection-driven cooling in porous media influences thermo-poro-mechanical stresses, thereby causing deformation. These processes are strongly influenced by the presence of fractures, which dominate flow and heat transfer. At the same time, the fractures deform and propagate in response to changes in the stress state. Mathematically, the model governing the physics is tightly coupled and must account for the strong discontinuities introduced by the fractures. Over the last decade, and motivated by a number of porous media applications, research into such coupled models has advanced modelling of processes in porous media substantially. Building on this effort, this work presents a novel model that couples flow, heat transfer, deformation, and propagation of fractures with flow, heat transfer, and thermo-poroelasticity in the matrix. The model is based on explicit representation of fractures in the porous medium, and discretised using multi-point finite volume methods. Frictional contact and non-penetration conditions for the fractures are handled through active set methods, while a propagation criterion based on stress intensity factors governs fracture extension. Considering both forced and natural convection processes, the numerical results show the intricate nature of thermo-poromechanical fracture deformation and propagation.


Introduction
For a porous medium, possibly containing fractures, the interplay between flow, thermal transport, and deformation can be strong. In particular, cooling of the medium induces thermal stress that can lead to deformation and fracturing. Furthermore, fractures deform and propagate as a result of the coupled dynamics. The result is coupled thermo-hydro-mechanical (THM) processes in the intact porous medium, interacting with flow and thermal transport in fractures as well as fracture deformation and propagation. Such coupled process-structure interaction is characteristic for a wide range of natural and engineered processes in natural and manufactured materials. For example, the structural and functional performance of concrete structures, like dams, bridges, nuclear and liquefied natural gas containers, and cement sheaths of subsurface well bore constructions, are affected by the time evolution of their properties under variable THM loads [5,14,44,18,19,45]. In the subsurface, THM processes interact with deformation and propagation of fractures in fluid injection operations [58,32,64]. The coupled dynamics is also hypothesised to be crucial in heat transfer from the deep roots of geothermal systems by deepening natural convection through evolving fractures [46,13,11,12]. Common to all these applications is that tight coupling in the dynamics limits the knowledge which can be gained from analysis of individual processes and mechanisms in isolation. This motivates development of simulation models that acknowledge the coupled nature of the physics.
Since its foundation by Biot [10], the theory of poroelasticity has successfully been applied to model coupled hydro-mechanical processes. The extention to thermo-poroelasticity [23] is also widely applied, including in geomechanics [56]. More recently, models accounting for discontinuities in the form of fractures in poroelastic and thermo-poroelastic media have been developed. Typically, the development has focused either on deformation of preexisting fractures or the mechanical fracturing of the materials. The models can be distinguished based on whether fractures are represented explicitly as discrete objects embedded in the porous medium, or represented as part of the porous medium itself. The latter incorporate the effect of the extent to which the material is fractured by use of smeared or distributed representations. Such models include phase-field and damage approaches for fracture [16] and continuum and multi-continuum approaches for flow models [9].
Approaches based on explicit representation of the fractures can further be distinguished by how the fractures are represented in discretisation, specifically on whether a conforming or non-conforming representation of the fractures is used in the grid [9]. Non-conforming methods represent the fracture through an enriched representation. For poromechanics, combinations of the embedded discrete fracture method, extended finite element methods and/or embedded finite element methods have been applied [52,24,35]. Such non-conforming approaches have also been extended to include tensile fracture propagation based on extended finite element [42] and embedded discrete fracture methods [27]. Conforming methods use a representation where the fractures coincide with matrix faces. Considering fractures that have a negligible aperture compared to the modelled domain, this representation can be combined with an approach where fractures are modelled as lower-dimensional structures [47,39] and discretised with elements of zero thickness [43,15,28,8]. For poroelastic media with fractures, this allows for the application of standard finite element [55], finite volume [62,7] and combined finite element/finite volume schemes [34,57,33], more recently also including fracture contact mechanics [34,31,30] and tensile fracture propagation [57,55]. Thermal effects on fracture deformation and propagation in poroelastic media are less studied, although some recent studies model deformation of existing fractures in thermo-poroelastic media [54,60,33].
Motivated by the development of increasingly sophisticated THM models for fractured porous media, our goal in the present paper is to extend numerical modeling of THM to also incorporate fracture propagation, and thereby contribute to bridge the gap between fracture mechanics models and coupled THM models for porous media. Specifically, we consider mathematical and numerical modelling of fracture deformation and propagation resulting from coupled THM-processes. Our focus is on convection-driven cooling in the subsurface, where forced or natural fluid convection induces thermo-poromechanical stress changes leading to fracture deformation and propagation. The dynamics is characterised by tight coupling between physical processes and strong interaction between the physical processes and the (evolving) geometry of the fracture network. Accordingly, our model and simulation approach is designed to faithfully represent these couplings, including fracture deformation and propagation. The fractured thermo-poroelastic medium is represented using a discrete fracture-matrix model, where fractures are represented as lower-dimensional discontinuities in an otherwise continuous thermo-poroelastic medium. Deformation of existing fractures is modelled through contact mechanics relations based on a Coulomb friction criterion for slip along the fractures and a nonpenetration condition [37,7]. This is combined with a simple criterion for fracture propagation based on the mode I stress intensity factor, which we compute directly from the displacement jump in the vicinity of the fracture tip using a variant [48] of the displacement correlation method [21]. To adjust the grid to an arbitrary fracture propagation path is highly technical [51,25], and we instead make the assumption that fractures propagate along existing faces in the matrix grid. This constrains the numerical representation of an evolving fracture and makes it difficult to preserve reasonable fracture geometries for general propagation scenarios, in particular for three-dimensional problems. We therefore further limit ourselves to tensile fracturing, where the possible propagation path is easy to predict and the grid can be constructed to accommodate the propagation. We discretise the model using a control volume framework for fracture contact mechanics in thermo-poroelastic media [60]. The control volume approach builds on a combination of the multi-point stress approximation method for Biot poroelasticity [49,41] with the multi-point flux approximation method for flow [2]. This combination is previously applied for numerical modelling of fractured poroelastic media [62] with a simplified model for deformation along fractures. The fracture contact mechanics builds on work by Berge et al. [7], who formulated the contact conditions on the fracture using Lagrange multipliers representing the contact tractions [63]. Using this approach, the variational inequality representing the contact problem can be rewritten using complementary functions, and the resulting system of equations solved by a semi-smooth Newton method [37,7]. Our model is implemented in the open-source simulator PorePy [40], which is designed for multiphysics problems in fractured porous media.
We assess the reliability of our simulation tool by tests that probe the approximations of both the onset of fracturing and the speed of fracture propagation. We then present two application-related simulations that both involve fracture propagation driven by convective cooling. The cases include respectively forced convection during production of geothermal energy and natural convection in vertical fractures in the presence of high thermal gradients. Taken together, the results show the importance of developing simulation tools that can accurately represent the tight couplings in THM processes, and also deal with deformation and propagation of fractures.
The paper is structured as follows. Section 2 presents the governing model equations for poroelastic media with deforming and propagating fractures. The discretisation schemes and numerical solution strategy is presented in Section 3. Section 4 presents simulation results, before concluding remarks are given in Section 5.

Governing equations
The conceptual model is based on explicit and conforming representation of fractures in the porous medium. Two modes of fracture deformation are considered: Deformation with fixed transverse extension governed by contact mechanics relations and deformation through irreversible fracture propagation. We also impose conservation of mass and energy in matrix and fractures and momentum balance in the matrix.

Geometrical representation of fractured porous media
The model and governing equations are posed in a mixed-dimensional framework arising from considering fractures as lower-dimensional objects. Hence, in a three-dimensional domain, fractures are represented as two-dimensional surfaces, and in a two-dimansional domain, they are one-dimensional lines. In a D-dimensional domain, we denote the matrix subdomain by Ω h and fractures are represented by subdomains Ω l of dimension D −1. The matrix and fractures are connected by interfaces denoted by Γ j , with the subscript pair j, k used to indicate the two interfaces on either side of a fracture, see Fig. 1. The boundary of Ω i is denoted by ∂Ω i , and the internal part of it corresponding to Γ j is ∂ j Ω i .
We also use subscripts i, h and l to identify the domain of the primary variables, which are displacement, pressure, temperature, contact traction (u, p, T and λ). Similarly, subscript j denotes the four interface variables defined in Sections 2.2 and 2.6. The subscripts are suppressed when context allows, as are the subscripts f and s denoting fluid and solid, respectively.
To model fracture deformation, it is necessary to decompose a vector into its normal and tangential components relative to a fracture. The fracture normal is defined to equal the outwards normal n h on the j side, i.e. n l = n h | ∂j Ω h . A vector ι l may now be decomposed as ι n = ι l · n l and ι τ = ι l − i n n l , where subscripts n and τ denote the normal and tangential direction, respectively.

Contact mechanics for fracture slip and opening
The contact mechanics relations are a traction balance between the two fracture surfaces and a nonpenetration condition, complemented by a Coulomb friction law governing the relative displacement when the surfaces are in contact. These relations are formulated in the displacement jump [[u]] and the contact traction λ l . The higher-dimensional THM traction, σ h · n h , is balanced by the contact traction and the fracture pressure on the two interfaces: Here the notation indicating that the variable is taken at the interface Gamma j or Γ k should be interpreted as the extension and projection of this variable to the respecitve interface. The displacement jump over the fracture is defined as with u j and u k denoting displacement at Γ j and Γ k , cf. Fig. 1. The gap function g is defined as the normal distance between the fracture surfaces when these are in mechanical contact. Following Stefansson [60], we set with ψ denoting the dilation angle [6], thus accounting for shear dilation of the fracture resulting from tangential displacement [[u]] τ of the rough fracture surfaces.
Given that fracture surface interpenetration and positive normal contact traction are prohibited, the following conditions have to be fulfilled: Hence, when a fracture is mechanically open and there is no mechanical contact across the fracture, the normal contact force, λ n , is zero. The friction law is imposed by enforcing where F and [[u ]] τ denote the friction coefficient and the tangential (shear) displacement increment, respectively. For simplicity, we consider a constant coefficient of friction in this work.

Fracture propagation
Fracture propagation occurs when the potential energy released by the extension exceeds the energy required to separate the fracture surfaces by breaking atomic bonds [36]. Using concepts of linear elastic fracture mechanics, we evaluate propagation based on computation of stress intensity factors (SIFs). The SIFs are computed directly from the displacement jump in the vicinity of the fracture tip using the variant of the displacement correlation method [48]. Referring to Fig. 1, the local geometry at a fracture tip is described using a coordinate system given by the orthogonal basis vectors e ⊥ , e n and e (or e ⊥ and e n if D = 2). We set e n = n l , and the tangential (τ ) vectors e ⊥ and e are respectively perpendicular and parallel to ∂Ω l at ∂Ω l (see Fig. 1). By use of the components of the displacement jump in the local coordinate system, the displacement correlation method gives the three SIFs Here, µ denotes the shear modulus and κ = 3 − 4ν the Kolosov constant, with ν being the Poisson ratio. R d is the distance from the fracture tip to the point at which the displacement jump is evaluated. The three stress intensity factors are related to tensile (K I ), shear (K II ) and torsional (K III ) forces. As stated in the introduction, we limit ourselves to tensile fracture in this work, and ignore contributions from (K II ) and (K III ). With this assumption, a tip propagates if the computed mode I factor exceeds a critical value, and the propagation angle φ illustrated in Fig. 1 is zero. Criteria for more sophisticated mixed-mode propagation, which can be highly relevant for subsurface applications, are reviewed by Richard et al. [53].

Fracture mass and energy balance
The thickness of a dimensionally reduced fracture is represented by the aperture, which changes as the domain deforms according to with a res denoting the residual hydraulic aperture in the undeformed state representing the effect of small-scale roughness of the two fracture surfaces. The tangential fracture permeability K l is chosen to depend on aperture by the nonlinear relationship K l = a 2 /12, which corresponds to setting the hydraulic aperture of the fracture equal to a [65].
On the assumption that the fractures are completely filled with fluid, the parameters of this subsection equal that of the fluid. We assume single-phase flow according to Darcy's law: where µ, ρ and g denote viscosity, density and the gravity acceleration. The total heat flux may be split into continuum scale heat diffusion modelled by Fourier's law and advection along the fluid flow field: where κ and C denote thermal conductivity and heat capacity, respectively. Following Stefansson et al. [60] (see also Brun et al. [20] and Coussy [23]), balance of mass for a fracture Ω l reads with c, β and q p denoting compressibility, thermal expansion coefficient and a fluid source or sink term. Next, assuming local thermal equilibrium between fluid and solid, neglecting viscous dissipation and linearising [60], the energy balance is where we assume thermal sources and sinks to satisfy aq T = aq p Cρ T0 (T − T 0 ) and T 0 denotes a reference temperature. In Eqs. (12) and (13) the last terms on the right hand sides represents the fluxes from matrix to fractures, which are defined in Section 2.6.
In deriving these equations, the following equations of state are assumed [23] for density and entropy

Matrix thermo-poroelasticity, energy and mass balance
The following section presents the balance equations and constitutive relations for the matrix problem. The model resembles that of the previous section, with the addition of a momentum balance equation for the thermo-poroelastic medium, yielding three balance equations for Ω h . For details on the derivations of the equations, we again refer to Coussy [23] and Brun et al. [20]. We first define the following effective parameters [22], arising through the assumption of local thermal equilibrium: φ and α denote porosity and the Biot coefficient, respectively. Neglecting inertial terms, the momentum balance is with q u denoting body forces and the linearly thermo-poroelastic stress tensor related to the primary variables by an extended Hooke's law The mass balance equation reads while the energy balance is On Ω h ∪ Γ j , the following internal boundary conditions ensure coupling from Ω h to the interface variables on Γ j : The conservation equations are complemented by appropriate boundary conditions on the domain bondary. This applies to both the matrix and fracture domains.

Interface fluxes between fractures and matrix
Interface flux relations close the mixed-dimensional system of mass and energy balance equations [47,38]: We set the normal permeability and thermal conductivity equal to their tangential counterparts, i.e. K j = K l and κ j = κ l .

Discretisation and solution strategy
Discretisation of the governing equations entails devising discrete representation of the conservation equations and of the contact mechanics relations on existing fractures. Moreover, when the propagation criteria are met, the fracture geometry must be modified and the discretisations updated accordingly. We make the following assumptions on the computational grid: Grids for the subdomains Ω h and Ω l and the interface Γ j are constructed so that faces on ∂ j Ω i match with cells in Γ j and Ω l . We make no assumptions on the cell types; for the simulations presented in Section 4 we mainly use Cartesian grids as these are most easily fit to a known, straight propagation path, but also consider simplex cells for one simulation.

Spatial discretisation
Pressure and temperature are represented by their cell centre values in Ω h and Ω l , as is the displacement in Ω h and contact force in Ω l . The discrete primary variables on Γ are displacements, mass flux and advective and diffusive heat fluxes.

Contact mechanics
The non-linear contact mechanics problem is represented by an active set approach implemented as a semi-smooth Newton method following [37,7]. The treatment of Eqs. (5) and (6) depends on whether the previous iterates were in an open, sticking or gliding state, with the states evaluated cell-wise in Ω l . Equation (2) is discretised by relating the cell centre pressures and contact force in Ω l to the discrete traction on Γ.

Discretisation of balance equations
For Ω h , the stress term in (17) and the diffusive fluxes in (19) and (20) are all discretised with a family of finite volume multi-point approximations termed MPxA [2,49,50]. The methods construct discrete representations of the constitutive relations, Hook's, Darcy's and Fourier's law, in terms of the cell centre variables. These relations are used to enforce conservation of THM traction, mass and (diffusive) heat flux over the cell faces. For faces on the fracture surfaces, the discrete traction enters the contact mechanics discretisation described above. The full heat flux is given by the sum of the discrete Fourier's law and the advective flux, where the latter is discretised by a single-point upstream method. For further information on the MPxA methods, we refer to [50].

Solution of non-linear system
The discretised system of equations is solved by Newton's method, with the terms from the contact conditions handled by a semi-smooth approach following [37,7]. The termination criterion for the Newton iterations considers the residuals and updates of each of the primary variables u h , u j , p and T. Within each non-linear iteration, the linearised system is solved using a direct sparse solver [26]. While simple, this approach is memory intensive and puts practical constraints on mesh resolution, in particular for three-dimensional problems. A more scalable method would involve iterative solvers with block preconditioners for the THM components of the linear system [17], with a tailored treatment of the contact conditions [29].

Solution algorithm
The temporal derivatives are discretised by a backward Euler scheme, and the THM contact mechanics problem is solved monolithically, using implicit in time evaluation of all spatial derivatives. When the non-linear solver has converged, we proceed to fracture propagation evaluation.
Stress intensity factors and the fracture propagation criterion are evaluated for each fracture tip faces using Eqs. (7) and (8). The displacement jump is evaluated at the neighbouring cell of the tip face, i.e. R d is the distance between Variables are initialised in the new cells using the reference values p 0 and T 0 , and new apertures are set to a res . This in effect adds mass to the system, cf. Eq. (14). To compensate, we prescribe an additional term on the right hand side of equation (12) equal to −a res /dt in newly formed fracture cells the subsequent time step, with dt denoting time step size. Since Eq. (13) is derived by considering s − s 0 , Eq. (15) implies that no right-hand side term arises with the chosen initialisation values.
Before the simulation proceeds to the next time step, all terms are rediscretised to account for modifications of the grids. This can be done locally, i.e. only for the faces and cells where the discretisation is affected by the grid update.

Simulation results
The results presented in this section serve to first verify the computational approach, and then to show application to two subsurface cases involving THM processes and fracture propagation. The PorePy toolbox [40,1] was used for all simulations and run scripts for geometry and parameter setup etc. are available on GitHub [61]. All parameters not specified in the text are listed in Table 1.

Verification
The verification of the computational approach entails first a test of the numerical stress intensity factors and next a convergence test of the fracture propagation speed.

Example 1: Stress intensity factors
To verify the SIF computation, we consider an analytical solution, first derived by Sneddon [59], for a single crack in an infinite medium with uniform internal pressure on the fracture surfaces. Boundary conditions for the finite simulation domain are computed using the boundary element method following Keilegavlen et al. [40], who also presents a thorough convergence study for the aperture using PorePy. Herein, we compare the SIFs as computed by the displacement correlation method to the analytical solution Here, p f denotes the internal pressure on the fracture and l denotes fracture length. We use a square domain of side length 50 m, l = 10 m and p f = 1 × 10 −4 Pa. We consider a sequence of four grids, the finest of which is shown in Fig. 3. To probe the method for different material parameters, we also use four different Poisson ratios. Based on displacement solutions on each with the j index running over the two fracture tips. The piecewise linear displacement representation of the MPSA discretisation does not capture the stress singularity at the fracture tips. Since the SIFs are computed from [[u]] in these very tip cells, the method does not converge with mesh refinement. Rather, the results presented in Fig. 4 demonstrate robustness with respect to mesh size and the Poisson ratio ν. While we do not consider K II in the subsequent simulations, we also present results demonstrating that the method indeed predicts tensile stresses (i.e. K II K I ) for this purely tensile problem. The results of this test indicate that the MPSA solution can be used to estimate K I in tensile problems, and thus form the basis of fracture growth evaluation.

Example 2: Propagation speed
We now consider a test case designed to evaluate the simulated propagation speed of a fracture in a tensile regime of stable propagation. The unit square domain contains two horizontal fractures Ω 2 and Ω 3 extending 1/4 from the left and right boundary, respectively, see We use four temporal refinement levels and three spatial refinement levels in addition to a highly refined reference solution. The finest (non-reference) mesh and the final fracture geometry are illustrated in Fig. 5. Figure 6 shows fracture size plotted against time for all refinement combinations. With one exception discussed below, the results group according to spatial resolution. While the propagation speed is fairly constant across all mesh sizes, propagation onset occurs earlier for the coarser meshes. We attribute this offset to the SIFs being evaluated on the basis of [[u]] at the centre of the fracture tip cell. The location of this cell centre is closer to the boundary for the coarser meshes, implying shorter travel time for the cooling front. As expected, the plot indicates convergence with mesh refinement.
The outlier is the smallest mesh size combined with the largest time step, for which the propagation speed is notably lower. The propagation speed is simply not resolved by the spatio-temporal discretisation, i.e. the propagation speed exceeds h/dt. In other words: Given a spatial resolution, an upper bound on the time step must be honoured in the explicit type of propagation solution algorithm used herein.

Applications
We present two simulations that involve THM processes coupled with fracture propagation. The first case resembles geothermal energy production, with convection forced by fluid injection and production. The second case involves natural convection that takes place mainly inside fractures. In both cases, convection   acts to alter thermal stresses and thereby cause fracture propagation.

Example 3: Thermal fracturing and forced convection
We consider two immersed fractures in a cube shaped domain of side length 1000 m centred 1500 m below the surface. Each fracture contains one well, implemented as a source or sink term in a single cell, with injection in the leftmost fracture, Ω 2 , and production in the rightmost fracture, Ω 3 . The domain, fracture geometry and spatial mesh is shown in Fig. 7. The flow rate is 5 L s −1 for both wells and the injection temperature is 30 K below the formation temperature. The anisotropic boundary tractions are based on lithostatic stress, with This background stress implies that Ω 3 , with normal vector n 3 = [1, 0, 0] T , is initially more favourably oriented for propagation than is Ω 2 . The Fig. 9 fracture size plot shows that Ω 2 grows at a steady speed after an initial phase of limited propagation. The growth is driven by elevated pressure due to injection and matrix cooling, which is most pronounced on the side of Ω 2 facing Ω 3 due to the advective component of the heat flow, cf. Fig. 8. Assuming the thermal driving force to dominate, which is reasonable given the relative size of injection pressure and background stresses, the relatively constant speed could be linked to the constant rate and temperature of injection. The fracture Ω 3 , where fluid is produced, does not propagate at all. Towards the end of the 2.5 yr simulation, the magnitude of normal traction on Ω 3 has increased considerably relative to the initial value of approx. 1 × 10 7 Pa, see Fig. 9. We attribute this to the contraction ensuing from matrix cooling surrounding Ω 2 , which leads to a larger proportion of the compressive forces being supported by the non-cooled surroundings, including Ω 3 . Figure 9 also shows temperature and pressure in the two wells throughout the simulation. Most notably, injection pressure gradually declines. This increased injectivity in Ω 2 is caused by the combination of an increased aperture in the pre-existing part of the fracture, and the increase in the geometric extension of the fracture. Thus, fracture deformation caused by thermal and hydraulic stimulation strongly affects the (flow) properties, providing a clear example of the two-way process-structure interaction characteristic of fractured porous media.
This simulation indicates that long-term cooling during geothermal energy production may alter the stress state to a stage where fractures propagate. It is thus important to develop simulation tools that can incorporate such changes to fracture geometry, in addition to handling multiphysics processes in the reservoir. Moreover, the injection pressure evolution shows the importance of also capturing deformation of existing fractures in the same model.

Example 4: Thermal fracturing and natural convection
As a final example, we consider fracture propagation driven by cooling that is mainly caused by convection cells inside vertical fractures. The process, known as convective downward migration, has been proposed as a mechanism for trans- port of heat in the deep roots of volcanic geothermal systems [46,13]. It is also predicted to have an important role in the source mechanism of hydrothermal activity in a more general perspective [13,4]. We consider five vertical fractures evenly spaced along the x direction and extending from the top boundary half-way through the cube-shaped domain with side length 400 m, see Fig. 10. The domain is centred 2800 m below the surface to mimic conditions within the earth crust where natural heat convection is likely to take place. Boundary and initial conditions are hydrostatic pressure and temperature according to a vertical gradient of −0.15 K m −1 and upper boundary temperature 500 K, considered to represent background temperature gradient close to the boundaries of geothermal areas. This value is estimated between −0.10 and −0.15 K m −1 within Iceland's active zone of volcanism and rifting, where many high temperature systems exist [3]. The boundary traction is the same as in the previous example and the simulation time is 70 years. The results are displayed in Figures 10, 11 and 12.
The vertical temperature gradient leads to instabilities in fluid density, which triggers convection cells inside the fracture, see Fig. 12. As shown by the temperature contour surfaces in Fig. 11, the resulting energy transport cools the rock surrounding the fractures, to the point where propagation occurs at the lower end of the fractures. This change in fracture geometry, together with changes in aperture in the existing fracture due to contraction of the surrounding rock, again gives feedback to the fluid convection, as is evident from the difference in flow patterns between the solutions at the two different times reported in Fig. 12. As in Example 3, we see evidence of tight process-structure interaction, with the convection-induced cooling altering thermo-poro-mechanical stress sufficiently for the fractures to open and propagate. Figure 10 displays size evolution for individual fractures. Propagation begins approximately half-way through the simulation, first for the fracture in the center of the domain. Even after all fractures have started propagating, the central fractures Ω 3 , Ω 4 and Ω 5 propagate significantly faster than the two outermost. This should be understood in the context of the compressive boundary conditions: The normal tractions on Ω 2 and Ω 6 , respectively, are not relieved by the cooling of any fractures lying between them and the left and right boundary.
After onset, propagation continues until the end of the simulation, but not at all time steps for all propagating fractures, and certainly not along the entire propagation front. This is because the matrix surrounding the new part of a fracture must be cooled before the fracture proceeds, and indicates that the fracture growth is stable as in Example 2 and that the propagation speed is resolved in the temporal discretisation. An approximate downward propagation speed for fractures 3-5 is obtained by dividing the estimated slopes from Fig. 10 by the initial lateral fracture length 200 m, yielding ∼ 2 m yr −1 . The setup for this test case is based on average properties in high temperature settings, and the results are in agreement with previous assessments of 0.3 m yr −1 to 5 m yr −1 [13,11], using a simple relation between the temperature difference sufficient for thermal stress to outweigh the hydrostatic force to keep the fracture closed, at approximatly 3 km depth in the crust with average properties of water and

Conclusion
While both numerical models considering flow and heat transfer in fractured media and models considering deformation of poroelastic media and fracture mechanics have separately been studied extensively, models which combine these fields are more recent. In the current work, we present a novel numerical model that couples fracture contact mechanics and propagation with deformation, flow and heat transfer in fractured thermo-poroelastic media. The methodology is built on a multi-point control-volume framework, combined with an active-set approach for fracture contact mechanics. The fracture propagation is based on stress intensity factors, and computed using a variant of the displacements correlation method. In the numerical model, fractures are restricted to propagate conforming to the existing grid. The numerical results show mesh convergence for computation of stress-intensity factors and fracture propagation speeds. Focusing on tensile fracture propagation, three-dimensional numerical test cases also show how the model can be used to investigate fracture propagation caused by forced and natural convection, exemplified by long-term thermal reservoir stimulation due to cooling and convective downward migration of fractures. The simulations demonstrate the need for coupled models accounting for both contact mechanics and fracture propagation as well as the coupled thermo-poroelasticity. Figure 11: Example 4: Solution and fracture geometry at the end of the simulation. p, a and λ n are shown on the fractures, while T is shown both on the fractures and as contour surfaces indicating where the matrix is significantly cooled (7.5 K and 15 K below initial formation temperature). The top and bottom row correspond to just before propagation onset (t = 32 yr) and the end of simulation (t = 70 yr), respectively. The scaling of the arrows indicating flux direction and magnitude is a factor five larger for the top row, when fluxes are smaller due to smaller apertures. The red rectangles indicate initial fracture geometry.