A fully coupled numerical model of thermo-hydro-mechanical processes and fracture contact mechanics in porous media

A range of phenomena in the subsurface is characterised by the interplay between coupled thermal, hydraulic and mechanical processes and deforming structures such as fractures. Modelling subsurface dynamics can provide valuable phenomenological understanding, but requires models which faithfully represent the dynamics involved; these models, therefore are themselves highly complex. This paper presents a mixed-dimensional thermo-hydro-mechanical model designed to capture the process-structure interplay using a discrete-fracture-matrix framework. It incorporates tightly coupled thermo-hydro-mechanical processes based on laws for momentum, mass and entropy in subdomains representing the matrix and the lower-dimensional fractures and fracture intersections. The deformation of explicitly represented fractures is modelled by contact mechanics relations and a Coulomb friction law, with particular attention on coupling of fracture dilation to the governing equations in both fractures and matrix. The model is discretised using multi-point finite volumes for the balance equations and a semismooth Newton scheme for the contact conditions and is implemented in the open source fracture simulation toolbox PorePy. Finally, simulation studies demonstrate the model's convergence, investigate process-structure coupling effects, explore different fracture dilation models and show an application of the model to a 3d geothermal pressure stimulation and long-term cooling scenario.


Introduction
Fluid injection operations into the subsurface are common in e.g. geothermal energy and petroleum production, wastewater disposal, CO 2 storage and groundwater management. Injection can severely alter subsurface hydraulic, mechanical, thermal and chemical conditions. These coupled processes are strongly affected by preexisting fractures, which represent extreme heterogeneities and discontinuities in the formation. The processes may in turn cause deformation of the fractures, giving rise to dynamic and highly complex process-structure interactions.
In some subsurface engineering operations, fracture deformation is deliberately induced, e.g. to enhance permeability through hydraulic stimulation, in which fluid is injected at elevated pressure to overcome a fracture's frictional resistance to slip [1,2,3]. There may also be interest in preventing deformation of fractures to, for example, avoid induced seismicity of unacceptable magnitude in disposal of wastewater [4,5,6,7] or during hydraulic stimulation of fractured geothermal reservoirs [8,9,10].
As data related to subsurface dynamics are limited, physics-based modelling can complement data analysis in understanding governing mechanisms for fracture deformation. This requires numerical simulation tools that can capture the governing structure of the fractured formation and relevant coupled processes as well as process-structure interactions, which necessitates explicit representation of both the matrix and dominant fractures in the model. Typically, major fractures or faults are represented explicitly while the rest of the domain is represented as a matrix continuum, possibly integrating effects of finer-scale fractures.
In a spatial grid, there are two alternatives for representing such a discrete-fracture-matrix (DFM) model: Resolving the width of the fractures in the grid in an equi-dimensional model imposes severe restrictions put on the spatial discretisation of the domain due to the high aspect ratio of the fractures, thereby limiting the number of fractures that can be included in the model. A geometrically simpler alternative, which was introduced for flow models, is a codimension one model, where fractures are represented as objects of one dimension lower than the surrounding domain [11,12,13,14]. In contrast to simulation models for coupled flow and mechanics that treat faults as equidimensional zones of different rheology resolved in the grid [15,16,17], the co-dimension DFM model facilitates modelling of fracture slip and dilation [18,19] and can be combined with full mechanical fracture opening [20,21]. A conceptually simpler alternative to co-dimension one DFM models is to incorporate only the dynamics in the fracture network and either disregard the dynamics in the matrix altogether or approximate them using semi-analytical methods. These approaches are based on Discrete-Fracture-Network (DFN) representations [22,23] and will be referred to as DFN methods.
Driven by the need to improve the result of injection operations and avoid unacceptable environmental impacts, intense focus has been placed on physics-based modelling. Early works by Willis-Richards et al. [24], Rahman et al. [25], Kohl and Mégel [26] and Bruel [27] developed DFNtype models considering only deformation and flow in the fractures and using a Coulomb friction law to model fracture slip due to changes in effective stress as a consequence of local change in fluid pressure. Later, Baisch et al. [28] improved on this type of model by including redistribution of shear stress along the fracture as a consequence of slip through a block-spring model. McClure and Horne [29] further developed the modelling of mechanical interaction between fractures with the boundary integral equation method and introduced a rate-and-state friction model. This type of method has been combined with fracture propagation [30,31]. As only the fracture is discretised when using the boundary integral equation method, models based on this approach can be classified as DFN-type models. Common to all of these approaches is use of semi-analytical approaches and sequential coupling of physical processes.
The last decade has seen developments in the inclusion of dynamics in the matrix as well as improved models and numerical solution schemes for coupling of different dynamics. Building on previously developed DFN-type models, McClure and Horne [29] and McClure [32] introduced a semi-analytical leakoff term to mimic fracture-matrix flow. Norbeck et al. [30] expanded on previous models developed by McClure and Horne [29] and accounted for the interaction between fracture and matrix flow through an embedded discrete fracture model, where flow in the fracture and in the matrix are discretised on non-conforming grids and connected through transfer terms. Hydro-mechanical simulation tools based on co-dimension one DFM models combined with Coulomb friction laws for fracture slip have also been introduced [19,33,34,35,36,37,38], motivated by applications related to CO 2 -storage [35], oil and gas production [36,37] and hydraulic stimulation of fractured geothermal reservoirs [38,19]. In these tools, the contact mechanics conditions are typically handled through Lagrange-multipliers [35,34,33] or penalty methods [36]; see e.g. Wriggers and Zavarise [39].
More recently, thermal effects have been taken into account in deformation of fractured porous media. Based on a DFN-type model, where the boundary integral equation method was used so that only the fracture is discretised, Ghassemi and Zhou [40] included thermo-poroelastic effects in the matrix. Based on a DFM conceptual model, Pandey et al. [41] and Salimzadeh et al. [42] have presented models with linear thermo-poroelasticity for the matrix combined with flow, heat transfer and deformation of a single fracture. However, none of these works included modelling of fracture slip or shear dilation when fracture surfaces are in contact. Gallyamov et al. [21] consider a conceptually similar model which includes multiphase flow and a fracture-contactmechanics model combined with opening and propagation of fractures, and present simulation studies with a large number of fractures. Their work considers the impact of the contact traction on the hydraulic aperture of closed fractures. In contrast to the majority of previously mentioned works, where simplifications that impact the solution are made in the solution of the coupled system, they solve the equations fully coupled, i.e. the flow, energy and mechanics equations are solved simultaneously building on the work by Garipov et al. [36].
Recent work by Garipov and Hui [43] combines several previous developments. Their work is based on a DFM model and considers a fully coupled thermo-poroelastic model for the matrix, flow and heat transfer in the fractures and contact mechanics for fractures based on a Coulomb friction law. Energy and mass conservation are discretised by a finite volume method, while momentum is discretised by a Galerkin finite element method. This work also presents robust treatment of couplings in the model. However, while the work accounts for permeability enhancement due to full opening of fractures as well as shear dilation, stress response due to dilation as a consequence of slip is not included in the model. This paper presents a mathematical model based on a mixed-dimensional DFM representation of coupled thermo-hydro-mechanical (THM) processes in a porous rock containing deforming fractures with an accompanying discretisation and numerical solution approach. The model fully couples fluid flow and transport in both matrix and fractures, linear thermo-poromechanics in the matrix and nonlinear fracture deformation. Fracture deformation is based on traction balance, nonpenetration and a Coulomb type friction law, and allows for shear slip and dilation as well as complete fracture opening. To the authors' knowledge, this is the first model that consistently and fully coupled represents stress redistribution due to slip-induced dilation of fractures. As demonstrated by the numerical results, the effect of this coupling can be significant.
Based on the modelling of fractures as lower-dimensional surfaces, the domain is decomposed into subdomains of different dimensions corresponding to matrix, fractures and intersections. Model equations, sets of variables and parameters are defined on each subdomain and the interfaces between them. The resulting mixed-dimensional model [44] facilitates systematic modelling on the decomposed structure while incorporating interaction between processes both within and between subdomains. The governing balance equations in each subdomain are discretised based on multi-point finite volume methods preserving local conservation, using the same spatial grid for discretisation of all processes. The contact conditions are formulated in terms of contact tractions at the interfaces between fractures and matrix blocks, similar to approaches based on Lagrange-multipliers [39]. The nonlinear fracture deformation equations are discretised using a semismooth Newton scheme formulated as an active set method. The model and its implementation extend the work presented by Berge et al. [34], who consider the poromechanical problem without flow in the fractures, shear dilation and thermal effects.
The model is presented in Section 2, and its discretisation is described in Section 3. In both sections, particular emphasis is placed on fracture deformation as well as its impact on the balance equations for the fractures and the back-coupling to the higher-dimensional momentum balance. Three examples are presented in Section 4: The first investigates governing mechanisms and coupling effects and verifies the model and its implementation in a convergence study. In the second, three different models for fracture dilation are compared. In the last example, the model is applied to a three-dimensional hydraulic stimulation and long-term cooling scenario for geothermal energy extraction. Finally, Section 5 provides some concluding remarks.

