A thermo-flow-mechanics-fracture model coupling a phase-field interface approach and thermo-fluid-structure interaction
Sanghyun Lee Henry von Wahl Thomas Wick





1]Sanghyun Lee
slee17@fsu.edu
2]Henry von Wahl
henry.von.wahl@uni-jena.de
3]Thomas Wick
thomas.wick@ifam.uni-hannover.de

[1]organization=Department of Mathematics, Florida State University,
            addressline=1017 Academic Way,
            city=Tallahassee,
            postcode=32306-4510,
            state=FL,
            country=USA

[2]organization=Friedrich-Schiller-Universität, Fakultät für Mathematik und Informatik,
            addressline=Ernst-Abber-Platz 2,
            city=Jena,
            postcode=07743,
            country=Germany

[3]organization=Leibniz Universität Hannover, Institut für Angewandte Mathematik,
            addressline=Welfengarten 1,
            city=Hannover,
            postcode=30167,
            country=Germany


§ ABSTRACT 
Geothermal energy, a promising renewable source, relies on efficiently utilizing geothermal reservoirs, especially in Enhanced Geothermal Systems (EGS), where fractures in hot rock formations enhance permeability. Understanding fracture behavior, influenced by temperature changes, is crucial for optimizing energy extraction. To address this, we propose a novel high-accuracy phase-field interface model integrating temperature dynamics into a comprehensive hydraulic-mechanical approach, aiming for a thermo-fluid-structure interaction representation. Therein, the key technical development is a four-step algorithm. This consists of computing the fracture width, reconstructing the sharp interface geometry, solving the thermo-fluid-structure interaction (TFSI) problem, and employing a phase-field approach coupled to the temperature and pressure from the TFSI problem. By coupling temperature-hydraulic-mechanical processes with our newly proposed high-accuracy phase-field interface approach, we investigate how temperature impacts fracture width values, which are crucial for permeability in EGS reservoirs. Through this model and three different numerical simulations, we aim to provide an approach to deepen understanding of the complex interplay between temperature, mechanical deformation, and permeability evolution. Therein, we substantiate our formulations and algorithms through mesh convergence results of crack width and total crack volumes for static fractures, and crack lengths in the case of propagating fractures.




phase-field fracture thermo-hydro-mechanics thermo-fluid-structure interaction smeared interface sharp interface







§ INTRODUCTION



Geothermal energy is a promising renewable energy source, offering sustainable power generation with minimal environmental impact. Understanding the behavior of geothermal reservoirs is crucial for efficient and sustainable exploitation of this resource. In addition, Enhanced Geothermal Systems (EGS) represent a frontier in geothermal energy development, where the creation and stimulation of fractures within hot rock formations are essential for enhancing permeability and facilitating fluid circulation <cit.>. In such systems, understanding the behavior of fractures and their response to changes in temperature is essential to optimize reservoir performance and energy extraction efficiency <cit.>.

The crack opening displacement or fracture width values within these fractured reservoirs serve as critical permeability indicators, directly influencing fluid flow rates and heat transfer capabilities <cit.>. However, accurately predicting these fracture width values requires comprehensive models that account for the coupled effects of temperature, hydraulic processes, mechanical deformation, and fracture propagation.

In numerous works over the last two decades, the phase-field fracture method has been shown to be one of the most effective approaches due to the existence of the diffusive zone <cit.>. Furthermore, monographs and extended papers of phase-field methods for fracture propagation with multi-physics applications include <cit.>. These methods are in contrast to sharp interface approaches such as XFEM <cit.>. However, coupling different physical phenomena may be of interest near or across the interface. The challenge with the phase-field approach lies in the accurate modeling of fundamental physical principles, due to the inherently diffusive nature of the fracture zone, which complicates precise crack boundary localization for modeling interface-related physics. As a result, the phase-field method, while powerful for propagation and representing fracture patterns in two and three dimensions, requires careful consideration when applied to problems where complex physics at the interface between the fluid-filled fracture and the surrounding solid are critical to the fracture width <cit.>. This is particularly the case when interactions between the fluid and the surrounding solid and temperature effects need to be included. This results in a classical fluid-structure interaction problem with thermal effects in our case, where a sharp interface is required. As a result, the phase-field method, while powerful, requires careful consideration when applied to problems where the exact fracture width plays a critical role.

To address this issue, we build upon a split approach introduced in <cit.>. This approach combines a diffusive interface phase-field model for fracture dynamics and an interface resolving fluid-structure interaction problem. To combine these two somewhat opposing approaches, we consider a geometry reconstruction approach <cit.>. To reconstruct the geometry of the open fluid-filled fracture, we use the crack opening displacements or the fracture width to describe the interface between the fluid and the solid. This gives a flexible method for switching between the interface-capturing phase-field method and the interface-tracking fluid-structure interaction approach.

In this paper, we extend the algorithm to consider temperature dynamics into a comprehensive hydraulic-mechanical model using a phase-field approach. By coupling the temperature-hydraulic-mechanical (THM) processes with the phase-field approach, we aim to provide a more realistic representation of fluid-structure interaction processes within geothermal reservoirs, particularly in the context of EGS. Other works considering THM related to phase-field and porous media, include <cit.>. Specifically, we model the flow-temperature part through a Boussinesq approximation <cit.>. The temperature then enters into the solid equation via the stress tensor; see, e.g., <cit.>, <cit.> and <cit.>. The THM process is modeled and simulated in the reconstructed domain, which fully resolves the fluid reservoir. However, this differs from the domain in which classical phase-field fractures are considered. The latter is a problem usually posed using a smeared zone of a very thin fracture. This contrasts with the resolved fluid-filled reservoir considered for the THM problem. Therefore, to take quantities, such as the fluid pressure and the temperature from the THM problem, and couple these to a phase-field model, we must derive a novel phase-field fracture model. Consequently, this allows us to switch back to the phase-fracture model on the reconstructed geometry and directly use information from the THM model in the phase-field model.

With this new model at hand, we present a detailed investigation into the influence of temperature on fracture width values using our coupled THM-phase-field model. We hypothesize that temperature variations play a significant role in altering the mechanical properties of the reservoir rock, thereby influencing the opening and closure of fractures and, consequently, permeability enhancement. By incorporating temperature effects into our model, we anticipate gaining deeper insights into the thermal behavior of geothermal reservoirs and its implications for production optimization and reservoir management in EGS.



The remainder of this paper is structured as follows. First, we introduce the system of equations that model the coupled non-isothermal fluid-solid system in <Ref>. Then, we propose our novel algorithm based on a coupled iteration between the thermo-fluid-structure interaction problem and a phase-field fracture approach. In <Ref>, we then present a detailed derivation of our fluid-filled fracture phase-field model, which incorporates both pressure and temperature effects at the fracture interface. We then describe the details of how we recover the geometric information from the diffusive phase-field model to reconstruct the sharp interface geometry in <Ref>. Furthermore, we present the weak formulation of the thermo-fluid-structure-interaction problem in this section. In <Ref>, we then compute several numerical examples to validate and demonstrate the capabilities of the proposed algorithm. Finally, we give some concluding remarks in <Ref>.




§ GOVERNING EQUATIONS AND OVERALL CONCEPT





 §.§ Modeling Overview



Consider the computational domain Ω⊂ℝ^d (d∈{2,3}). 
We assume that this is subdivided into the fluid domain Ω_f and the solid domain Ω_s, respectively, as shown in <Ref>. In addition, the solid domain Ω_s is considered a porous medium, while the fluid domain Ω_f is considered a fracture. The concept of our model is to deal with a sharp interface between Ω_f and Ω_s, while that interface is moved with the help of a phase-field approach.




  §.§.§ Thermo-Fluid-Structure-Interaction problem


In our domain, we start with a thermo-fluid-structure-interaction (TFSI) problem. For the fluid and its temperature within the fluid domain Ω_f, we employ the Boussinesq equation. This is a widely adopted approximation to address nonisothermal flow phenomena such as natural convection, circumventing the need to solve the complete compressible formulation of the Navier-Stokes equations. This approximation holds true when density fluctuations are minor, thereby diminishing the problem's nonlinearity. It assumes that density fluctuations minimally influence the flow field, except for their contribution to buoyancy forces. 

The displacement is modeled by linear elasticity in the solid domain Ω_s while accounting for thermal effects. Consequently, the temperature is solved in the fluid and the solid domain Ω_s ∪Ω_f. In total, we search for the vector-valued velocity v: Ω_f→ℝ^d, the vector-valued displacement u : Ω_s →ℝ^d, the scalar-valued fluid pressure p : Ω_f →ℝ, and the scalar-valued temperature θ : Ω→ℝ. These are determined through the following set of equations: Find the velocity, pressure, temperature, and displacement (v,p,u,θ), such that


    ρ (v·∇) v -
      ∇·_Bou
      - α_θ (θ-θ_0) _f    = 0        in Ω_f, 
    ∇·v   = 0       in Ω_f ,
    
      -∇·_R(u,p,θ)    = _s       in Ω_s, 
    
      -∇·  (κ∇θ)  + v·∇θ   = f_θ      in Ω_f∪Ω_s.



Equations (<ref>),(<ref>), and (<ref>) correspond to a Boussinesq approximation <cit.> with the stress tensor

    _Bouρν (∇v + ∇v^T)  - (p-p_0)I,

where the fluid density is given as ρ >0 and  the fluid viscosity is ν>0. Here, α_θ is the thermal expansion coefficient, and the external force per unit of mass is given as _f.