Mixed-dimensional governing equations
This section describes the model for THM processes in a porous medium with contact mechanics at the fractures. It relies on a DFM model in which the matrix, the fractures and fracture intersections are explicitly represented by individual subdomains. To avoid resolving the small geometric distances introduced to the fracture network geometry by the high aspect ratio of the fractures, the dimensions of the fracture and intersection subdomains are reduced. The subdomains are collected in a hierarchical structure and connected by interfaces to yield the full mixed-dimensional model.
Decomposition into subdomains facilitates tailored modelling of processes in distinct subdomains, while interactions between subdomains take place on the interfaces. Specifically, separate sets of variables, equations and parameters are defined on each subdomain and interface. This procures the flexibility needed to model the highly complex system arising from the coupled THM system posed in both matrix and fractures.
The model consists of balance equations for momentum, mass and energy and relations governing the fracture deformation posed on the subdomains. These are supplemented by constitutive laws and equations for coupling over the interfaces. The equations are formulated in terms of the primary variables displacement, pressure, temperature and contact traction on the fractures.
The governing equations for THM-processes in a mono-dimensional porous medium are introduced succinctly in Section 2.1, followed by a more elaborate presentation of the lowerdimensional scalar equations for deforming fractures and intersections emphasising the effect of volume change in Section 2.2. Section 2.3 describes the model for fracture deformation and its relation to volume change.

Matrix thermo-poromechanics
We first consider the governing equations in the matrix domain consisting of a solid and a fluid phase. For the remainder of this subsection, let Ω denote the matrix domain, n the outer normal vector on the boundary ∂Ω and dx = ndx. Neglecting inertia, the momentum balance equation reads with q u denoting body forces and the thermo-poroelastic stress tensor for infinitesimal deformation modelled as linearly elastic obeying an extended Hooke's law Here, D denotes the stiffness tensor, α the Biot coefficient, β s the volumetric thermal expansion and K s the bulk modulus of the solid, while u, p, and I are displacement, pressure and identity matrix. The subscript 0 indicates the reference state of a variable. The assumption of local thermal equilibrium between fluid and solid leads to a single temperature unknown T and to the definition of effective coefficients, which are computed as weighted sums [45]: Subscripts s and f indicate that all quantities in the parentheses are for the solid and fluid, respectively. Herein, the relations D 2 (∇u + ∇u T ) = µ(∇u + ∇u T ) + K s tr(∇u)I and q u = ρ s g are used, with µ denoting the shear modulus, tr(·) the trace of a tensor, ρ s the solid density and g the gravitational acceleration vector.
We assume the density of the slightly compressible fluid to follow with β f and K f denoting fluid thermal expansion coefficient and bulk modulus, which is the inverse of the compressibility. Balance of mass reads [46] with porosity φ and effective thermal expansion β e . Fluid flux relative to the solid is denoted by v and volume sources and sinks by q p . With K denoting the permeability and η the viscosity, the flux is modelled according to Darcy's law: Neglecting viscous dissipation on an assumption of small velocities [46], the energy balance equation is [47,48,46] Here, the internal energy of the solid is U s = ρ s c s T, where c s is the specific heat capacity of the solid. Assuming a simplified low-enthalpy description of the fluid, we approximate the fluid internal energy as In the computation of the effective heat conductivity κ e by Eq. (3), dispersion due to microscale tortuous flow in the porous medium is neglected. The source term is assumed to equal the internal energy of the fluid of the volume source and sink terms, i.e. q T = ρ f c f Tq p . The energy equation can thus be written where (ρc) e denotes the effective heat capacity of the porous medium and is calculated by Eq. (3). The first term of Eq. (10) can be expanded, giving where effective parameters are calculated using Eq. (3). Figure 1: Schematic representation of a throughgoing one-dimensional fracture in a two-dimensional matrix (left) and three fractures meeting at an intersection point (right). All subdomains, internal boundaries and interfaces are indicated, as are select projection operators. Ω h and Ω l are separated by the interfaces Γ j and Γ k corresponding to the internal boundaries ∂ j Ω h and ∂ k Ω h . The intersecting fractures are indexed counterclockwise from 1 through 3 and intersect at the point Ω 4 . The interfaces are numbered so that Γ i matches ∂ i Ω i . In the model, Ω l , Γ j , Γ k , ∂ j Ω h and ∂ k Ω h coincide geometrically, as do all zero-dimensional points in the right figure.