To consider non-isothermal effects in porous media in the solid domain Ω_s, we assume thermo-poroelasticity <cit.> with the effective solid stress tensor _R defined as

    _R   2 μ() + λ tr(()) I
                - α_B (p-p_0)I
                - 3α_θ K_dr (θ- θ_0)I,

with the Lamé parameters μ,λ>0, the identity matrix
I∈ℝ^d× d and the linearized strain tensor

    () 1/2 (∇ + ∇^T).

Moreover, α_B∈ [0,1] is Biot's coefficient, K_dr = 2/3μ + λ is the bulk modulus, and _s is an external forcing source term. The values p_0 and θ_0 are reference values for the pressure and the temperature, for instance obtained at some initial time or a background state. We remark that specifically, θ_0 will play an important role: θ - θ_0 < 0 means that we inject colder fluid than the existing fluid in the porous media, and θ - θ_0 > 0 means that we inject warmer fluid. The reference value p_0 plays a relatively less important role in this work (since we only consider the injection) and is set to p_0 = 0.

Finally, for the temperature, we consider steady state convection-diffusion heat transfer in the entire domain Ω where κ is the heat conductivity coefficient, with κ_| Ω_f = κ_f and κ_| Ω_s =κ_s, and f_θ is an external forcing source term.






  §.§.§ Interface and Boundary Conditions

The coupling of Ω_s and Ω_f involves integrating the solid domain equations and the fluid domain equations on the interface , resulting in a thermo-fluid-structure interaction (TFSI) problem. Here, the sharp interface between Ω_s and Ω_f is often referred as the fracture interface in our setup and defined as

    Ω_f ∩Ω_s.

First, we  define the following notations

    p_s = pΩ_s,   p_f = pΩ_s,   and  θ_s = θΩ_s, θ_f = θΩ_f,

to specify the pressure and the temperature values for each subdomains. Furthermore, we propose 
continuity of the temperature and pressure, as well as continuity of the normal stresses interface conditions

    p_s = p_f,   θ_s = θ_f,   κ_s ∇θ_s   =  κ_f ∇θ_f ,   _R = _,

where the normal vector  points into the fluid domain Ω_f (fracture region).
We note that the stress tensor _ in Ω_f is defined as

    _  -(p-p_0)I - α_θ(θ-θ_0)I,

where we neglect the displacement in the fluid (fracture) domain. The phase-field formulations of the interface conditions are further discussed  in <Ref>.

Finally, we assume the domain boundary ∂Ω∂Ω_s ∖ contained in the solid boundary, and the system is supplemented by the following (outer) boundary conditions

    u = u_D,  p = p_D,  θ   = θ_D       on ∂Ω,      
                κ_s ∇θ·   = 0       on ∂Ω,

where u_D, p_D, θ_D are corresponding Dirichlet boundary conditions for the displacement, the pressure, and the temperature, respectively. 







  §.§.§ Phase-Field Fracture Problem

As we consider the domain Ω_f as a propagating fracture, we need a method
to compute its changing width and changing length.
In this work, we account for dynamic changes in the subdomains Ω_s and Ω_f by considering fracture propagation due to variations in pressure and temperature.
Here, we employ the phase-field fracture (PFF) approach to track the propagation of fractures. Consequently, the fracture (fluid) domain Ω_f evolves due to the propagation of fractures and the variation in their width and length; see <Ref> for an illustration.

The PFF problem not only help to propagate the fracture and tracks the change of Ω_f, but also provides the fracture width (or crack opening displacement) values to create a sharp geometry representation of the fracture. Fracture width, also known as fracture aperture or crack opening displacement (COD), refers to the perpendicular distance between the two opposing faces of a fracture. In the context of geological formations and fluid-structure interactions, it is a critical parameter that influences the flow of fluids through the fracture. The width of a fracture can change over time due to various factors such as stress, pressure/temperature changes, and the propagation of the fracture itself. In modeling and simulations, accurately determining the fracture width is essential for predicting the behavior of fluids within the fractured medium and for understanding the mechanical properties of the fractured solid.
Thus, one of the main goal of this work is to couple the phase-field approach with the thermo-fluid-structure interaction approach to accurately assess the fracture width,  the deformation of the fracture, and related physics across the fracture interface.

We reconstruct the fractured domain Ω_f from the phase-field variable φ to fully resolve the fracture interface  between the fluid and the intact solid.
Here, the scalar-valued phase-field function, φ: Ω→ [0,1], acts as an indicator function. For example, the fracture domain is defined where φ = 0, and the intact domain is defined where φ = 1. The sharp fracture interface  becomes a diffusive area/volume domain because the phase-field has a diffusive zone where φ∈ (0,1) with a characteristic length scale ε; see <Ref> (a). The phase-field approach solves the fracture problem to propagate the fracture by tracking the phase field values to simulate the fracture propagation as illustrated in <Ref> (b).








  §.§.§ Representation of Fractures

For the coupling of the TFSI and PFF problems, careful definitions of the fracture are required. The TFSI problems involve a sharp interface  between Ω_f and Ω_s, whereas the PFF problems consider a diffusive interface.
First, the classical diffusive phase-field fracture (PFF) is illustrated in <Ref>. This is the diffusive fracture obtained by solving the classical PFF problem. We note the diffusive zone where φ∈ (0,1) around the thin fracture zone (where φ = 0).
Secondly, <Ref> presents the sharp interface, ellipse-shaped fracture. This fracture is obtained by computing the COD values from the PFF shown in <Ref>, and it has no diffusive zone. We utilize this fracture to solve the TFSI problem.
Finally, the sharp interface ellipse fracture is converted back to a diffusive interface ellipse fracture by employing the PFF problem.






 §.§ Overall Coupled Algorithm



In this section, we discuss our proposed algorithm, which considers coupling the non-isothermal TFSI problem to the classical phase-field fracture (PFF) problem. The overall concept can be summarized as follows <cit.>:
enumi-1
    
  * Initialization. This step is only done for the initial time, where we obtain the classical diffusive phase-field fracture φ^0 (initial phase field) by solving the PFF problem with given initial pressure p^0 and temperature θ^0.
    
  * Step 1. Compute the fracture width (crack opening displacement (COD)) values to create a sharp geometry representation of the sharp interface ellipse fracture.
    
  * Step 2. Reconstruct the geometry of the open fluid-filled fracture Ω_s and Ω_f, based on the previously computed COD values.
    
  * Step 3. Solve coupled thermo-fluid-structure interaction (TFSI) problem to get (v, p, θ) in Ω_f and (u, θ) in Ω_s.
    
  * Step 4. Given the pressure p and the temperature θ, solve the PFF problem to obtain the displacement and the phase field (u, φ).  This step is considered to be the prediction of the phase-field fracture, and provides the new fracture domain Ω_f.
    The difference between the initialization and and Step 4 is that the phase-field fracture representation in Step 4 is diffusive interface ellipse fracture whereas the initialization is the thin classical PFF.

A sketch of the algorithm can be seen in <Ref>.
The above algorithm considers multiple couplings between the different variables. First, the temperature and velocity are coupled through the convective and buoyancy terms. Secondly the temperature couples to the solid stress. Finally, the displacement couples to the fluid and temperature equations through the boundary between the fluid and the solid domains.








§ GOVERNING SYSTEM FOR A PHASE-FIELD FRACTURE MODEL



In the previous section, we introduced our overall solution concept. Here, we derive the phase-field fracture approach (PFF problem) with the governing equations used in
this work to model the fluid-filled non-isothermal fracture (Initialization and Step 4).
This PFF problem solves the system to obtain the unknown vector-valued displacement field  and the scalar-valued phase-field function φ. 


In this section begin with a classical formulation and then present the linearized and regularized formulation used in our implementation. 
Crucially, the latter takes into account that the fluid-pressure driving the fracture is only available inside the open crack from an FSI problem. In the following, starting from quasi-static phase-field fracture modeling based on the original work from <cit.>, we recapitulate how pressure interface conditions are included such 
that a pressurized phase-field fracture model is obtained <cit.>.

Let us assume we have a fracture 𝒞⊂ℝ^d-1 in the domain Ω, then note we have Ω_s = Ω\𝒞. To get the correct contributions from traction boundary forces, it is convenient to start from the energy level. Here, traction forces can be described as

    ∫_∂Ω_s_R(,p,θ) · s.

These forces then form the starting point to derive the driving contributions to the phase-field model.


For the following, we note that 𝒞⊂Ω_f in our current setting. However, Ω_f ⊂ℝ^d will later approximate 𝒞 by utilizing the phase-field function. Thus, = ∂Ω_f will also approximate 𝒞, and we obatin 𝒞≈ in the global formulation (<Ref>) with the phase-field function.




 §.§ Interface Conditions





The boundary conditions ∂Ω may be chosen appropriately for a given problem under consideration. For example, we assume homogeneous Dirichlet boundary condition for . However, care must be taken to obtain the correct interface conditions is necessary for , since 𝒞≈, and = Ω_s ∩Ω_f becomes the interface between the fluid (fracture) and the
solid (intact) domain. We follow <cit.> to model the pressure interface
conditions between the surrounding medium and the fracture, and we refer to <cit.> and <cit.> regarding the interface conditions for both the pressure and temperature. The resulting model is a non-isothermal, pressurized phase-field fracture approach. In the following, we provide the mathematical details from prescribing the integral interface conditions and their equivalent formulation as domain integrals.