Lower-dimensional flow and heat transfer
This section derives balance equations for mass and energy for fluid-filled fractures and intersections which may undergo significant relative deformation and volume change, giving rise to an additional term compared to equations for static domains. Along with outlining dimension reduction for the mass and energy equations, the connection between the subdomains of the mixed-dimensional model is presented. Some notation is needed for a unified description of the mixed-dimensional model of a fractured porous domain of dimension D = 3 or D = 2. The domain is split into connected subdomains corresponding to the rock matrix, the co-dimension one fracture planes and co-dimension two fracture intersections. In the case D = 3, the model also generalises to account for intersections of fracture intersection lines, i.e. zero-dimensional points. A subdomain is denoted by Ω i and its boundary by ∂Ω i . The subscript i is also used to identify variables defined within Ω i , but suppressed as context allows. Each part ∂ j Ω i of the internal boundary is associated with an interface Γ j to an immersed lower-dimensional domain Ω l (see Fig. 1). All lower-and higher-dimensional interfaces of a subdomain are collected in the setsŠ andŜ; in particular, the interfaces corresponding to surfaces of fracture i constituteŜ i . Where convenient, the higherand lower-dimensional neighbours of an interface are denoted by Ω h and Ω l , respectively.
Finally, four types of projection operators are needed to transfer variables between interfaces and the neighbouring higher-and lower-dimensional subdomains. As illustrated in Fig. 1, projection from the interface to the subdomains is performed by Ξ h j and Ξ l j , respectively, whereas Π h j and Π l j project from the part of a subdomain geometrically coinciding with the interface to the interface.
The thickness of a fracture is characterised by the aperture a [m], which will be related to the fracture deformation in Section 2.3. The aperture of an intersection is taken to be the average of the intersecting higher-dimensional neighbours, i.e.
where the projection operators transfer the higher-dimensional aperture first to the interface and then to the intersection, see Fig. 1. The specific volume V i = a D−d accounts for the dimension reduction from the deforming equi-dimensional Ω to the corresponding spatially fixed d-dimensional Ω i . With V = 1 for d = D, dimension reduction for a scalar quantity ζ and a vector quantity ι satisfy and respectively. Here, ι i denotes the tangential d-dimensional flux, while the interface flux into the domain ι j satisfies Ξ h j ι j = ι h · n h on ∂ j Ω h , with n h denoting the outwards normal vector of Ω h . The weighting with V h ensures that the interface flux matches the dimension of fluxes of the higher-dimensional neighbour, which are scaled by specific volumes as seen by the expression for the tangential flux. Note that all differentials in reduced integrals should be interpreted as relative to the domain of integration; i.e. dx is two-dimensional for a fracture with d = 2. Furthermore, dx = n i dx, where n i denotes the outwards normal at ∂Ω i lying in the tangent plane of Ω i . For d = 1, the boundary integral equals evaluation of the integrand at the boundary points.
Assuming unitary fracture porosity, the fluid mass balance equation for a deforming equidimensional domain Ω is Averaging in the normal direction for the tangential flux and replacing the normal part of the boundary by Γ j according to Eq. (14), the fluid flux becomes Thus, the interdimensional coupling between Ω h and Ω l takes the form of interface fluid fluxes v j , which also appear as a Neumann condition for Ω h : The interface flux is modelled using a Darcy type law extended from Martin et al. [12] to account for gravity Both the weighting by a l and V h and the normal permeability K j arise through dimension reduction. The remaining terms of Eq. (15) are averaged in the normal direction using Eq. (13), and Eq. (4) is inserted for the fluid density. Collecting terms, dividing by ρ f and assuming Darcy's law with tangential permeability K yields the dimensionally reduced mass balance with the third term accounting for changes in fracture volume. The dimension reduction is now performed for the energy balance, which for an equi-dimensional domain reads where w and q are given by Eqs. (8) and (9). Utilising Eq. (14), the dimension reduction of the flux terms is and the internal boundary conditions on ∂ j Ω h are The advective interface flux is defined according to the upstream direction of the interface fluid flux: The Fourier-type conductive interface flux is with the normal heat conductivity modelled as κ j = Π l j κ f,l since it originates from the dimension reduction of a fluid-filled domain.
The dimension reduction of the remaining terms of Eq. (20) is performed by again invoking Eq. (13), yielding where the second term in the first integral is decomposed similarly to Eq. (11), with φ = 1 in the calculation of effective coefficients.

Fracture contact mechanics
The traction balance, nonpenetration condition and friction law posed on a fracture l are formulated in terms of interface displacements and fracture contact traction. The interface displacements on the two surfaces Γ j and Γ k are u j and u k , and the jump between the two sides is Since the fracture deformation depends on traction caused by the contact between the two surfaces, the contribution from p l should be subtracted on the fracture surfaces to yield the traction balance posed on the interfaces: The fracture contact traction λ l will for notational convenience be referred to as λ in the following, and is defined according to the normal of the fracture, which is defined to equal n h on the j side, i.e. n l = Ξ l j Π h j n h . When there is no mechanical contact between the interfaces, λ is 0, implying that the higher-dimensional thermo-poromechanical tractions projected to the interfaces on the right-hand sides of Eq. (27) are balanced by the fracture pressure.
A vector ι defined on a fracture may be decomposed into the normal and tangential components ι n = ι · n l and ι τ = ι − ι n n l .
With this notation, the nonpenetration condition reads with the gap function g defined to equal the distance between the two fracture interfaces when in contact. The Coulomb friction law is with F denoting the friction coefficient and [[u ]] τ denoting the tangential displacement increment. In addition to enforcing the traction balance of Eq. (27) and the conditions of Eqs. (29) and (30), a Dirichlet condition is assigned on ∂ j Ω h so that The aperture introduced in Section 2.2 is a function of displacement jump, a = a ([[u]]). Due to roughness of the fracture surfaces, tangential displacements may induce dilation [49] as illustrated in Fig. 2. The relationship between the dilation and the magnitude of tangential displacement is assumed to be linear and described by the dilation angle ψ following Rahman et al. [25]. As modelled herein, the dilation is not merely a hydraulic effect impacting e.g. the fracture permeability, but a mechanical effect in the sense that the normal distance between the fracture surfaces increases. As such, the dilation must be coupled back to the normal interface displacements and the matrix deformation through Eq. (31), which is achieved by choosing the gap function The update is reversible; if the tangential displacement is reversed, g takes on its initial value. Small-scale fracture roughness may provide a volume for the fluid to occupy even when the fractures are in an undeformed state. This leads to the following relation between aperture and displacement: where a 0 denotes the residual aperture in the undeformed state.
In addition to entering the equations as a result of dimension reduction, a governs the tangential permeability of a fracture or intersection line i according to the cubic law [50], where I i denotes the identity matrix of the fracture dimension. Equation (34) constitutes a strongly nonlinear coupling, especially as K i is multiplied by V in Eq. (19). Finally, the normal permeability of an interface is inherited from the lower-dimensional neighbour:

Discretisation
This section describes the discretisation of the model presented in the previous section. The system is discretised in time using implicit Euler and solved monolithically using a direct solver [51]; a more scalable option would be to use iterative methods with block-preconditioners, following ideas in e.g. [52,53]. The spatial grids are simplicial, and are constructed such that the lower-dimensional cells coincide with higher-dimensional faces; grids are generated by Gmsh [54]. The model is implemented in the open-source fracture simulation toolbox PorePy presented by Keilegavlen et al. [44]. The mixed-dimensional framework gives rise to a two-level block structure as the equations are discretised. The outer level corresponds to the subdomains and interfaces, with entries internal to the subdomains on the diagonal and entries for the interdimensional coupling on the off-diagonals. The inner level corresponds to the primary variables, with coupling effects between different variables on the off-diagonals. The block structure is illustrated in Fig. 3, which will be used in the following description of discretisation of individual terms by referring to the block in row r and column c as A (r,c) with r and c ranging from 1 to 14. The two-level block structure represents the left-hand side of the linear equation system for a matrix Ω h , a fracture Ω l and interfaces Γ j and Γ k with corresponding equation numbers shown to the left. Bottom left: Spatial discretisation and spatial location of degrees of freedom for a domain corresponding to the equation system. Matrix, fracture and interface grids are shown in black, green and orange, respectively, and the corresponding degrees of freedom are shown as squares, triangles and diamonds. Black represents displacement, blue pressure and red temperature; the relation between all markers and unknowns is shown at top right. Bottom right: Subgrid around a node xn of the primary grid, which is shown in solid black lines. The O-shaped interaction region forming the stencil for the local systems is constructed by connecting the surrounding cell centres xν and face centres x f as indicated by the dotted lines. Continuity of primary variables is enforced in the points xc; the shaded area indicates a subcell.

Matrix thermo-poromechanics
The spatial discretisation of the diffusive terms of the balance equations is achieved using a family of cell-centred finite volume schemes. The approach is based on the multi-point flux approximation (MPFA) [55] defined for diffusive scalar problems and the multi-point stress approximation (MPSA) for vector problems [56] and their combination for THM problems [57,58]. The scheme is formulated in terms of discrete displacement (D vectors), pressure and temperature unknowns and is locally momentum, mass and energy conservative.
The scheme's construction is based on a subdivision of the spatial grid as illustrated in Fig. 3, with the gradients of displacement, pressure and temperature defined as piecewise constants on the subdivision. The fluxes of the conserved quantities momentum, mass and energy are discretised via Hooke's, Darcy's and Fourier's law, respectively. Continuity is enforced for traction and mass and energy fluxes over faces of the subgrid and for the primary variables in the continuity points x c , leading to one local system for the node of the primary grid. Each local system is partially inverted to express gradients in terms of the cell-centre values in nearby cells. A global system is constructed by collecting for each cell all face fluxes as expressed in terms of the cell-centred primary variables. For details, see Nordbotten and Keilegavlen [58].
The coupling between the three equations is achieved by using the thermo-poroelastic stress for the local traction balances, which directly yields the contributions A (1,2) and A (1,3) representing the scalar variables' effect on the momentum balance. A (2,1) and A (3,1) , which represent the displacement effects on the scalar balances, are constructed by assembly of the discrete divergence based on the local systems for the displacement gradients.
The standard finite volume implicit Euler discretisation is applied to all time derivatives; that is, both the TH coupling blocks A (2,3) and A (3,2) and the accumulation terms of A (2,2) and A (3,3) . The advective term of Eq. (10) is discretised using a first-order upwind scheme, i.e. the temperature flux between cells k and l is with the fluid flux from cell k to cell l v k,l and ρ f computed from the solution at the previous iteration.