To include the traction forces (<ref>) in the phase-field model, we assume continuity of the normal stresses on . It then follows that <cit.>:

    - ∫_𝒞_R· s
         = -∫_𝒞_· s 
       = ∫_𝒞 (p-p_0)· s
       + ∫_𝒞α_θ (θ-θ_0) · s 
       =    ∫_Ω_s∇· ((p-p_0) ) 
        -∫_∂Ø (p-p_0)· s 
        +∫_Ω_sα_θ∇· ((θ-θ_0) )
        -∫_∂Øα_θ(θ-θ_0)· s 
       = ∫_Ω_s (∇ (p-p_0) + (p-p_0)∇·)
        + ∫_Ω_sα_θ (∇ (θ-θ_0)
        + (θ-θ_0)∇·),

where Gauss' divergence theorem is applied in (<ref>), and the
homogeneous Dirichlet conditions = 0 on ∂Ω are employed
in (<ref>).




 §.§ 



To transform integrals from subdomains to the global domain Ω, we follow the standard technique in phase-field fracture, and introduces a degradation function, given by

    g(φ) (1-κ) φ^2 + κ,

with the (small) bulk regularization parameter κ>0. We note that g(φ)≈ 0 in the fracture domain (i.e., Ω_f) and g(φ)≈ 1 in the intact domain (i.e., Ω_s).

Next, phase-field models start from lower-dimensional fractures, described by
the Hausdorff measure of the fracture 𝒞. This is then approximated by an Ambrosio-Tortorelli type functional <cit.>:

    G_c H^d-1(𝒞) = G_c ∫_Ω12ε (1-φ)^2 + ε2(∇φ)^2 ,

with the critical energy release rate G_c> 0. Moreover, the classical phase-field fracture model does not allow an open crack to reseal. Thus, the phase field is subject to the crack irreversibility constraint ∂_t φ≤ 0. In our phase-field model, this continuous irreversibility constraint is approximated through a difference quotient by φ≤φ^old. Thus, the phase-field fracture problem is often referred to be in a quasi-static regime. Particularly, the formulation does not contain any time derivatives. Nevertheless, temporal dependence may enter the system through factors such as time-dependent pressure and temperature, and to satisfy the irreversibility constraint.
Due to the quasi-static nature of this problem formulation, we apply some further approximations to arrive at the time-discretized problem.

Let the iteration steps be denoted by τ^m, with index m∈ℕ_0. Then, we have  φ^m ≈φ(τ^m). In other works, this iteration is also known as incremental steps, pseudo-time steps, or time steps. We choose our nomenclature to distinguish this from (time) steps t^n with index n used to advance the overall coupled system of phase-field mesh reconstruction and fluid-structure interaction.


Finally, to formulate the  weak form of global phase-field formulation,
we consider the function spaces W H^1(Ø),
[H^1_0(Ø)]^d and the convex set

    K{w∈ H^1(Ø) |  w≤φ^old≤ 1  a.e. on Ø}∩ L^∞(Ø).

Then, our proposed non-isothermal, pressurized phase-field fracture problem is given be the following definition.


Let the pressure p∈ W^1,∞(Ø), the temperature θ∈ W^1,∞(Ø),
Dirichlet boundary data _D on ∂Ω, and the initial condition
φ(0)φ^0 be given. Furthermore, let the phase-field
regularization parameter >0 and the critical energy release rate G_c> 0
be given. We define the interface driven coupled thermo-phase-field fracture
problem as follows. For the iteration steps m=1,2,3,…,M, find
(,φ)(^m,φ^m) ∈{ + _D}× K,
such that


    [b]
      (g(φ)  _R(), ( ))_Ø
        + (g(φ) (p-p_0), ∇·)_Ø
        + (g(φ) ∇ (p-p_0), )_Ø
    
        + (α_θg(φ) (θ-θ_0), ∇·)_Ø
        + (α_θg(φ)   ∇(θ-θ_0), )_Ø   = 0   ∀∈,
    [b]
        (1-κ) (φ _R():(), ψ-φ)_Ø
        +  2(1-κ) (φ   (p-p_0) ∇·,ψ-φ)_Ø
    
        +  2(1-κ) (φ ∇ (p-p_0) ,ψ)_Ø
        + 2(1-κ)(α_θφ  (θ-θ_0) ∇·, ψ)_Ø
    
        + 2(1-κ)(α_θ  (φ-θ_0)  ∇θ , ψ)_Ø
        + G_c  ((∇φ, ∇ (ψ - φ))_Ø
          -1/(1-φ,ψ-φ)_Ø)
         ≥  0   ∀ψ∈ K.



This formulation uses the above interface law formulated as a domain integral
using the Gauss divergence theorem <cit.>, as derived in
<Ref>.




 §.§ 

Now, the above formulation assumes that the temperature and pressure are given in Ω=(Ω_s ∖𝒞). However, our aim is to couple the temperature and pressure from a thermo-fluid-structure interaction problem to this phase-field model. Consequently, the pressure will only be available in Ω_f and the temperature will be defined in Ω=Ω_f∪Ω_s. Consequently, we need to derive a formulation involving the fracture boundary =∂Ω_f∩∂Ω_s.


To this end, we recall _R from (<ref>), and split it into

    _R = _s - α_B (p-p_0)I - 3α_θ K_dr (θ - θ_0)I     in Ω_s,

where _s = 2 μ() + λ tr(()) I is the linear elasticity part with the displacement.
As in (<ref>), we do not consider the stress
contributions, such that only the pressure and temperature components interact
from Ω_s to , i.e

    - α_B (p-p_0)I - 3α_θ K_dr (θ - θ_0)I   in Ω_s

As shown from (<ref>) to (<ref>), which
transforms the interface integrals to the domain integrals, we perform the similar procedure.

Here, we transform (<ref>) into an interface integral.
To this end, we work again on the energy level with  as variation,
we go backwards the chain and obtain

    - α_B ∫_  (p-p_0)s - 3α_θ K_dr∫_ (θ- θ_0)s.

We note that we are now employing  instead of 𝒞 due to the
given phase-field fracture domain. To transition from the sharp 
interface to the diffusive phase-field representation, we have to include the 
phase-field variable in (<ref>). As for the case φ = 0 the entire integral would vanish,
we add the regularization κ >0 such that the discrete system matrices remain well-posed.


We then have

    - α_B ∫_ g(φ) (p-p_0)s - 3α_θ K_dr∫_ g(φ)(θ- θ_0)s.

Differentiating in  in the direction  and in φ in the direction ψ yields

    - α_B ∫_ g(φ) (p-p_0)s - 3α_θ K_dr∫_ g(φ)(θ- θ_0)s,

and

    - 2(1-κ)α_B ∫_φ (p-p_0)ψs
      - 2(1-κ) 3α_θ K_dr∫_φ(θ- θ_0)ψs,

respectively. We notice that _s remains as the domain integral. Thus,
the pressure and temperature contributions from _R enter in the formulation
as interface integrals on  and the solid stress _s enters as domain
integral contribution. With the above derivation, we finally have the following
non-isothermal, pressurized interface phase-field problem.


  Let the data from <Ref> be given, and  denote the unit normal
  vector pointing into the crack. We define the semi-linearized
  interface driven coupled thermo-phase-field problem as follows. For the
  iteration steps m=1,2,3,…,M, find (,φ)(^m,φ^m)
  ∈{ + _D}× W, such that
  
  
    (g(φ^m-1) _s(), ())_Ø+∫_(1 - α_B)(p - p_0) ·s+∫_(α_θ - 3 α_θ K_dr)  (θ - θ_0) ·s   = 0
          ∀∈,
    [b]
          (1-κ) (φ _s():(), ψ)_Ø
            +  G_c  ((∇φ, ∇ψ)_Ø -1/(1-φ,ψ)_Ø)
            + (γ(φ - φ^m-1)^+,ψ)_Ø
    
          +  2 (1-κ)∫_(1 - α_B)(p - p_0) ·ψs
            +  2 (1-κ) ∫_(α_θ - 3 α_θ K_dr)  (θ - θ_0) ·ψs   =0  ∀ψ∈ W.

  

To the best of our knowledge, a phase-field formulation based on
interface couplings, such as <Ref> has only been used in
<cit.>. This is because the interface formulation seems to contradict
the phase-field concept, where the interface is not known exactly.
In this paper, we extend the idea in <cit.> to consider thermal effects.




In <Ref> we have relaxed the non-linear behavior in the
first term in (<ref>) by using
the approximation φ≈φ^m-1, i.e.,

    g(φ) _R()   ↦   g(φ^m-1) _R().

This follows the extrapolation introduced in <cit.> and is numerically
justified for slowly growing fractures <cit.>.
In the case of fast-growing fractures, this is known to fail <cit.> due
to the time lagging errors. In addition, this approximation could introduce the
fix point iteration error which vanishes with the number of iterations.
Alternatively, fully monolithic schemes <cit.> or an additional
iteration <cit.> must be introduced, to avoid this error.



The second approximation in <Ref> addresses 
irreversibility constraint. We relax this inequality
constraint by considering a simple penalization, see <cit.>
or <cit.>), i.e.,

    φ≤φ^m-1  ↦  γ(φ - φ^m-1)^+.

Here, (x)^+ = x for x>0 and (x)^+ = 0 for x≤ 0,
and where γ>0 is a penalty parameter.




Moreover, we have the following relation for the temperature
interface terms. We have K_dr>1, and so it holds

    (α_θ - 3 α_θ K_dr) < 0.

This means that the injection of colder fluid than the existing fluid in the
porous media, i.e., (θ - θ_0)<0, will cause the fracture to
increase in width and length. On the other hand, for warmer water injection,
i.e., (θ - θ_0)>0, the fracture width will decrease <cit.>.