Mixed-dimensional flow and heat transfer
All terms of the scalar equations for the lower-dimensional subdomains are discretised using lower-dimensional versions of the corresponding D-dimensional discretisations. For the Darcy and Fourier fluxes, this implies that we use the MPFA scheme, while the advective fluxes are again treated by first-order upwinding. The interdimensional coupling relations are discrete analogues to Eqs. (18), (24) and (23). Thus, they involve reconstruction of p and T on ∂ j Ω h , which we base on discretisation matrices pertaining to the MPFA discretisations. For the matching grids used herein, the discrete projections are straightforward bijective mappings between faces of Ω h and cells of Γ j (Π h j and Ξ h j ) and between the cells of Ω l and Γ j (Π l j and Ξ l j ). The nonlinearities arising through the products involving a and V are solved iteratively within the Newton scheme for fracture deformation described below. Specifically, the time derivatives are computed as additional right hand side terms based on values from the previous iterate and time step. However, the linear volume-change terms in the fractures are coupled fully implicitly to u j so that the contribution for each fracture is the jump between the neighbouring higherdimensional interfaces, as illustrated by the off-diagonal blocks A (5,7) , A (5,11) , A (6,7) and A (6,11) . Densities are computed from the solution at the previous iteration. In some simulations involving strong advection and high temperature gradients, the density dependence in the gravity term of Darcy's law may lead to oscillatory fluxes between Newton iterations. This may result in convergence problems related to the upstream discretisation of the advective term. In these situations, convergence was achieved by damping the updates of the fluid flux of the advective term.

Fracture contact mechanics
Fracture deformation discretisation is based on the approach presented by Hüeber et al. [59] and Wohlmuth [60] with the frictional contact problem formulated as a variational inequality. The formulation is expanded to account for the [[u]] τ dependency of g. Deformation constraints are reformulated as complementary functions C = C(X), with X being the unknowns. The constraints are imposed by solving C = 0 through application of the semismooth Newton method where δX k = X k+1 − X k . Correspondingly, the increment of a function f between successive iterations k and k + 1 is δf . D is the generalised Jacobian of C, i.e. the convex hull of the standard Jacobian wherever C is differentiable.
To facilitate imposition of different constraints based on the deformation states defined in Eqs. (29) and (30), three disjoint sets describing the deformation state as open, sticking or gliding are defined: Here,c denotes a numerical parameter, the friction bound is b = −F (λ n +c ([[u]] n − g)) and [[u ]] τ denotes the increment from the previous time step. Replacing [[u ]] τ by [[u]] τ in the above definition yields the cumulative fracture state sets, which are denoted by subscript c.
The normal and tangential complementary functions are and The corresponding generalised Jacobians are and 13 Here, χ is the characteristic function of a set for a fracture cell ν, while the increment of the friction bound is with g u denoting the derivative of g with respect to [[u]]. Hence, sorting each cell according to Eq. (38) and imposing Eq. (37) results in the following constraints: The coefficients L, ν and r are functions of [[u]] k τ and λ k , and can thus be computed from the previous iterate. For the exact expressions and further details of the discretisation and implementation of the fracture deformation equations, see Berge et al. [34].
The effect of letting g depend on [[u]] τ only appears in the normal condition in the two terms involving the derivative g u . The two cases g = 0 and Eq. (32) will be considered below. The former obviously gives g u = 0 while the latter gives which may be inserted into Eq. (45) to finally yield A (4,4) , A (4,7) and A (4,11) of Fig. 3.

Numerical results
This section presents three sets of simulations aimed at demonstrating the model's representation of complex process-structure interactions. In the first example, a convergence study is presented and coupling mechanisms investigated. The second example explores different modelling choices for the relationship between displacement jumps and apertures. Finally, the model is applied to a geothermal scenario with a pressure stimulation phase and long-term cooling during a production phase. Run scripts for the example simulations and animations showing temporal evolution of the solutions may be found in a dedicated GitHub repository [61].