In this work, we assume α_B = 0. If the full coupling in poroelasticity
with α_B =1 is assumed, one needs to solve for the pressure by utilizing
the poroelastic coupling, then the interface integrals
∫_(1 - α_B)(p - p_0) ·s
and 2(1-κ) ∫_ (1 - α_B)(p - p_0) ·ψs
vanish. However, in our derivation, we still have the pressure contributions
(from the TFSI problem) at the interface as p still contributes to _R.




§ INTERFACE RECONSTRUCTION, REMESHING, AND COUPLED PFF-FSI FRAMEWORK



With the phase-field fracture model presented above, we may compute a fracture 
if the pressure and temperature data are given. However, to obtain these 
quantities from the considered TFSI problem (<ref>), we must first 
obtain the geometry in which the problem is posed.
Following our previous work in <cit.>, we use the geometry
reconstruction approach of the fluid-solid interface presented in <cit.>.
In this approach, we construct a fitted mesh of the open crack and the
surrounding solid are then able to pose our TFSI
problem on this geometry using an interface tracking approach discussed next.



 §.§ Step 1 and Step 2: Computing COD and Remeshing

The geometry reconstruction is based crack opening displacement (COD), or aperture
of the crack <cit.>. This can be computed by

    () = ·≃∫_ℓ^,v() ·∇φ()  s,

where ℓ^,v is a line through  along the vector v
<cit.>. See also <cit.> for a
simplification of the formula when the crack is aligned with a Cartesian axis.

We assume that the centerline of the crack is known, i.e., the line such that
half the fracture width lies on either side of this line, c.f. the dashed line
on the left of <Ref>. With the knowledge of
this line, the crack opening displacements give a set of points on the boundary
of the open crack domain, c.f the blue points on the left of <Ref>.
These can be connected by line segments or higher-order splines, c.f. 
the green line segments on the left of <Ref>. This then forms an
approximation of the crack interface. This geometry can the be remeshed
using an automated meshing tool, resulting in an appropriate mesh for the
a finite element based fluid-structure interaction solver, c.f. the right of 
<Ref>.
We further assume that the number of CODs computed is sufficiently large, such
that the geometry of the open crack is sufficiently well resolved. This can be
achieved by computing O(h^-1) crack opening displacements
along the centerline of the crack.

This re-meshing approach has the advantage
that we can consider propagating and merging cracks, while driving the
crack using accurate quantities from a fluid-structure interaction model
with a resolved interface between the solid and the fluid-filled crack.




  Computing the COD can become numerically unstable near the tips of the crack, c.f. <Ref> below. This is especially the case when iterating between the phase-field fracture and thermo-fluid-structure interaction problems. To avoid this issue, it can become necessary to preprocess the COD values to smooth out oscillations in the COD. Details of this are given below in <Ref>.





 §.§ Step 3: Stationary Thermo-Fluid-Structure Interaction in the Reference Configuration



To obtain a weak formulation suitable for a finite element simulation,
we model the thermo-fluid-structure interaction (TFSI)
problem from <Ref>
in arbitrary Lagrangian-Eulerian (ALE) coordinates. This uses
variational monolithic coupling in a reference configuration
<cit.>.
For our model, we assume a stationary flow inside the fluid-filled crack and,
therefore, consider the stationary Navier-Stokes equations for the fluid.
This goes beyond the previously considered Stokes flow in
<cit.>.

For the TFSI problem, let the domain Ø⊂^d be divided into
a d-dimensional fluid domain , a d-dimensional solid domain
 and a d-1-dimensional interface  between the two, such that
Ø=∪̇∪̇. Furthermore, we also require these domains
in a reference configuration, which we denote by Ø̂, , and . Similarly, we denote by , , and
 the velocity, pressure, deformation and coordinates in the reference
configuration.
In the present setting of a fluid-filled crack, the fluid domain
is the interior of the crack =, the solid is the intact medium
 and the interface is the crack boundary =.


Formulating the problem in the domains  and  leads to the
well-established formulation in ALE coordinates <cit.>. To
obtain a monolithic formulation, we need a transformation _f from the
reference configuration to the physical domain in the fluid-domain. This
transformation is given on the interface  by the structure displacement:

    _f(,t)|_ = +_s(,t)|_.

On the outer boundary of the fluid domain
∂∖, it holds _f=id.
Inside , the only requirement on the transformation is that it should
be as smooth and regular as possible. To this end, we use a harmonic
extension of _s|_ to the fluid domain  and
define _fid+ on . That is,
id() =, i.e., _fid+
= _f(,t) = id()+(,t),
such that

    (∇̂_f,∇̂)_ = 0,  _f=_s on ,  _f=0 on ∂∖.

Consequently, we define a continuous deformation  on all Ω,
which coincides with the solid deformation in  and gives the
appropriate transformation in . Skipping the subscripts and
because _f coincides with the solid transformation _s,
we define on the entire domain :

    (, t) + (, t),  F̂(, t) ∇̂= I+∇̂(, t),  Ĵdet(F̂).


With this transformation into the reference configuration at hand, we
present weak formulation of the stationary thermo-fluid-structure
interaction problem, see also <cit.> for the formal derivation.



Let  be a subspace H^1() with trace zero on
Γ̂^DΓ̂_f^D∪Γ̂_s^D and
L̂ L^2()/ℝ. Furthermore, let
_D, _D∈H^1() be prolongations of the Dirichlet data
for the velocity and deformation, a right-hand side fluid force
∈ L^2() be given. We define the stationary 
thermo-fluid-structure interaction problem as follows.
Find ∈{ + _D}, ∈{+_D}, p̂∈L̂, and θ̂∈, such that

  
    [b]
        (Ĵ_BouF̂^-T,∇̂)_
        + ρ((ĴF̂^-1·∇̂) ,)_
        - (Ĵα_θ (θ̂- θ̂_0) _f,)_
    
          + (Ĵ_RF̂^-T,∇̂)_   = (ρ_f Ĵ,)_   ∀   ∈,
           - (,)_ +
          (α_u ∇̂,∇̂)_   =0    ∀   ∈,
          ( (ĴF̂^-1_f),)_   =0   ∀   ∈L̂,
          (Ĵσ̂_θF̂^-T,∇̂)_
          + (ĴF̂^-1·∇̂θ̂,)_   = (Ĵ f_θ,)_   ∀   ∈,


with the harmonic mesh extension parameter α_u>0, the stress
tensor in the solid as defined in (<ref>) and
Ĵ_RF̂^-T_R.
The thermal stress is σ_θ = k(θ)∇θ with
Ĵσ̂_θF̂^-Tσ_θ.
The temperature dependent ALE fluid stress tensor _Bou is given by

    _Bou -p̂_fI +ρ_f ν(∇̂_f F̂^-1
      + F̂^-T∇̂_f^T),

with the kinematic viscosity ν>0 and the fluid's density ρ_f>0.


Let us comment on the above system in more detail. In (<ref>), 
we combined the momentum equations of the fluid and the solid into one single equation.
This possible with variational-monolithic coupling
in which the interface conditions _f = _s (Dirichlet) and
_Bou = _R (Neumann) are fulfilled in an exact fashion
on the variational level. Moreover, the geometric condition
_f = _s (Dirichlet) is fulfilled as well. The Dirichlet type conditions
are built into the function spaces as usual. The Neumann type condition
cancels out on the interface; see e.g., <cit.>.
In the second equation in (<ref>), the ALE mapping is realized
and for implementational reasons by using globally defined functions, we also work
with = 0 in Ω_s. The third equation (<ref>)
is the mass conservation of the fluid. The last equation (<ref>)
is the weak form of the temperature equation.



§ NUMERICAL TESTS




In this section, we provide several numerical examples to validate and demonstrate the capabilities of the proposed algorithm. The numerical realisation is performed using Netgen/NGSolve[See also <https://ngsolve.org>] <cit.> and the add-on package ngsxfem <cit.>.



 §.§ Numerical Approximation



In total, our algorithm requires four different steps.


    
  Phase-field Fracture A PFF problem must be computed both during initialization and in Step 4 of our algorithm. Specifically, we need to numerically solve the weak formulation of the phase-field fracture problem as described in <Ref>. For this purpose, we use a finite element discretization, where both the phase-field and displacement field spaces are discretized using continuous piecewise linear finite elements. The pressure and temperature are provided as external parameters for this problem.

    To initialize each phase-field computation, we set φ=0 in Ω_f and φ = 1 else, and then solve (<ref>) without the coupling terms, applying a homogeneous Neumann boundary condition. This ensures that the initial condition satisfies the phase-field equation, preventing artificial strength at the crack tips.
    
    
  COD computation 
    In Step 1 of our algorithm, to compute the crack opening displacements (CODs) from the approximated phase-field function and displacement field, we use the unfitted finite element technology provided by ngsxfem to evaluate (<ref>)  over an arbitrary line defined by a level set function. Notably, this level set does not need to be aligned with the mesh.

    
  Geometry reconstruction To construct the geometry from the computed CODs, we use Netgen's OCC (OpenCascade) interface to create a piecewise linear approximation of the interface, which is then meshed. Consequently, this step can also be viewed as an automated CAD model generation process. This is Step 2 of our algorithm.
    
    
  Thermo-fluid-structure-interaction 
    In Step 3 of our algorithm, we numerically solve the weak formulation of the thermo-fluid-structure interaction problem in ALE coordinates as defined in <Ref>. We use the given mesh and discretize the spaces with inf-sup stable elements. Specifically, the velocity space is discretized using continuous piecewise quadratic elements, the pressure with continuous piecewise linear elements, and both the displacement and temperature with continuous piecewise quadratic elements.




 §.§ Example 1:  Temperature Sensitivity Study



As a first example, we consider a basic test inspired by Sneddon's test <cit.>. Here, we do not consider the full model, but only the phase-field problem given in <Ref>, to validate our model as used in Step 0 of our algorithm.



  
Set-up
The set-up for this is as follows. We consider Ω=(-2,2)^2 and an the initial crack (where φ^0=0) is given in  (-0.2, 0.2)×(-,). The material parameters are E=1, ν_s=0.3, G_c=1. The Lamé parameters are then obtained by μ = E/2(1+ν_s) and λ = ν_s E/(1+ν_s)(1-2ν_s). 
The pressure is p=4×10^-2 and the reference pressure is p_0=0. Both the temperature and reference temperature are θ=θ_0=0. The discretizations parameters are κ=10^-10, ε= 2 h, where h is the local mesh size, and γ = 100 h^-2. 



  
Convergence Results
We consider this set up over a series of eight meshes, constructed such that the local mesh size of the crack =  / 100. On the coarsest mesh level, we take =1.28. Each subsequent mesh is constructed by halving both mesh parameters. We compute the crack opening displacement (COD) in the center of the crack and the total crack volume (TCV) and compare the results with those obtained the finest mesh as the reference solution. This reference solution was computed with a total of 2.3× 10^6 degrees of freedom, and therefore close to the limit of the direct solver used to solve the resulting linear systems within each Newton-Step. The results for this can be seen in <Ref>. We see that both the COD and TCV converge with a rate between 1 and 1.5.





  
Temperature effects
To study the effect of the temperature difference on the crack opening displacement, we consider the previous set up on mesh level six with temperatures θ∈{240, 160, 80, 0, -80, -160, -240}. Note that with our chosen material parameters we obtain K_dr = 5/6 and α_θ - 3α_θ K_dr = - 1.5× 10^-5. Therefore, as we consider both constant pressure and temperature in this example, an increase in temperature of 2/3× 10^5 corresponds to a pressure decrease of 1.

We present the results for the change of temperature in <Ref>. As expected, an increase in temperature relative to the reference temperature causes the crack width to shrink, while a decrease in temperature causes the crack to open further <cit.>.





 §.§ Example 2: Fully Coupled Stationary Test Case



As our second example, we consider the fully coupled algorithm in a stationary setting.



  
Set-up
The basic set up is the same as in [sec.numerical_approx:subsec.ex1]Example 1. 
That is, we consider Ω=(-2,2)^2, the initial crack is (-0.2, 0.2)×(-,), the material parameters are E=1, ν_s=0.3, G_c=1, the pressure is p=p^0=4×10^-2, the reference pressure is p_0=0, and both initial temperature and reference temperature are θ=θ_0=0. For the thermo-fluid-structure interaction problem we have κ_f = 0.01, κ_s=1.0 and the force per unit mass is 1_z. The fluid is driven by the body force = 0.2 exp(-1000||^2)_x and the temperature is driven by the external forcing term f̂_θ = 100 exp(-10||^2). The maximal mesh size chosen is h=0.16 and the crack has a local mesh size of = 0.0016. This corresponds to mesh level 3 in the previous example. The remaining discretisations parameters are as in [sec.numerical_approx:subsec.ex1]Example 1.

In Step 4, the temperature from the TFSI problem is then used as the temperature in the PFF problem, i.e., θ^n = θ^n, TFSI. As the pressure in the TFSI problem is normalized to be mean zero, we add this to the initial temperature as the driving pressure in the PFF problem, i.e., p^n= p^0 + p^n, TFSI.



  
Results
The resulting crack opening displacements for eight iterations between the phase-field-fracture and thermo-fluid-structure-interaction problem can be seen in <Ref>, and the resulting total crack volume for every iteration in <Ref>. A visualization of the phase-field, temperature and pressure after the final iteration can be seen in <Ref>. As the FSI temperature is positive (relative to the reference temperature), we have the expected reduction in the fluid/crack volume as seen through the COD plot and TCV values. Furthermore, the FSI pressure is negative in the left half of the domain and positive in the right half of the crack, resulting in the observed skewness of the crack. Furthermore, we see that after eight iterations, there is very little change in the COD.










  
Domain reconstruction
For fine-meshes and and near the singularity of the crack tip, the computation of the COD becomes numerically unstable. Due to the increasing number of necessary CODs computed with decreasing mesh size, this causes a rough boundary of the approximated crack and with each iteration of the coupled loop, this effect increases. To smooth out the crack boundary for the FSI (and subsequent PFF) computations, we process the COD data before the domain is reconstructed. To this end, we use the  function  to compute a least-squares polynomial approximation of the crack boundary. The resulting polynomial values at the O(h^-1) points where the COD was originally computed, plus the roots of this polynomial, are then used to define the crack boundary.

In <Ref>, we illustrate the process on the COD data resulting after four iterations of our coupling loop where the smoothing was not applied using a mesh with =0.04 and =0.0004. Here, we see the visible oscillations towards the crack tip are smoothed out effectively while the shape of the crack is maintained.







 §.§ Example 3: Propagating Crack



We test our novel phase-field model in the case of a propagating crack in two situations: A spatially constant and increasing pressure and finally with a time-dependent temperature and pressure resulting from a TFSI problem. The first example aims to study pressure-driven crack propagation for our novel interface phase-field approach, while the latter aims to study the effect of the temperature coupling to the crack's propagation.



  §.§.§ Example 3a: Pressure driven propagation without TFSI coupling



To test the present interface phase-field model regarding crack propagation, we consider a series of loading steps (pseudo time-steps) and apply an increasing pressure in each iteration, rather than taking the pressure and temperature from a TFSI problem. As the temperature and pressure enter the phase-field model through the same interface integrals, it is sufficient to just consider the pressure in this example.



  
Set-up
The basic set-up is as before. We have, Ω=(-2,2)^2, the initial crack is (-0.2, 0.2)×(-,), material parameters are E=1, ν_s=0.3, G_c=1, the initial pressure is p=4×10^-2 and the temperature is θ=0. The pressure in each iteration of our loop is chosen as p^n = 4× 10^-2 + n × 10^-4 and, we consider a total of 100 loading steps.

We consider a series of three meshes. The coarsest mesh is constructed with = 0.5 and = 0.005. Each subsequent mesh is constructed by halving both mesh parameters. The penalization parameter is again γ = 100 h^-2, but the phase-field regularization parameter is fixed to ε=0.01 for all meshes. The latter corresponds to the previous choice on the coarsest mesh in this example.



  
Results

The resulting horizontal position of the crack tips and the total crack volume from each iteration are shown in <Ref>. We can see the consistent propagation of the crack over each of the four meshes, thereby showing that crack propagation is also feasible for our phase-field model. In <Ref>, we show the phase-field at initialization and for four iterations. Here, we see that the phase-field is well behaved and that the crack increases both in the horizontal and vertical directions.








  §.§.§ Example 3b: Propagation with TFSI coupling



In this example, we consider the fully coupled scheme. 



  
Set-up
The initialization step is as in <Ref>. In the reconstructed domain, the TFSI problem is driven by the fluid and temperature forcing terms

    = 0.02 exp( -400(x-0.1)^2 )_x  and  f̂_θ= -800exp(-||^2).

The heat conductivity parameters are κ_f = 0.005, κ_s=1.0 and the force per unit mass is 0.5_z. The remaining material parameters are as in [sec.numerical_approx:subsec.ex3a]Example 3a.

In each iteration of our coupling loop, the phase-field is then driven by the pressure p^n = 4× 10^-2 + 5n×10^-5 + p^n,TFSI and temperature θ^n=θ^n, TFSI. We again consider 100 loading steps. The remaining phase-field parameters are again chosen as in [sec.numerical_approx:subsec.ex3a]Example 3a.




  
Results
We consider a series of three meshes. First, we see in <Ref> the resulting temperature and pressure from the TFSI computation in the first iteration. Here we see that the temperature is negative and that the TFSI pressure is positive at the right tip of the crack and negative at the left tip of the crack. Due to the cooling of the medium and the higher pressure at the right tip of the crack, we expect the fracture to grow faster than in [sec.numerical_approx:subsec.ex3a]Example 3a. Furthermore, we expect it to grow faster towards the left than the right.



In <Ref>, we see the position of the left and right tips of the crack for each considered mesh. First we note, that the crack does indeed grow faster than in [sec.numerical_approx:subsec.ex3a]Example 3a and that the results are consistent over the series of meshes. Furthermore, we see that on finer meshes, the trajectory becomes smoother. After twenty iterations, the left and right positions of the crack tips are -0.2591, 0.2656, -0.2609, 0.2651 and -0.2571, 0.2598 for the three meshes, respectively. We, therefore, initially observe the expected faster growth toward the right than the left. However, this difference is not very large and after 100 iterations, the left tip has moved further than the right in some cases. We attribute this to numerical error, since this inconsistency occurs earlier on coarser meshes, which we observe to be less stable in <Ref>.







 §.§ Example 4: Two orthogonal cracks with TFSI coupling



In this final example, we consider two orthogonal, connected cracks, illustrating that this approach is applicable to multiple cracks. While this situation is not challenging for phase-field computations, it is more involved with respect to the geometry reconstruction.

To study the effects of including the temperature in our model, we compare two cases here. First, we only couple the TFSI pressure back to the phase-field computation, and denote the resulting phase-field fracture deformation by _h^p. Secondly, we couple both the temperature and pressure back to the phase-field computation, and denote the resulting phase-field fracture deformation by _h^p, θ.

The following example is an extension of the example presented in <cit.>.



  
Set-up
The background domain is Ø=(-2,2)^2, and as before, we consider homogeneous Dirichlet boundary conditions for the displacement and homogeneous Neumann conditions for the phase-field. The initial phase field is given by a flipped "T", i.e.,

    (-0.2, 0.2)×(-h, h)∪(0.2-h, 0.2+h)×(-0.2, 0.2).

As before, the material parameters are as above E=1, ν_s=0.3, and G_c=1, while the initialization pressure and temperature are p=4× 10^-2 and θ=0, respectively.

The thermo-fluid-structure interaction problem in the reconstructed domain is then driven by the forcing terms

    = 0.2exp(-1000( (x - 0.2)^2 + y^2))_x   and  f̂_θ= 100 exp(-10 (x-0.2)^2 -5 y^2).

The reference temperature is set as θ_0 = 0, the heat conductivity parameters are κ_f = 0.01, κ_s=1.0 and the force per unit mass is 1_y.



  
Results

The resulting TFSI temperature in the reconstructed domain can be seen in <Ref>. This temperature is positive and larger inside the crack than the surrounding material. We, therefore, expect the crack to open less when the temperature is coupled to the phase-field fracture model 
in addition to the pressure. Furthermore, we note that while the maximum temperature gets smaller with each smaller mesh size, the results are overall comparable and consistent. This difference appears to be driven by the different of domain from to the geometry reconstruction with changing CODs on each mesh.

The difference between the resulting deformations _h^p - _h^p,θ is shown in <Ref>. We see that this difference points away from the crack, indicating that the deformation _h^p resulting from just the pressure coupling is indeed larger than the deformation _h^p,θ from both pressure and temperature coupling.









§ CONCLUSIONS



In this work, both a new mathematical model and a new numerical approach for thermo-flow-mechanics-fracture are derived. The key idea is to utilize a phase-field approach for fracture opening and fracture propagation. Having the fracture subdomain at hand, a geometry reconstruction approach is employed. This then results in a mesh with which resolves the boundary between the fluid filled crack and the intact solid domains. This allows us to formulate sharp interface problems for the fracture subdomain and the surrounding medium. The resulting framework is a mixture of interface-capturing and interface-tracking approaches that are conveniently combined. It is substantiated for thermo-flow-mechanics-fracture, which is on the one hand thermo-flow-mechanics (THM) phase-field fracture coupled with thermo fluid-structure interaction (TFSI). The latter is prescribed on moving domains (as the fracture moves) for which the arbitrary Lagrangian-Eulerian (ALE) technique is employed. The governing physics are newly developed, specifically the interface conditions for the phase-field sub-problem. These ingredients yield an overall coupling algorithm with four principle steps after an initialization step: fracture width computation (step 1), re-meshing of the reconstructed subdomains (step 2), solving the TFSI problem (step 3), solving the THM phase-field fracture problem (step 4). The algorithmic details and resulting sub-problems are carefully worked out. 

In order to substantiate our new model and new algorithms, we conducted several numerical experiments. Therein, a key component are mesh refinement studies in which the computational robustness is investigated. This was done for the total crack volume (TCV) and the crack opening displacements (COD) as well as the fracture length. All yielded satisfactory findings in view of the complexity of the problem statement. It should be mentioned that specifically the re-meshing of the crack tips required additional algorithmic developments and obtaining computational convergence was a challenge. From the physics point of view, we emphasized the temperature's influence, where there is agreement in the literature that fractures open due to cold water injection and close due to warm water injection. Single fracture as well as two joining fracture were considered. 

The beauty of our overall framework is that different physics can be easily exchanged as long as the interface conditions are correctly modeled. Therefore, our framework presents the opportunity for future extensions, which could include, for example, two-phase flows, or more complicated mechanics. Furthermore, the extension to three spatial dimensions will remain a challenge due to the geometry reconstruction and requires more extensive future work. 




§ DATA AVAILABILITY STATEMENT

The code used to realize the presented results, as well as the resulting raw data, is freely available on github under <https://github.com/hvonwah/repro-tfsi-pff> and is archived on zenodo <https://doi.org/10.5281/zenodo.13685486>.




§ ACKNOWLEDGEMENTS

This material is based upon work supported by the National Science Foundation under Grant No. DMS-1929284 while the HvW was in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the Numerical PDEs: Analysis, Algorithms, and Data Challenges program.
SL acknowledges support within his research stay in May 2024 at the Leibniz University Hannover.
The work of S. Lee was partially supported by the U.S.
National Science Foundation Grant DMS-2208402 and by the U.S. Department of Energy, Office of Science,
Energy Earthshots Initiatives under Award Number DE-SC-0024703.


10
url<#>1urlprefixURL href#1#2#2 #1#1

moore2019utah
J. Moore, J. McLennan, R. Allis, K. Pankow, S. Simmons, R. Podgorney,
  P. Wannamaker, J. Bartley, C. Jones, W. Rickard, The Utah Frontier
  Observatory for p in Geothermal Energy (FORGE): An international
  laboratory for enhanced geothermal system technology development, in: 44th
  Workshop on Geothermal Reservoir Engineering, Stanford University, 2019, pp.
  11–13.

olasolo2016enhanced
P. Olasolo, M. C. Juárez, M. P. Morales, I. A. Liarte, et al., Enhanced
  geothermal systems (EGS): A review, Renewable Sustainable Energy Rev. 56
  (2016) 133–144.
https://doi.org/10.1016/j.rser.2015.11.031
  doi:10.1016/j.rser.2015.11.031.

lu2018global
S.-M. Lu, A global review of enhanced geothermal system (EGS), Renewable
  Sustainable Energy Rev. 81 (2018) 2902–2921.
https://doi.org/10.1016/j.rser.2017.06.097
  doi:10.1016/j.rser.2017.06.097.

mcclure2014investigation
M. W. McClure, R. N. Horne, An investigation of stimulation mechanisms in
  Enhanced Geothermal Systems, Int. J. Rock Mech. Min. 72 (2014)
  242–260.
https://doi.org/10.1016/j.ijrmms.2014.07.011
  doi:10.1016/j.ijrmms.2014.07.011.

caulk2016experimental
R. A. Caulk, E. Ghazanfari, J. N. Perdrial, N. Perdrial, Experimental
  investigation of fracture aperture and permeability change within Enhanced
  Geothermal Systems, Geothermics 62 (2016) 12–21.
https://doi.org/10.1016/j.geothermics.2016.02.003
  doi:10.1016/j.geothermics.2016.02.003.

zhang2022thermal
X. Zhang, Z. Li, X. Wang, H. Wang, B. Li, Y. Niu, Thermal effect on the
  fracture behavior of granite using acoustic emission and digital image
  correlation: An experimental investigation, Theor. Appl. Fract. Mec. 121
  (2022) 103540.
https://doi.org/10.1016/j.tafmec.2022.103540
  doi:10.1016/j.tafmec.2022.103540.

donahue1972crack
R. J. Donahue, H. M. Clark, P. Atanmo, R. Kumble, A. J. McEvily, Crack opening
  displacement and the rate of fatigue crack growth, Int. J. Fract. Mech. 8
  (1972) 209–219.
https://doi.org/10.1007/BF00703882
  doi:10.1007/BF00703882.

burdekin1966crack
F. M. Burdekin, D. E. W. Stone, The crack opening displacement approach to
  fracture mechanics in yielding materials, J. Strain Anal. Eng. Des. 1 (2)
  (1966) 145–153.
https://doi.org/10.1243/03093247V012145
  doi:10.1243/03093247V012145.

zhou2022thermal
L. Zhou, W. Gao, L. Yu, Z. Zhu, J. Chen, X. Wang, Thermal effects on fracture
  toughness of cracked straight-through Brazilian disk green sandstone and
  granite, J. Rock Mech. Geotech. Eng. 14 (5) (2022) 1447–1460.
https://doi.org/10.1016/j.jrmge.2022.02.016
  doi:10.1016/j.jrmge.2022.02.016.

BourFraMar00
B. Bourdin, G. A. Francfort, J.-J. Marigo, Numerical experiments in revisited
  brittle fracture, J. Mech. Phys. Solids 48 (4) (2000) 797–826.
https://doi.org/10.1016/S0022-5096(99)00028-9
  doi:10.1016/S0022-5096(99)00028-9.

KuMue10
C. Kuhn, R. Müller, A continuum phase field model for fracture, Eng. Fract.
  Mech. 77 (18) (2010) 3625–3634, Computational Mechanics in Fracture and
  Damage: A Special Issue in Honor of Prof. Gross.
https://doi.org/10.1016/j.engfracmech.2010.08.009
  doi:10.1016/j.engfracmech.2010.08.009.

MieWelHof10a
C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field
  models of fracture: Variational principles and multi-field FE
  implementations, Internat. J. Numer. Methods Engrg. 83 (10) (2010)
  1273–1311.
https://doi.org/10.1002/nme.2861
  doi:10.1002/nme.2861.

MieWelHof10b
C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent
  crack propagation: Robust algorithmic implementation based on operator
  splits, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2765–2778.
https://doi.org/10.1016/j.cma.2010.04.011
  doi:10.1016/j.cma.2010.04.011.

BoVeScoHuLa12
M. J. Borden, C. V. Verhoosel, M. A. Scott, T. J. R. Hughes, C. M. Landis, A
  phase-field description of dynamic brittle fracture, Comput. Methods Appl.
  Mech. Engrg. 217 (2012) 77–95.
https://doi.org/10.1016/j.cma.2012.01.008
  doi:10.1016/j.cma.2012.01.008.

AmGeraLoren15
M. Ambati, T. Gerasimov, L. De Lorenzis, A review on phase-field models of
  brittle fracture and a new fast hybrid formulation, Comput. Mech. 55 (2)
  (2015) 383–405.
https://doi.org/10.1007/s00466-014-1109-y
  doi:10.1007/s00466-014-1109-y.

ARRIAGA201833
M. Arriaga, H. Waisman, Stability analysis of the phase-field method for
  fracture with a general degradation function and plasticity induced crack
  generation, Mech. Mater. 116 (2018) 33–48, iUTAM Symposium on Dynamic
  Instabilities in Solids.
https://doi.org/10.1016/j.mechmat.2017.04.003
  doi:10.1016/j.mechmat.2017.04.003.

SARGADO2018458
J. M. Sargado, E. Keilegavlen, I. Berre, J. M. Nordbotten, High-accuracy
  phase-field models for brittle fracture based on a new family of degradation
  functions, J. Mech. Phys. Solids 111 (2018) 458–489.
https://doi.org/10.1016/j.jmps.2017.10.015
  doi:10.1016/j.jmps.2017.10.015.

WheWiLee20
M. F. Wheeler, T. Wick, S. Lee, IPACS: Integrated Phase-Field Advanced Crack
  Propagation Simulator. An adaptive, parallel, physics-based-discretization
  phase-field framework for fracture propagation in porous media, Comput.
  Methods Appl. Mech. Engrg. 367 (2020) 113124.
https://doi.org/10.1016/j.cma.2020.113124
  doi:10.1016/j.cma.2020.113124.

BourFraMar08
B. Bourdin, G. A. Francfort, J.-J. Marigo, The variational approach to
  fracture, J. Elasticity 91 (1–3) (2008) 1–148.
https://doi.org/10.1007/s10659-007-9107-3
  doi:10.1007/s10659-007-9107-3.

WNN20
J.-Y. Wu, V. P. Nguyen, C. T. Nguyen, D. Sutula, S. Sinaie, S. P. A. Bordas,
  Phase-field modeling of fracture, Advances in Applied Mechanics, Elsevier,
  2020, Ch. 1, pp. 1–183.
https://doi.org/10.1016/bs.aams.2019.08.001
  doi:10.1016/bs.aams.2019.08.001.

Wi20
T. Wick, Multiphysics Phase-Field Fracture, Vol. 28 of Radon Series on
  Computational and Applied Mathematics, De Gruyter, Berlin, Boston, 2020.
https://doi.org/10.1515/9783110497397
  doi:10.1515/9783110497397.

HEIDER2021107881
Y. Heider, A review on phase-field modeling of hydraulic fracturing, Eng.
  Fract. Mech. 253 (2021) 107881.
https://doi.org/10.1016/j.engfracmech.2021.107881
  doi:10.1016/j.engfracmech.2021.107881.

DiLiWiTy22
P. Diehl, R. Lipton, T. Wick, M. Tyagi, A comparative review of peridynamics
  and phase-field models for engineering fracture mechanics, Comput. Mech. 69
  (2022) 1259–1293.
https://doi.org/10.1007/s00466-022-02147-0
  doi:10.1007/s00466-022-02147-0.

LeeWheWi16
S. Lee, M. F. Wheeler, T. Wick, Pressure and fluid-driven fracture propagation
  in porous media using an adaptive finite element phase field model, Comput.
  Methods Appl. Mech. Engrg. 305 (2016) 111–132.
https://doi.org/10.1016/j.cma.2016.02.037
  doi:10.1016/j.cma.2016.02.037.

lee2016phase
S. Lee, A. Mikelić, M. F. Wheeler, T. Wick, Phase-field modeling of
  proppant-filled fractures in a poroelastic medium, Comput. Methods Appl.
  Mech. Engrg. 312 (2016) 509–541.
https://doi.org/10.1016/j.cma.2016.02.008
  doi:10.1016/j.cma.2016.02.008.

MDB99
N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth
  without remeshing, Internat. J. Numer. Methods Engrg. 46 (1) (1999) 131–150.
https://doi.org/10.1002/(sici)1097-0207(19990910)46:1<131::aid-nme726>3.0.co;2-j
  doi:10.1002/(sici)1097-0207(19990910)46:1<131::aid-nme726>3.0.co;2-j.

zi2003new
G. Zi, T. Belytschko, New crack-tip elements for XFEM and applications to
  cohesive cracks, Internat. J. Numer. Methods Engrg. 57 (15) (2003)
  2221–2240.
https://doi.org/10.1002/nme.849 doi:10.1002/nme.849.

duarte2001generalized
C. A. Duarte, O. N. Hamzeh, T. J. Liszka, W. W. Tworzydlo, A generalized finite
  element method for the simulation of three-dimensional dynamic crack
  propagation, Comput. Methods Appl. Mech. Engrg. 190 (15-17) (2001)
  2227–2262.
https://doi.org/10.1016/S0045-7825(00)00233-4
  doi:10.1016/S0045-7825(00)00233-4.

FB10
T.-P. Fries, T. Belytschko, The extended/generalized finite element method:
  An overview of the method and its applications, Internat. J. Numer. Methods
  Engrg. 84 (3) (2010) 253–304.
https://doi.org/10.1002/nme.2914
  doi:10.1002/nme.2914.

LWW17
S. Lee, M. F. Wheeler, T. Wick, Iterative coupling of flow, geomechanics and
  adaptive phase-field fracture including level-set crack width approaches, J.
  Comput. Appl. Math. 314 (2017) 40–60.
https://doi.org/10.1016/j.cam.2016.10.022
  doi:10.1016/j.cam.2016.10.022.

YNK20
K. Yoshioka, D. Naumov, O. Kolditz, On crack opening computation in variational
  phase-field models for fracture, Comput. Methods Appl. Mech. Engrg. 369
  (2020) 113210.
https://doi.org/10.1016/j.cma.2020.113210
  doi:10.1016/j.cma.2020.113210.

WaWi24_RINAM
H. von Wahl, T. Wick, A coupled high-accuracy phase-field fluid–structure
  interaction framework for stokes fluid-filled fracture surrounded by an
  elastic medium, Results Appl. Math. 22 (2024) 100455.
https://doi.org/10.1016/j.rinam.2024.100455
  doi:10.1016/j.rinam.2024.100455.

WaWi23_CMAME
H. von Wahl, T. Wick, A high-accuracy framework for phase-field fracture
  interface reconstructions with application to Stokes fluid-filled fracture
  surrounded by an elastic medium, Comput. Methods Appl. Mech. Engrg. 415
  (2023) 116202.
https://doi.org/10.1016/j.cma.2023.116202
  doi:10.1016/j.cma.2023.116202.

NoiiWi19
N. Noii, T. Wick, A phase-field description for pressurized and non-isothermal
  propagating fractures, Comput. Methods Appl. Mech. Engrg. 351 (2019)
  860–890.
https://doi.org/10.1016/j.cma.2019.03.058
  doi:10.1016/j.cma.2019.03.058.

HEIDER2018116
Y. Heider, S. Reiche, P. Siebert, B. Markert, Modeling of hydraulic fracturing
  using a porous-media phase-field approach with reference to experimental
  data, Eng. Fract. Mech. 202 (2018) 116–134.
https://doi.org/10.1016/j.engfracmech.2018.09.010
  doi:10.1016/j.engfracmech.2018.09.010.

NgHeiMa23
C.-L. Nguyen, Y. Heider, B. Markert, A non-isothermal phase-field hydraulic
  fracture modeling in saturated porous media with convection-dominated heat
  transport, Acta Geotech. 50 (6) (2023) 821–833.
https://doi.org/10.1007/s11440-023-01905-5
  doi:10.1007/s11440-023-01905-5.

SUH2021114182
H. S. Suh, W. Sun, Asynchronous phase field fracture model for porous media
  with thermally non-equilibrated constituents, Comput. Methods Appl. Mech.
  Engrg. 387 (2021) 114182.
https://doi.org/10.1016/j.cma.2021.114182
  doi:10.1016/j.cma.2021.114182.

Dai2024
Y. Dai, B. Hou, S. Lee, T. Wick, A thermal–hydraulic–mechanical–chemical
  coupling model for acid fracture propagation based on a phase-field method,
  Rock Mech. Rock Eng. (2024).
https://doi.org/10.1007/s00603-024-03769-x
  doi:10.1007/s00603-024-03769-x.

LIU2024117165
Y. Liu, K. Yoshioka, T. You, H. Li, F. Zhang, A phase-field fracture model in
  thermo-poro-elastic media with micromechanical strain energy degradation,
  Comput. Methods Appl. Mech. Engrg. 429 (2024) 117165.
https://doi.org/10.1016/j.cma.2024.117165
  doi:10.1016/j.cma.2024.117165.

lee4920845phase
S. Lee, M. Wheeler, T. Wick, https://ssrn.com/abstract=4920845A
  phase-field diffraction model for thermo-hydro-mechanical propagating
  fractures (Aug. 2024).
<https://ssrn.com/abstract=4920845>

LoBo96
S. A. Lorca, J. L. Boldrini, Stationary solutions for generalized Boussinesq
  models, J. Differ. Equ. 124 (2) (1996) 389–406.
https://doi.org/10.1006/jdeq.1996.0016
  doi:10.1006/jdeq.1996.0016.

FARHAT1991349
C. Farhat, K. C. Park, Y. Dubois-Pelerin, An unconditionally stable staggered
  algorithm for transient finite element analysis of coupled thermoelastic
  problems, Comput. Methods Appl. Mech. Engrg. 85 (3) (1991) 349–365.
https://doi.org/10.1016/0045-7825(91)90102-C
  doi:10.1016/0045-7825(91)90102-C.

Coussy2004
O. Coussy, Poromechanics, Wiley, 2004.
https://doi.org/10.1002/0470092718
  doi:10.1002/0470092718.

mayeli2021buoyancy
P. Mayeli, G. J. Sheard, Buoyancy-driven flows beyond the Boussinesq
  approximation: A brief review, Int. Commun. Heat Mass 125 (2021) 105316.
https://doi.org/10.1016/j.icheatmasstransfer.2021.105316
  doi:10.1016/j.icheatmasstransfer.2021.105316.

MiWheWi13a
A. Mikelić, M. F. Wheeler, T. Wick,
  www.oden.utexas.edu/media/reports/2013/1315.pdfA phase-field approach
  to the fluid filled fracture surrounded by a poroelastic medium, iCES Report
  13-15 (Jun. 2013).
<www.oden.utexas.edu/media/reports/2013/1315.pdf>

MWW19
A. Mikelić, M. F. Wheeler, T. Wick, Phase-field modeling through
  iterative splitting of hydraulic fractures in a poroelastic medium, GEM -
  Int. J. Geomath. 10 (1) (Jan. 2019).
https://doi.org/10.1007/s13137-019-0113-y
  doi:10.1007/s13137-019-0113-y.

AmTo90
L. Ambrosio, V. M. Tortorelli, Approximation of functionals depending on jumps
  by elliptic functionals via γ-convergence, Comm. Pure Appl. Math.
  43 (8) (1990) 999–1036.
https://doi.org/10.1002/cpa.3160430805
  doi:10.1002/cpa.3160430805.

AmTo92
L. Ambrosio, V. M. Tortorelli, On the approximation of free discontinuity
  problems, Boll. Un. Mat. Ital. 6 (1992) 105–123.

HWW15
T. Heister, M. F. Wheeler, T. Wick, A primal-dual active set method and
  predictor-corrector mesh adaptivity for computing fracture propagation using
  a phase-field approach, Comput. Methods Appl. Mech. Engrg. 290 (2015)
  466–495.
https://doi.org/10.1016/j.cma.2015.03.009
  doi:10.1016/j.cma.2015.03.009.

Wic17
T. Wick, An error-oriented Newton/inexact augmented Lagrangian approach for
  fully monolithic phase-field fracture propagation, SIAM J. Sci. Comput.
  39 (4) (2017) B589–B617.
https://doi.org/10.1137/16m1063873
  doi:10.1137/16m1063873.

KMW23
L. Kolditz, K. Mang, T. Wick, A modified combined active-set Newton method
  for solving phase-field fracture into the monolithic limit, Comput. Methods
  Appl. Mech. Engrg. 414 (2023) 116170.
https://doi.org/10.1016/j.cma.2023.116170
  doi:10.1016/j.cma.2023.116170.

MWT15
A. Mikelić, M. F. Wheeler, T. Wick, A quasi-static phase-field approach
  to pressurized fractures, Nonlinearity 28 (5) (2015) 1371–1399.
https://doi.org/10.1088/0951-7715/28/5/1371
  doi:10.1088/0951-7715/28/5/1371.

TrSeNg13
D. Tran, A. T. Settari, L. Nghiem, Predicting growth and decay of
  hydraulic-fracture width in porous media subjected to isothermal and
  nonisothermal flow, SPE J. 18 (4) (2013) 781–794.
https://doi.org/10.2118/162651-PA
  doi:10.2118/162651-PA.

CHUKWUDOZIE2019957
C. Chukwudozie, B. Bourdin, K. Yoshioka, A variational phase-field model for
  hydraulic fracturing in porous media, Comput. Methods Appl. Mech. Engrg. 347
  (2019) 957–982.
https://doi.org/10.1016/j.cma.2018.12.037
  doi:10.1016/j.cma.2018.12.037.

HrTu06a
J. Hron, S. Turek, A monolithic FEM/Multigrid solver for ALE formulation of
  fluid structure with application in biomechanics, Vol. 53, Springer, Berlin,
  Heidelberg, 2006, pp. 146–170.
https://doi.org/10.1007/3-540-34596-5_7
  doi:10.1007/3-540-34596-5_7.

Du07
T. Dunne, Adaptive finite element approximation of fluid-structure interaction
  based on Eulerian and arbitrary Lagrangian-Eulerian variational
  formulations, Ph.D. thesis, University of Heidelberg (2007).
https://doi.org/10.11588/heidok.00007944
  doi:10.11588/heidok.00007944.

Wi11_phd
T. Wick, Adaptive Finite Element Simulation of Fluid-Structure
  Interaction with Application to Heart-Valve Dynamics, Ph.D. thesis,
  University of Heidelberg (2011).
https://doi.org/10.11588/heidok.00012992
  doi:10.11588/heidok.00012992.

Ri17_fsi
T. Richter, Fluid-structure interactions: Models, analysis, and finite
  elements, Springer, Cham, 2017.
https://doi.org/10.1007/978-3-319-63970-3
  doi:10.1007/978-3-319-63970-3.

HuLiZi81
T. J. R. Hughes, W. K. Liu, T. Zimmermann, Lagrangian-Eulerian finite
  element formulation for incompressible viscous flows, Comput. Methods Appl.
  Mech. Engrg. 29 (1981) 329–349.
https://doi.org/10.1016/0045-7825(81)90049-9
  doi:10.1016/0045-7825(81)90049-9.

DoGiuHa82
J. Donea, S. Giuliani, J. P. Halleux, An arbitrary Lagrangian-Eulerian
  finite element method for transient dynamic fluid-structure interactions,
  Comput. Methods Appl. Mech. Engrg. 33 (1982) 689–723.
https://doi.org/10.1016/0045-7825(82)90128-1
  doi:10.1016/0045-7825(82)90128-1.

RiWi10
T. Richter, T. Wick, Finite elements for fluid-structure interaction in ALE
  and fully Eulerian coordinates, Comput. Methods Appl. Mech. Engrg. 199
  (2010) 2633–2642.
https://doi.org/10.1016/j.cma.2010.04.016
  doi:10.1016/j.cma.2010.04.016.

Sch97
J. Schöberl, NETGEN an advancing front 2D/3D-mesh generator based on
  abstract rules, Comput. Vis. Sci. 1 (1) (1997) 41–52.
https://doi.org/10.1007/s007910050004
  doi:10.1007/s007910050004.

Sch14
J. Schöberl, C++11 implementation of finite elements in NGSolve, Tech.
  Rep. ASC Report No. 30/2014 (Sep. 2014).

LHPvW21
C. Lehrenfeld, F. Heimann, J. Preuß, H. von Wahl,
  github.com/ngsxfem/ngsxfem: Add-on to ngsolve for
  geometrically unfitted finite element discretizations, J. Open Source Softw.
  6 (64) (2021) 3237.
https://doi.org/10.21105/joss.03237
  doi:10.21105/joss.03237.
<github.com/ngsxfem/ngsxfem>

Sne46
I. N. Sneddon, The distribution of stress in the neighbourhood of a crack in an
  elastic solid, Proc. R. Soc. A 187 (1009) (1946) 229–260.
https://doi.org/10.1098/rspa.1946.0077
  doi:10.1098/rspa.1946.0077.

SneddLow69
I. N. Sneddon, M. Lowengrub, Crack problems in the classical theory of
  elasticity, SIAM series in Applied Mathematics, John Wiley and Sons,
  Philadelphia, 1969.

morrow2001permeability
C. A. Morrow, D. E. Moore, D. A. Lockner, Permeability reduction in granite
  under hydrothermal conditions, J. Geophys. Res. Solid Earth 106 (B12) (2001)
  30551–30560.
https://doi.org/10.1029/2000JB000010
  doi:10.1029/2000JB000010.

yasuhara2006evolution
H. Yasuhara, A. Polak, Y. Mitani, A. S. Grader, P. M. Halleck, D. Elsworth,
  Evolution of fracture permeability through fluid–rock reaction under
  hydrothermal conditions, Earth Planet. Sc. Lett. 244 (1-2) (2006) 186–200.
https://doi.org/10.1016/j.epsl.2006.01.046
  doi:10.1016/j.epsl.2006.01.046.

hardin1982measuring
E. Hardin, N. Barton, M. Voegele, M. Board, R. Lingle, H. Pratt, W. Ubbes,
  Measuring the thermomechanical and transport properties of a rockmass using
  the heated block test, in: ARMA US Rock Mechanics/Geomechanics Symposium,
  ARMA, 1982, pp. ARMA–82.

rutqvist2008analysis
J. Rutqvist, B. Freifeld, K.-B. Min, D. Elsworth, Y. Tsang, Analysis of
  thermally induced changes in fractured rock permeability during 8 years of
  heating and cooling at the yucca mountain drift scale test, Int. J. Rock
  Mech. Min. 45 (8) (2008) 1373–1389.
https://doi.org/10.1016/j.ijrmms.2008.01.016
  doi:10.1016/j.ijrmms.2008.01.016.