Example 1 -Convergence study
While different components of the implementation have been verified in previous studies [44,34], analytical solutions probing the full model presented herein are not available. The first example is designed as a validation of the implementation: Starting from a coarse grid of 398 2d cells, 38 1d cells and two 0d cells, a sequence of six grids is produced by nested conforming refinement. The finest grid, which has 407 551 2d cells and a total of 1 647 254 unknowns, is used as the reference solution for the convergence study and forms the basis of the discussion of coupling mechanisms. The geometry of the two-dimensional domain with eight fractures is a modified version of a geometry presented in Berge et al. [34] and is shown in Fig. 4. It contains a kink formed by two fractures, an intersection formed by two other fractures and nearly intersecting fractures, as well as both immersed fractures and one fracture extending to the boundary. These features can be expected to challenge the accuracy of numerical simulations.
Simulating three different phases allows us to distinguish between the influence of mechanical, hydraulic and thermal driving forces. The three phases are defined through the boundary conditions as follows: Fixing the bottom and setting homogeneous stress conditions on the left and right boundary, a Dirichlet displacement value of (5 × 10 −4 , −2 × 10 −4 ) T m is applied at the top throughout the simulation and is the only driving force during phase I. Phase II begins when a pressure gradient of 4 × 10 7 Pa is applied from left to right. Once the solution has reached equilibrium, a boundary temperature 15 K lower than the initial temperature is prescribed at the left boundary marking the onset of phase III. The initial values are p = 0 Pa, T = T 0 = 300 K and a 0 = 5 × 10 −4 ; no gravity effects are included in this example. Figure 5 shows convergence results for the end of the three phases. For each phase, we plot the errors for three primary variables on individual subdomains for the different refinement levels. The variables are displacement jumps, contact tractions and the variable related to the main driving force of the phase. The error is computed by projecting the cell-centre value of the coarse grids onto the reference grid and then computing the L 2 norm of the difference between coarse and fine solution. Errors are normalised by the number of reference cells in the subdomain multiplied by a weight k representing the magnitude of the global range of the variable in question. The weights are obtained from the boundary conditions and are k u = ||(5 × 10 −4 , −2 × 10 −4 )|| m, k p = 4 × 10 7 Pa, k T =15 K and k λ = Ek u with E denoting Young's modulus.
In general, the expected first order convergence is observed. The exception is traction on some of the fractures (1, 2, 3 and 6). These local errors may be attributed to the geometrical challenges posed by those fractures: Fractures 1 and 2 meet in a kink, which seems to lead to relatively large errors compared to the remaining fractures as discussed by Berge et al. [34]. Fracture 3 intersects fracture 4, while the error on fracture 6 is concentrated around the leftmost tip, which is close to the neighbouring fracture 5. However, the traction solutions converge for all fractures without small transition regions and the challenging geometrical features have no discernible effect on the convergence in the other primary variables. Therefore, taken together,   the presented results serve as a verification of the model. Figure 6 shows the fracture deformation for each of the three phases, and thus demonstrates the effect of each of the three driving forces. The richness in physical processes and the complexity of coupling in the fractured THM problem is well illustrated by a phenomenon observed towards the end of phase III around fractures 5 and 6 (see Fig. 7). The role of fractures as preferential flow pathways leads to high flow rates and cooling in the region where fluid leaves fracture 5, both at the tip closest to the right boundary and in the area closest to fracture 6. This, in turn, leads to local contraction of the matrix and fracture dilation -in this particular case both through shear displacement and normal opening as seen from the final deformation state (bottom left in Fig. 6). The dilation further increasing the fracture conductivity can be expected to enhance the effect, which is also observed at the tip of fracture 4 somewhat earlier in the simulation. This phenomenon of enhanced cooling-induced aperture increase in regions where the fluid enters or leaves a fracture can be expected to be of a general character.

Example 2 -Fracture dilation models
The second example is a study of different models for fracture dilation based on simulation of the case described in Section 4.1 with two simplified aperture models. In the first simplification, model 0, there is no coupling between shear displacement and dilation, i.e. g = 0 and a = a 0 +[[u]] n . In the second simplified model (1)  A comparison in terms of the final spatial distribution of aperture increase and tangential displacement jump on each of the closed fractures is shown in Fig. 8. As the dilation relations are irrelevant for open fractures, analysis is based on the mostly closed fractures 5 through 7. Fractures 6 and 7 clearly demonstrate how the dilation coupling in model 2 reduces tangential displacement compared to the simplified methods, as the induced normal displacement increases the normal traction on the fractures. Interestingly, the apertures displayed in Fig. 8 show overestimation for the one-way coupling due to the above-mentioned overestimation of the tangential jumps. Model 0 obviously yields no shear dilation. The model 0 aperture increase of fracture 5 is thus related to the fracture being open. Note that part of this region is closed for model 2 (cf. the bottom left illustration of Fig. 6) demonstrating how inconsistency affects the results beyond the prediction of a.
The magnitude of the discrepancy must be expected to depend on rock and fault parameters, particularly the dilation angle. The results demonstrate the strong coupling in the problem, and that compromising this coupling in numerical modelling may lead to significant error in the results.

Example 3 -Hydraulic stimulation and long-term cooling of a geothermal reservoir
The third example shows hydraulic stimulation of a geothermal reservoir, followed by an injection and production phase leading to long-term reservoir cooling for the three-dimensional geometry in Fig. 9. .81 m s −2 denoting the gravitational acceleration. After letting the system reach equilibrium under the mechanical boundary conditions representing an anisotropic background stress in phase I, we simulate a pressure stimulation phase (II) and a production and long-term cooling phase (III). In the 10 hour stimulation phase, the flow rates of the injection and production wells are 75 and 0 L s −1 , respectively. During the 15 year production phase, both rates are 20 L s −1 . The injection temperature is 70 K below the reservoir temperature. The wells are incorporated as source terms in the fracture cells intersected by the well paths, with upwind discretisation for the energy equation in the production cell. Hydrostatic Dirichlet boundary conditions p = p H apply for the pressure. An anisotropic compressive background stress is imposed with the following non-zero stress tensor values While the remaining parameters listed in the Table 1 are plausible for geothermal reservoirs, they do not correspond to a specific site. The number of Newton iterations required for convergence is below ten for all time steps. The results are summarised through the temporal evolution of the norm of the displacement jumps on the three fractures shown in Fig. 9. Significant stimulation effects appear in both phase II and phase III, with the magnitude of the jumps somewhat larger during the cooling; dynamics initiate latest on the injection fracture due to its orientation relative to the background stress. Figure 10 shows spatial plots of pressure, temperature, aperture, displacement jumps and deformation state. The plots demonstrate the model's cell-wise spatial resolution of the dynamics both in fractures and matrix. For all three phases, displacement jumps are oriented in agreement with the background stress field and are very closely aligned. The only cells in O are around the intersection in phase III. For the remaining cells, the (relatively small) normal components of the orientation arrows are due solely to shear dilation.
During phase II, aperture increments are most pronounced on fracture 2, which has no wells. However, its intersection with fracture 1 where injection occurs leads to a significant pressure increase. Despite negligible pressure perturbation in fracture 3, some slip is observed due to stress redistribution following the deformation of fracture 2. The location of the slip in fracture 3, away from the stress shadow of fracture 2, highlights the complex mechanical interplay between fractures in a network.
During phase III, some displacement jumps are induced in fracture 1 close to the intersection, whereas there is significant aperture increase throughout fracture 2 as a result of cooling of the surrounding rock. Along fracture 3, the deforming region is different from the previous phase, with displacement occurring in the region closest to fracture 2, where the surrounding matrix has been cooled the most. This conforms with the observations in Section 4.1 of aperture increases in regions of fluid entry or departure from the fractures.

Conclusion
A model for fully coupled thermo-hydro-mechanical processes in porous media with deforming fractures is presented. Using the discrete-fracture-matrix approach, the matrix, the fractures and the fracture intersections are represented by subdomains of different dimensions connected by interfaces in a mixed-dimensional model. Balance equations for energy and mass in all subdomains are coupled by fluxes on the interfaces, while the momentum balance in the matrix and traction balance and non-penetration for the fracture surfaces are coupled through interface displacements. These governing equations are supplemented with constitutive laws, including a Coulomb type friction law and a linear shear dilation relation for the fractures. For the latter, a novel model consistently coupling slip and shear dilation of the fractures with the stress response of the matrix is presented. The resulting set of model equations is discretised using cell-centred multi-point finite volume schemes and a semismooth Newton method for fracture deformation and solved fully coupled.
The model and its implementation are verified through a convergence study displaying firstorder convergence for all primary variables and subdomains, except for the expected local reduction of convergence in the transition between contact regimes. An exploration of three different shear-dilation models reveals significant discrepancies, demonstrating the importance of accurate and consistent modelling of the underlying physical mechanisms and their couplings.
Investigations of 2d and 3d examples identify a mechanism by which cooling-induced shear dilation preferentially occurs in regions where fluid leaves or enters a fracture. The investigations also show the complexity of the process-structure interactions which may arise: in particular, how fracture deformation and resulting fracture dilation is induced by both mechanical, hydraulic and thermal driving forces. This confirms the need for models which explicitly incorporate all relevant processes and structural features as well as the resulting process-structure interactions. Furthermore, it demonstrates the proposed model's prowess in capturing such highly complex interactions and identifying their governing mechanisms. Extensions such as chemical processes and more advanced friction models and dilation relations could readily be accommodated in the applied mixed-dimensional framework.

Acknowledgements
The authors thank two anonymous reviewers for comments and suggestions which contributed to improving the paper.
Funding: This work was supported by the Research Council of Norway and Equinor ASA through grant number 267908. I. Stefansson and I. Berre also acknowledge funding from the VISTA programme, The Norwegian Academy of Science and Letters.