A non-isothermal phase-field hydraulic fracture modeling in saturated porous media with convection-dominated heat transport

This research aims to extend the isothermal continuum mechanical modeling framework of hydraulic fracturing in porous materials to account for the non-isothermal processes. Whereas the theory of porous media is used for the macroscopic material description, the phase-field method is utilized for modeling the crack initiation and propagation. We proceed in this study from a two-phase porous material consisting of thermomechanically interacting pore fluid and solid matrix. The heat exchange between the fluid in the crack and the surrounding porous environment through the diffusive fracture edges is carefully studied, and new formulations here are proposed. Besides, temperature-dependent solid and fluid material parameters are taken into account, which is of particular importance in connection with fluid viscosity and its effect on post-cracking pressure behavior. This continuum mechanical treatment results in strongly coupled partial differential equations of the mass, the momentum, and the energy balance of the thermally non-equilibrated constituents. Using the finite element method, two-dimensional initial-boundary-value problems are presented to show, on the one hand, the stability and robustness of the applied numerical algorithm in solving the emerged strongly coupled problem in the convection-dominated heat transport state. On the other hand, they show the capability of the modeling scheme in predicting important instances related to hydraulic fracturing and the role of the temperature field in this process. Additionally, they show the importance of using stabilization techniques, such as adding an artificial thermo-diffusivity term, to mitigate temperature fluctuations at high flow velocity.


Introduction
The subject of hydraulic fracturing and the existence of fracture networks in porous materials are of major importance in geothermal energy systems with naturally or artificially fractured underground rock layers.In geothermal systems with rock layers at high temperatures, a significant amount of thermal energy can be gained and used.Hence, depending on the depth of the rock layers and their temperature, this energy can be used for heating or for electricity generation using geothermal power plants.One of the major challenges in this regard is the low permeability of the underlying rock layers, which deteriorates the efficiency of geothermal systems.This challenge can be overcome by applying hydraulic fracturing leading to enhanced geothermal systems (EGS), see, e.g., [81] for review.In this connection, the focus of the underlying work is on the modeling of hydraulic fracturing to enhance the reservoir permeability and on the EGS operation that includes thermal energy exchange between the low-temperature injection fluid and the high-temperature surrounding deep rock layers.In particular, the hydraulic fracturing within non-isothermal conditions will be modeled using the phase-field method (PFM), see, e.g., [19,48,52,79,92], whereas the heat exchange with the surrounding porous rocks will be modeling using the fundamentals of the theory of porous media (TPM) for nonisothermal problems, see, e.g., [29,57,94,95].

Literature review for isothermal phase-field hydraulic fracture
While few research works addressed the PFM-based hydraulic fracture modeling under non-isothermal conditions, many works studied this under isothermal conditions as reviewed in [48].The starting point in these studies is the definition of a scalar-valued and space-time-dependent phase-field variable in the range [0, 1], which allows for differentiation between the cracked and the intact regions of the domain.With the PFM, the sharp edges of the crack are approximated by diffusive ones, whereas the width of this region is controlled via an internal length scale, see, e.g., [6,12,16,17,21,49,53,58,62,63,70,73,74] among others.In this context, it is worth mentioning that although the PFM is intensively applied for fracture modeling, it is also widely used in the simulation of many other engineering applications, as phase-change materials, see, e.g., [4,[94][95][96]104], for references.
In the development of the PFM, [45] presented pioneering studies for brittle materials' fracture, including the definition of the crack surface energy and the critical energy release rate.Later, [55] extended Griffith's approach to ductile materials and proposed the concept of the strain energy release rate and crack toughness.Subsequent studies by [12,13,40] proposed an energetic framework to capture the relationship between the diffusive interface approach, utilized in physics to realize the boundaries between different phases of heterogeneous materials, and fracture mechanics.Significant advancements to the PFM for fracture modeling, including the mathematical formulations within a thermodynamical framework, large deformations, solution algorithms, applications, and model validation, can be found in, e.g., [6,58,70,77,82,83,103] among others.
In talking about hydraulic fracture modeling, the work of [14] counts as one of the first efforts to utilize the PFM to simulate hydraulic fracturing within a variational framework.The PFM was thereafter used by many research groups to capture fluid-induced cracking in porous and solid materials as can be found in, e.g., [30,33,49,67,72,98,100,101], among others.Many of these works combined the PFM together with the continuum theory of porous media (TPM) [28] or Biot's poroelasticity [9] to capture the onset of fracture and the accompanying hydro-mechanical processes.For example, Markert and Heider [67] proposed a novel coupled, multifield PFM-TPM fracture simulation framework, which was advanced in [49] to capture irreversible changes in the fractured porous materials, like the volume fractions and the permeability.In line with this, Pillai et al. [84] implemented a PFM for dynamic hydraulic fracture in an Abaqus user-element subroutine, which captured also the domain heterogeneity.In this connection, a model validation via comparison with experimental data on hydraulic fractures was presented in [52] and an extension to unsaturated porous media fracture is found in [17,51].[30] introduced a thermodynamical framework within the TPM for modeling hydraulic fracturing under dynamic conditions using the PFM.This approach was extended in [31], in which they proposed a crack-opening indicator function in connection with the fluid flow in the crack.An extension of the TPM-PFM approach towards ductile fracture can be found in [85].Important advancements in the context of PFM of hydraulic fractures in connection with Biot's theory can be dedicated to [2,3,60,72,75,90,91,101] among others.In these research works, different characteristics of hydraulic fracturing in two-and three-dimensional settings could be captured.To mention some, important advancements were done with regard to the crack width computation.Additionally, advanced numerical algorithms, like the localglobal approach or the fixed-stress iterative coupling, were also implemented.More details and references will be provided in later sections.

Literature review for non-isothermal phasefield fracture and heat exchange processes
An extension of the PFM for hydraulic fracture modeling in poroelasticity to account for non-isothermal conditions can be found in [79], where they proceeded from quasistatic and small strain assumptions.A particular focus in this work was on the thermal exchange and force balance at the interface between the crack and the domain.Additionally, the numerical solution considered a robust semismooth Newton approach and a crack-oriented predictorcorrector spatial adaptivity method.Extension of the TPM-PFM approach for hydraulic fracture within non-isothermal, thermal non-equilibrium settings can be found in [78], where the heat exchange between the cold fluid in the crack and the hot ambient could be captured.Suh and Sun [92] developed an asynchronous FE solver to capture the crack propagation and fluid flow in porous media at the mesoscale under non-isothermal conditions.They considered thermally non-equilibrium states of the constituents.An operator-split time integration algorithm was also applied to compute the displacement, pressure, phase field, and temperature of each phase in an asynchronous way.
Apart from the PFM approaches to model crack propagation and the accompanying processes, the heat exchange between the fluid in the crack or in a tube and the surrounding domain has intensively been studied in the literature, both experimentally and numerically.For instance, Ramey [87] provided an approximate, analytical solution to the wellbore heat-transmission problem.The proposed solution proceeded from a steady-state heat transfer in the wellbore, whereas the heat transfer to the surrounding ambient is unsteady radial.The solution allows estimating the temperature of flowing fluid, tubing, and casing as a function of depth and time.Within the thermo-poroelasticity theory, Wang and Papamichos [99] studied the coupled thermo-hydro-mechanical processes in a circular wellbore or spherical cavity, in which a constant fluid flow rate and temperature were considered.They developed analytical solutions to estimate the temperature, pressure, and stresses in the area around the wellbore and cavity for low permeability ambient and negligible heat convective processes.They showed in their study that injection of warm fluid delays the fracture initiation around the cavity due to emerging of tensile effective tangential stresses.Ghassemi et al. [43] carried out an analytical study within poro-thermoelasticity to investigate in EGS the effect of injection fluid temperature (mostly cold water) on the thermal energy exchange and crack opening.They studied the effect of effective stresses and thermoelastic effects in the porous ambient on the reduction of the crack aperture.They showed that these effects are maximum near the inject point and minimum near the extraction point.In line with this, an alternative analytical solution for the crack width within a continuum thermo-poroelastic theory was presented [97].Recently, Frank et al. [41] carried out an experimental study on a single cracked rock with flowing fluid at a high rate to understand the heat transfer process.They investigated different geometrical and material parameters that affect the transient and steadystate heat transfer and proposed analytical formulation to describe the aforementioned processes.Proceeding with a TPM continuum mechanical model for porous materials under non-isothermal conditions, [57] derived a thermodynamically consistent biphasic model to describe the functionality of EGS systems.In this macroscopic model, the cracks were not explicitly modeled and different temperatures of the elastic solid and the viscous fluid were considered.They finally solved an initial boundary-value problem that relied on a circulation test in a geothermal plant in Soultz-sous-Fore ˆts (France).
In the current work, the modeling of the fully saturated porous biphasic material with materially incompressible solid phase and compressible pore fluid under nonisothermal conditions is based on the theory of porous media (TPM) principles, see, e.g., [8,22,23,32,47].The continuum-mechanical description of the porous media fracture problem results in a set of coupled partial differential equations (PDEs), which class and characteristics have been investigated in, e.g., [37,65,68,80].In connection with the numerical implementation, the finite element method (FEM) together with the implicit time integration scheme in a staggered manner to be applied.Of particular importance in the implementation are the stability challenges in the nearly incompressible fluid limit or at the crack-domain interface, as discussed in, e.g., [51].

Highlights and content overview
With the aim of deeply understanding the hydraulic fracturing process in saturated porous materials under nonisothermal conditions, the outlines of the underlying work can be summarized in the following points: (1) Extending the isothermal phase-field porous media fracturing approach, developed within the theory of porous media (TPM), by considering the non-isothermal processes.( 2) Presenting a mathematical model and numerical simulation of the heat exchange between the fluid in the crack and the surrounding intact porous ambient with diffusive fracture edges.(3) Studying the effect of the temperature contrast between the injection fluid and the residing porous material on the pressure needed to initiate the crack and on postcrack pressure response.(4) Discussing and overcoming the stability challenges in the convection-dominated heat transport state.To give an overview of mathematical modeling, Sect. 2 briefly describes the continuum mechanical treatment, the kinematics, balance, and constitutive relations.It also demonstrates the phase-field approach and the essential energy functions.In the end, this section summarizes the governing system of coupled nonlinear partial differential equations (PDEs), which are needed to solve IBVPs of non-isothermal phase-field hydraulic fracturing in saturated porous media.The numerical treatment of the resulting set of PDEs using the FEM is briefly presented in Sect.3.This includes defining the Dirichlet and Neumann boundary conditions, explaining the staggered-monolithic solution scheme, and addressing the numerical challenges and the instability sources of the coupled problem.The verification of the proposed modeling framework and proof of model applicability are carried out in Sect. 4 through solving a number of two-dimensional numerical IBVPs and comparison with reference solutions from the literature.Finally, a short summary, conclusions, and an outlook of future works are given in Sect. 5.

Homogenization and densities
In the underlying work, we consider a two-phase immiscible porous material consisting of a materially compressible pore-fluid phase u F and a materially incompressible solid-phase u S .The compressibility of the solid material is ignored in comparison with that of the solid matrix.In this, bulk compressibility is possible due to the presence of compressible pore fluid and under drained conditions, where the fluid can escape through the boundaries.Within the TPM-based continuum mechanical treatment, homogenization is applied to the multiphase material so that each spatial point consists of overlapping and interacting solid and fluid constituents, i.e., u ¼ u S S u F .The volume fraction n a is defined as the ratio of the volume element dv a to the total volume element dv , i.e., n a :¼ dv a =dv with a ¼ S for the solid phase and a ¼ F for the fluid phase.Following this, the saturated condition is expressed as: The continuum mechanical TPM applied to intact saturated porous materials demands that all constituents exist at each material point of the homogenized domain, i.e., 0\fn S ; n F g\1, see, e.g., [28].In the case of a fractured saturated porous material, extended phase-field-dependent porosity and solidity formulations need to be applied so that 0 n S \1 and 0\n F 1. Thus, for a complete crack region that is filled with the fluid, n F ¼ 1 and n S ¼ 0 as discussed in [52].In particular, the deterioration of n S is connected with the crack evolution, which is realized in our treatment using the phase-field method (PFM) as will be discussed in later sections.We assume in this study that the material is found under non-isothermal conditions, so that the material densities can be expressed as functions of the current temperature h a of the solid or fluid phases.For the materially compressible fluid constituent, the fluid material density is additionally expressed as a function of the pore pressure p as will be discussed in Sect.2.6.1.In particular, the material densities of u S and u F can be expressed as Herein, q aR is the current material density of the constituent u a , q SR 0 is the initial solid density, a a h is the thermal expansion coefficient of u a , and h a 0 is the initial temperature.For more details, we refer to, e.g., [57,95] and the references therein.To complete the discussion, the initial partial density q a 0 of a constituent u a and following its motion can be expressed as with n a 0S being the initial volume fraction of u a .In particular, n S 0S is the undeformed solid volume fraction at t ¼ t 0 and n F 0 S ¼ 1 À n S 0 S follows the saturation condition.

Kinematics
In describing the motion of multiphasic continua within the TPM and presenting the kinematic formulations, we start by defining the position vector x of all constituents at the current configuration and X a as the position of the phase u a at the reference configuration, see, e.g.[28].For each constituent, a Lagrangian motion function v a and a velocity vector v a are defined as In this, v À1 a represents the inverse of the motion function, which allows for an Eulerian (spatial) description of the motion.Following this, the material time derivative of an arbitrary vector-valued quantity ðÞ and a scalar-valued quantity ðÞ with respect to the motion of constituent u a can be expressed, respectively, as ðÞ 0 a :¼ The motion of the solid phase u S is described via the solid displacement u S and velocity v S , while the motion of the fluid phase is described either via Eulerian description using the fluid velocity v F or via modified Eulerian settings using the seepage velocity w F : For the case of small deformations, the strain tensor e S can be expressed as This includes an elastic contribution e e S and a thermal contribution e H S :¼ a S h ðh S À h S 0 Þ I with a S h as the thermal expansion coefficient of u S .
The formulation of the PFM relies on the definition of the phenomenological scalar-valued phase-field variable d S as: In this, d S depends on the position vector x and time t and allows to approximate the sharp edges of the crack by diffusive ones.More details can also be found in, e.g., [30,31,50,52,60,61,72,76]. Taking into account C as the set of sharp cracks that are found in a homogenized domain X and dA as the crack surface unit, the total crack surface A C as a volume integral rather than a surface integral can be expressed as In this, C d is the crack surface density function per unit volume, which is formulated in terms of d S , its gradient r d S , and the diffusive crack regularization parameter l c [ 0, see, e.g., [12,71] for details about the origin of Eq. ( 9) and [30,52,72] for applications within hydraulic fracturing.

Energy formulations
Proceeding with G c ðh S Þ as the temperature-dependent Griffith's critical fracture energy release rate along with C d defined in Eq. ( 9) 2 , the fracture energy can be expressed as Presenting the definition of W S frac helps in defining the overall energy function in addition to the derivation of the phase-field evolution equation as will be discussed in Sect.2.6.5.Moreover, a proper relationship for G c ðh S Þ, such as according to [25], is presented in Sect.2.6.5.
For a saturated porous material with a linear elastic solid phase and subjected to a thermal field, the total free energy density W is defined within the small strain assumption as the sum of the free energy densities of the solid and fluid constituents, i.e.
In this, W S el , W S th; el are, respectively, the elastic stored energy and the energy term due to the thermoelastic coupling of the fractured solid.The dependency on d S in these terms is connected to the energy dissipation or solid stiffness degradation, as the fracture occurs and propagates in the solid phase.The term W S th in (11) is the purely solid thermal energy.For the fluid phase, W F el is the elastic stored energy of the fluid phase presented as a function of the material density q FR ðp; h F Þ, see, e.g., [29].However, alternative formulations to this can be found within Biot's theory, see, e.g., [69,92].Moreover, W F th is the purely fluid thermal energy.The current energy functions of the fluid phase assume that the thermoelastic coupling for the fluid phase is negligible.Moreover, in the numerical examples in this contribution, we consider nearly incompressible fluids, such as water.In this, the dependency of the fluid material density on the pressure p can be neglected.Thus, for simplicity, the dependency of W F el on p in the following treatment will be neglected.
Proceeding with the method that was introduced by [45] to model crack onset in brittle solid materials and following the scheme explained in, e.g., [49], the total potential energy F of a cracked fully saturated and isotropic linear elastic porous solid can be formulated as To capture the fracture-related deterioration of the material stiffness, a degradation function, such as ð1 À d S Þ 2 , can be applied in the total energy formulation, see, e.g., [58,71].
In this, it is common to use a residual stiffness factor 0\g ( 1 in the degradation function that enters the momentum balance to prevent the elements on the crack path from having zero stiffness and, thus, avoid numerical stability problems.In particular, the degradation function becomes In this regard, it is assumed that the cracking occurs under tension and shear stresses, but not under compression.Thus, the related free energy can additively be split into a positive part ð. ..Þ þ , associated with tension and shear, and a negative part ð. ..Þ À , associated with compression.Therefore, the total local energy W can be expressed in an abstract form as More details about the energy formulation for the considered material will be presented in the following subsections in connection with constitutive modeling.

Model assumptions
In preparing the governing balance relations in the underlying study, we apply a number of assumptions, which are summarized in the following: • Non-polar solid constituent, which implies a symmetric stress tensor.• The deformation of the solid skeleton takes place under quasi-static conditions and within the small strain framework.• Under non-isothermal conditions, some of the material parameters can be formulated as functions of the temperature.These include the Young's modulus of solid matrix E S ¼ E S ðh S Þ, the critical solid tensile strength r S crit ¼ r S crit ðh S Þ, the effective fluid dynamic viscosity l FR ¼ l FR ðh F Þ, and the material densities.

Constituent balance relations of saturated porous media
Under the aforementioned assumptions within the TPM, the constituent balance equations of the two-phase saturated porous media can be expressed as • Partial mass balance: • Partial energy balance: Herein, qa is the mass supply or production term, which fulfills qS þ qF ¼ 0 for a closed system and vanishes if no mass exchange between the constituents exists.r a is the symmetric partial Cauchy stress tensor for the nonpolar solid material and b is the mass-specific body force.Additionally, pa is the direct momentum production, which represents the volume-specific local interaction force between the fluid and the solid phases and fulfills the relation pS þ pF ¼ 0. In the energy balance Eq. ( 16), e a is the mass-specific internal energy, L a is the spatial velocity gradient, q a is the heat influx vector, r a is the mass-specific external heat supply, and êa represents the constituent energy production.In (16), the constituent Helmholtz free energy w a can be incorporated using Legendre transformation between the entropy g a and the temperature h a as conjugate variables.In particular, using in ( 16) allows to have the free energy instead of the internal energy in the constituent energy balance equations.In addition to Eqs. ( 14)-( 16), the entropy inequality is usually used to ensure that the constitutive and evolution relation are thermodynamically admissible, which is beyond the scope of this work.

Effective stresses and permeability formulation
Following the concept of the effective stresses for multiphase porous media, the total solid and fluid stresses and the direct momentum production term are split into an 'extra' or 'effective' part, indicated by ðÞ E and a weighted pore-pressure part (p), see, e.g., [23].In the presence of thermal fields, the constituent temperatures also enter the effective stress formulations as has been discussed in, e.g., [27,95].In particular, these formulations are derived in a thermodynamically consistent framework based on the evaluation of the entropy inequality.Thus, under nonisothermal conditions and having a biphasic fluid-saturated porous material, the following relations hold: In the current work, we proceed from a fractured linear elastic solid phase, where the effective stress tensor r S E is derived based on the solid free energy function q S 0 w S ¼ W S , which includes also the pure solid thermal energy W S th ¼ q S 0 w S th .In particular, W S ðe S ; h S ; d S Þ can be written as Herein, hÁi AE ¼ ðÁ AE j Á jÞ=2 represents the Macaulay brackets, tr ðe S Þ :¼ e S Á I and e D S :¼ e S À 1 3 tr ðe S Þ I are the trace and the deviator of the strain tensor e S , respectively.Additionally, j S is the bulk modulus of the porous solid matrix, l S is the shear modulus, and a S h is the thermal expansion coefficient of the solid phase.Moreover, Dh S : ¼ h S À h S 0 with h S 0 as the initial or reference temperature of the solid phase and C S V is the specific heat capacity of the solid at constant volume.With the definition of the initial density in (3), i.e., q S 0 ¼ n S 0S q SR 0 , the stored thermal energy of the solid phase will vary in connection with the solid volume fraction variation due to the crack onset or in connection with predefined cracks as will be discussed in Sect.2.6.1.The dependency of W S on d S is due to the ad hoc added degradation function Gðd S Þ, where d S will not be considered when applying the time derivative of W S .Moreover, the material parameters of the solid matrix might also be considered as functions of the temperature in this non-isothermal treatment, as will be discussed in the numerical examples.
Following this, the solid effective stress tensor can be derived as For the fluid effective stress tensor formulation, we proceed from the assumptions of Newtonian fluid and negligible average normal viscous stress (i.e., Stokes's hypothesis), see, [52].Additionally, as discussed in Sect.2.4 and in [92], we neglect for the fluid phase the thermoelastic coupling.Thus, the fluid effective stress tensor under nonisothermal conditions can be expressed as with l F :¼ n F l FR ðh F Þ being the fluid dynamic viscosity and l FR ðh F Þ [ 0 is the effective fluid viscosity.A widely applied empirical formulation for the temperature-dependent dynamic fluid viscosity is that proposed by [15], which can be expressed as In this, a l and b l are model parameters, which can be adjusted based on available experimental data of the considered fluid.
Proceeding from an intact porous material (d S ¼ 0) with interacting solid and fluid constituents and following the discussion presented in [52], the effective momentum production term pF E can be expressed as with K poro being the diagonal entry of the isotropic intrinsic permeability tensor K S poro :¼ K poro I.In this, pF E represents a local interaction force and takes into account the mechanisms of momentum transfer between the deforming solid and the permeating viscous fluid.Thus, it includes information about the fluid viscosity and the deformation-dependent intrinsic permeability.This term enters the fluid momentum balance, which results in Darcy's-like filter law.
The momentum production term is additionally multiplied by a phase-field-based degradation function ð1 À d S Þ to capture the gradual decrease of the interaction forces pF E from the intact region to the crack region, i.e., pF E ! 0 as d S ! 1 .
Apart from fractured porous media, the elastic deformations of the solid skeleton might lead to changing of the intrinsic permeability as has been discussed in, e.g., [64].This effect can be captured by introducing an adjusting function c u , which is related to the permeability as In this, K S 0 is the initial intrinsic permeability, and j !0 is a model parameter, which can be adjusted based on reference or experimental data.Additionally, n F 0S and n S 0S are the initial porosity and solidity, respectively.In dealing with brittle geomaterials, such as rocks, this effect can be neglected as fracturing might take place for small deformations.In the context of fluid flow in porous media and in the absence of fractures, it was also shown in, e.g., [24,34] that the fluid effective stress is much smaller than the momentum production pF This argument, discussed in detail in [24,34], is based on a micro-to macroscale dimensional analysis, which includes a comparison between the macroscopic Darcy's law that results from the fluid momentum balance and the microscopic Poiseuille-type fluid flow that includes information on the microtopology (like pore diameters and numbers) together with the pore-fluid viscosity.It is shown under the assumption of a large number of pore channels with small diameters that the fluid effective (friction) stress can be neglected, whereas it is sufficient to consider the effect of fluid viscosity in the fluid momentum balance through the effective momentum production term.The onset of cracking in the solid matrix leads to anisotropy in the effective medium.This is associated with anisotropic flow, i.e., enhanced flow in the crack region, as thoroughly discussed in [51,52,72].To describe the occurring anisotropic flow within the framework of porous media mechanics, a permeability measure that realizes the crack opening and direction needs to be derived.In particular, for the flow in the crack region with pF E ! 0 as d S ! 1 , the fluid momentum balance equation (Eq.( 15) with a ¼ F) together with the definition of fluid effective stress (18) 2 and n F ¼ 1 yields the Stokes' free flow equation, expressed as The flow in the crack within the assumption of laminar flow can approximately be considered to be similar to the flow between two parallel plates with a distance w c as illustrated in Fig. 1.Through integration in the x 1 -direction devoted to Eq. ( 26) and considering the simplified boundary conditions results in the volume flux per unit area of the crack cross-section q F1 as For an arbitrary crack orientation and using the similarity to Darcy's law of porous media, Eq. ( 27) can be defined in tensorial form as In this, K c is defined as the fracture intrinsic permeability tensor.This can be expressed according to [51] as and In this, w h represents the hydraulic crack width, which can be found in closed or opened states.For this, w c represents the opened-crack aperture, w r is the residual closed-crack aperture, f c is a correction factor that captures the roughness of the fracture surfaces [54,88,102], and v d is a levelset function that determines the cracked region using a Heaviside step function H as v d ¼ Hðd S À 0:5Þ.In Eq. ( 29), the orientation of the crack is determined using the unit vector n d , which is perpendicular to the fracture edge and defined via the gradient of the phase-field variable.In the current treatment, w c is estimated at each quadrature point of the FE model based on the stretching of a 1D line element h e normal to the fracture path.Thus, having the strain tensor e S , this can be computed as Fig. 1 Illustration of the unidirectional Poiseuille-type flow between two plate with a distance w c and depth L c in x 2 direction, where v F1 ðx 3 AE w c =2Þ % 0, see, [52] With the definitions of the permeability in ( 24) and ( 29), the following unified intrinsic permeability formulation can be proposed and covers the whole fractured domain: where b ¼ 1 for a linear permeability transition and b ¼ 2 for a quadratic transition.An important outcome of the aforementioned discussion is the ability to describe the fluid flow in the fractured domain using a unified Darcy's law.Thus, proceeding with the fluid partial momentum balance (15) for u F and considering the relations ( 18), ( 23), (25), and (31), the unified Darcy's law can be expressed as In this, the presence of the crack is associated with n F ! 1, w F ! v F , and q F ! q FR .More details about the derivation of Eq. ( 32) can be found in [52].

Heat flux, energy production, and internal energy
In the underlying work, a local thermal non-equilibrium model is considered, which allows capturing the heat transfer between the solid and the fluid phases.Therefore, two distinct energy balance equations for u S and u F with interphase heat transfer laws need to be considered.
Recalling the energy balance in local form (16), we discuss in the following the related constitutive equations.Proceeding with isotropic thermal conduction for u a of the biphasic material, the heat influx q a can be specified according to Fourier's law of heat conduction as q S :¼ Àn S H SR grad h S ; where H SR and H FR are the real or intrinsic heat conductivity of constituents.The variation of n S and n F in accordance with the change of d S , i.e., n S !0 and n F ! 1 for d S ! 1, allows for a continuous variation of the heat influx and temperatures across the crack diffusive interface.It is worth mentioning here that other approaches exist in the literature to treat the heat transfer at the crack-ambient interface as discussed in, e.g., [1,79].For instance, it is shown in [1] within the sharp crack interface approach and under different conditions for the fluid flow that the continuity in the heat flux in the model leads to an accurate approximation of the temperature change between the porous medium and the fluid layer that represents the crack in our treatment.In other works, like [93] within a pure solid fracture, a degradation function is multiplied by the thermal conductivity to capture the crack effect on the heat transfer.
In analogy to the effective stress formulations presented in Eq. ( 18), the evaluation of the entropy inequality for the underlying non-isothermal process results in a relationship between the total entropy (g a ) and the effective entropy g a E of the constituents.Following the derivations in [57,95,96], the effective and total entropy can be expressed for the solid phase, respectively, as In the aforementioned relations, the definition of the solid Helmholtz free energy w S in Eq. ( 19) has been used, whereas m h :¼ À 3 j S a S h is the stress-temperature coupling modulus.The neglection of the mechanical contribution in the definition of g S E is in analogy to [57].This analogy is also extended to the formulation of the fluid entropy in the following.Thus, we start with the definition of the fluid thermal energy w F th as with C F V being the specific heat capacity of fluid at constant volume and q F 0S is the initial fluid partial density.Following this, the effective and total entropy can be expressed for the fluid phase, respectively, as Within the assumption of non-equilibrium thermal condition, the heat exchange through the interface between the solid and fluid phases is realized through the direct energy production terms êF .Following [7,29,86,92,94], this can be expressed based on Newton's law of cooling as where x d is the specific heat transfer surface area and k I h is the heat transfer coefficient thoroughly discussed in [41].
The second term in Eq. ( 37) results from the total energy production term under assumptions like no mass exchange between the constituents as derived in detail in, e.g., [94].

While the parameter k I
h is approximately considered constant in this study and independent of the temperature or flow state, i.e., the flow takes place at a low Reynolds number, x d is introduced as a function of d S .In particular, the onset of the crack is associated with distortions at the crack edges as can be observed, for instance, in l-CT images of fractured Sandstone in [56].These distortions lead to increasing the heat transfer surface area between the solid and the fluid phases.This was also experimentally proved in [41], which showed that the crack surface roughness has a major effect on the heat transfer coefficient.Assuming that the crack is free of the solid phase, then x d is considered to be zero in the completely cracked region.Thus, utilizing the Heaviside level-set function v d ¼ Hðd S À 0:5Þ, the phase-field-dependent formulation of x d can be expressed as In this, x 0 is the specific heat transfer surface area for the intact material, and x I captures its increase in the diffusive interface area.An alternative formulation can be found [92], which assumes that fragmented solid particles are found in the crack, and thus, x d increases in the crack region.This effect can be captured using Eq. ( 38) after removing ð1 À v d Þ.

Fluid mass balance
In the simulation of hydraulic fracturing, a fluid mass supply term qExt in Eq. ( 14) as an external injection source can be utilized, see, [52].For this, a location function QðxÞ is specified to determine the location the fluid injection supply source, i.e., QðxÞ ¼ 1 for the source location and QðxÞ ¼ 0 otherwise.Hence, the supply term can be expressed as qF :¼ qExt QðxÞ.Proceeding from Eq. ( 14) for u F , the fluid mass balance can be rewritten as In the reformulation of Eq. ( 39), we consider the following material time derivatives Moreover, a proper formulation for the term ðn S Þ 0 S can be derived based on Eq. ( 14) for u S and having no solid mass production, i.e., qS ¼ 0. In particular, we start with the material time derivative Integration of ðn S Þ 0 S and reformulation of the resulting equation within geometrically linearized settings according to [29,95] yield an explicit relation for n S as with n S 0S being the initial solidity.With the assumption that the preexisting or emerging cracks are free of the solid phase, a modified initial solidity n S d as a function of d S can be considered in Eq. ( 42) in analogy to [49,51].In particular, the reformulated solidity can be expressed as Moreover, under the assumption of a barotropic fluid phase, a simplified relationship among the fluid material density q FR ðp; h F Þ, the fluid compressibility j F , and the effective pore pressure p can be considered in analogy to [11].Additionally, the dependency on the pressure field h F as shown in Eq. ( 2) has to be considered in the material derivative of q FR .Having this, one gets Having this relation together with (40) and the solidity derivation in (41), the governing fluid mass balance can be rewritten as In this, ð1 À d S Þ can additionally be multiplied by div v S to eliminate this term in the fracture region.Moreover, to reduce the number of primary variables, Eq. ( 32) can be incorporated in Eq. ( 45) so that no need to explicitly compute n F w F in the simulation.

Overall momentum balance
The derivation of the overall momentum balance proceeds from the constituent momentum balance (15) for u S , u F , and the effective stress formulations (18) as The definition of the solid effective stress r S E ¼ r S E ðe S ; h S ; d S Þ is presented in (20), whereas the argument presented in ( 25) is also implemented.Following this, the overall momentum balance related to the quasi-static treatment can be expressed as

Energy balance of u F
Starting with the local energy balance ( 16) for u F and utilizing the Legendre transformation in (17) to define the fluid internal energy, i.e., e In this, L F :¼ grad v F is the fluid velocity gradient.The different terms in Eq. ( 48) have to be specified in order to get the governing equation that is used in the finite element treatment.In particular, we derive the material time derivative of the fluid thermal energy w F th in (35) as well as of the fluid entropy in (36) 2 yielding We also reformulate the term r F Á L F based on the fluid mass balance (39) for Q ¼ 0 together with ( 40), (41), and (44) yielding In this, we also neglect the effective part of the total fluid stress formulation in analogy to [57], and we apply the assumption of grad ðn S Þ % 0 for the geometrically linear framework [46].With ( 36), ( 49), ( 50), the fluid energy balance (48) can be reformulated as In connection with mass supply term qF ¼ qExt QðxÞ presented in Eq. ( 39), injecting fluid with a different temperature from the domain needs to be considered in the modeling framework.Therefore, an external energy supply term ŵF has additionally been considered in the fluid energy balance equation.In particular, this can be expressed as with h F s being the temperature of the supplied fluid.The relation (52) relies on the transferred heat per unit time under steady-state condition as discussed, e.g., by [41].

Energy balance of u S
Starting with the local energy balance (16) for u S and utilizing the Legendre transformation in (17) to define the solid internal energy, i.e., e S ¼ w S þ h S g S , yields In this, L S ¼ grad v S is the solid velocity gradient, which can additively be split into a symmetric rate of deformation Having a symmetric stress tensor r S yields r S Á L S ¼ r S Á D S ¼ r S Á ðe S Þ 0 S , where r S Á W S ¼ 0 and the definition of the strain tensor is given in Eq. ( 7).As we did for the fluid energy balance, the different terms in Eq. ( 53) are to be specified for the numerical finite element implementation.We carry out first the material time derivative of the solid Helmholtz free energy w S in (19) as well as of the fluid entropy in (34) 2 yielding Having these derivations together with the total solid stress (18) 1 , the solid entropy (34) 2 , and the heat exchange term in (37) 2 , the solid energy balance can be expressed as

Phase-field evolution equation
Several approaches can be found in the literature to derive the governing equation of the phase-field variable evolution within hydraulic fracturing [48].In the current work, we follow the generalized Ginzburg-Landau approach in analogy to [67,72], which relies on Allen-Cahn equation [5].Proceeding from the local definition of the total strain energy (13) together with the fracture energy (10) and the specification of the solid energy in (19), the Allen-Cahn equation can be expressed as with M [ 0 being a scalar-valued interface mobility parameter.Considering a quasi-static formulation to the phase-field evolution as in [52], i.e ð1=MÞ ðd S Þ 0 S % 0, the phase-field equation ( 56) can be reformulated as where Herein, W Sþ max is the continuously increasing elastic strain energy density, presented to prevent crack healing upon possible crack closure [70].Within the non-isothermal treatment, a proper relationship for G c as a function of the solid temperature h S needs to be specified.For instance, [25] proposed the following empirical relation: where G c0 represents the critical energy release rate at the reference temperature h S 0 and a c [ 0 is a model parameter, which can be adjusted based on the considered material.
The formulation for the crack deriving force in (58) has the disadvantage of underestimating the solid crack resistance as discussed in [72].Therefore, an alternative threshold-based formulation for the crack driving force can be utilized.Specifically, a threshold-based crack driving force Dc is proposed in [72], where this force remains equal to zero as long as the positive thermoelastic energy or the related principal stresses r Sþ E; pi with i 2 f1; 2; 3g smaller than a certain critical value r S crit .Following this, the crack driving force can be expressed as In analogy to Eq. ( 59), where G c is expressed as a function of h S within a non-isothermal framework, a proper formulation for r S crit that considers the dependency on h S will be discussed in the following.In this, the one-dimensional solution of the phase-field evolution equation under quasistatic condition, as, e.g., in [59], shows a proportionality between the critical stress needed to initiate the crack r and G c as r / ffiffiffiffiffi ffi G c p .Based on this, we propose the following empirical relation for r S crit under non-isothermal conditions: where r S crit;0 represents the critical stress needed to initiate the crack at the reference temperature h S 0 and z c [ 0 is a model parameter, which can be adjusted based on experimental data of the considered material.

Summary of the governing balance relations
Following the previous discussion, the governing system of coupled nonlinear partial differential equations (PDEs), which are needed to solve IBVPs of non-isothermal phasefield hydraulic fracturing in saturated porous media, are summarized in Box 1.

Numerical treatment
In the current work, we use the finite element method (FEM) in the numerical solution of initial-boundary value problems (IBVP) of hydraulic fracturing under nonisothermal conditions.In the following, we briefly discuss the main steps in the numerical treatment of the governing coupled partial differential equations (PDEs), summarized in ( 62)- (66)  with n being an outward-oriented unit surface normal.In deriving the governing weak formulations, the method of weighted residuals is applied.In this, time-independent weight (or test) functions, as variations of the real primary fields, are specified, i.e., dnðxÞ :¼ ½du S dp dh F dh S dd S T .These test functions are then multiplied by the strong form Eqs. ( 62)-( 66) and integrated over the spatial domain.Moreover, the boundary terms are derived by utilizing the product rules and the Gauss integral theorem to the integrands.For simplicity, the variational formulations of the governing equations will be referred to as For the time discretization, the system of coupled equations is solved in the FE package FlexPDE 7.16 [39], which uses the second-order implicit backward difference formula (BDF2), discussed in, e.g., [42], for the time integration.In our implementation of the first IBVP in Sect.4.1, a three-step staggered procedure is applied, in which the phase field, the hydro-mechanical, and then the thermal problems are sequentially solved.
For the spatial discretization within the FEM and following the Bubnov-Galerkin approach, the primary variables nðx; tÞ and the weighting functions dnðxÞ are approximated over the discrete spatial domain using proper shape functions and nodal values, i.e.,

Weak form equation Abbreviation Test function
Overall momentum balance In the implementation in FlexPDE, an adaptive timestep strategy is considered so that the root-mean-square (RMS) errors remain below an assigned limit, whereas the initial time step is set to Dt ¼ 1 Â 10 À4 s.

Overcoming numerical stability challenges
In the FE implementation, equal-order, quadratic shape functions are applied for all variables.Hence, in solving the coupled mechanical problem in a monolithic manner (step 3 in Box 2), one has to overcome a number of stability challenges in the limits of very small permeability in the porous domain, high fluid viscosity, near crack-domain interfaces, and next to drained boundaries with p ¼ 0. In particular, the existence and uniqueness of the solutions require the fulfillment of the inf-sup condition as discussed in, e.g., [66,80].Therefore, within the approach of mixed quasi-compressibility, stability can be enforced by applying a pressure projection stabilization procedure.Specifically, a stabilization term of the form rÁ ða st r _ pÞ is added to the fluid mass balances, where a st [ 0 is a small stabilization parameter, chosen as small as possible to maintain the stability but not deteriorate the accuracy of the computation.It is to note here that the stabilization term within the FE approach is only applied in the element interiors and not on the boundary elements.
Another source of instability arises from the fluid energy balance (step 4 in Box 2) for the case of convectiondominated heat transport, i.e., with high fluid velocity and heat transported mainly by the flow of the fluid, as explained in, e.g., [57].In the current treatment, the fluid velocity increases significantly in the crack region in comparison with the surrounding ambient, which presents a source of instability in the numerical treatment.To understand this challenge, the fluid energy balance can be rewritten as a simple convection-diffusion equation as where the scalar-valued terms A; B; D can be determined based on comparing equation ( 69) with (64).The stability of the numerical simulation of such problems using the FEM with the Bubnov-Galerkin approach is sensitive to both the mesh size and the fluid velocity.The stability challenge, which is particular for the convective-dominated process, is characterized by a high value of the Fig. 2 The geometry and boundary conditions considered in the numerical modeling of heat transmission in the fractured porous domain problem.The preexisting cracked zone is colored blue and has a width of 0.28 m.For the boundary conditions, the upper and lower boundaries are fixed in the vertical direction (u S2 ¼ 0), the lateral boundaries are fixed in the horizontal direction (u S1 ¼ 0), zero fluid flux (q F ¼ 0) is applied to all boundaries except for the inlet and the outlet of the crack, zero pressure (p ¼ 0) at the outlet, and zero heat flux (q h F ¼ 0; q h S ¼ 0) is applied at all boundaries except for the inlet of the crack dimensionless Pe ´clet number P e .This number can be computed in the one-dimensional stationary convectiondiffusion process with an element size of h e as For P e [ 1, the convective part in the fluid energy balance (69) becomes dominated and a suitable stabilization technique needs to be applied.If P e \1, then the diffusive part becomes dominated and no special stabilization strategies are required [44].
In the underlying work, a stabilization technique based on adding an isotropic artificial thermal diffusivity term D in the domain is followed.By doing so, the fluid energy balance (69) can be expressed as In the literature, several approaches can be found for the determination of the stabilization term D in convectiondiffusion problems, see, e.g., [26].Many of the approaches rely on the Pe ´clet number P e to assign a suitable value for D as can be found in [57].For simplicity in the current treatment, the value of D has been chosen as small as possible and the same in the domain to get a stable solution.

Numerical examples
As numerical applications of the discussed formulations and algorithms, two numerical examples are introduced in the following and implemented in the FE package FlexPDE 7.16 [39].In the first example, the focus is on the modeling and numerical challenges of the heat transmission process in a fractured medium with reference to an analytical solution, whereas the second example investigates the effect of temperature-dependent material properties on the cracking onset and propagation as well as the effect on the fluid flow in the crack.

Heat transmission in fractured medium with reference to an analytical solution
An analytical solution that estimates the heat transmission in a wellbore system is developed by [87].This model considers the fluid flow through a pipe in a wellbore and the heat transmission between an injected cold fluid and an ambient with a higher temperature.The output of this model is the time-dependent temperature change of the outlet fluid, which also depends on the depth of the outlet  point, i.e., h F ¼ h F ðx; tÞ.In this study, we rely on the analytical solution in [87], as provided in a MATLAB code in [20], in order to evaluate the performance of our nonisothermal phase-field saturated porous media model in capturing the heat transfer while the fluid flows in a fractured medium.The geometry and boundary conditions for the finite element solution of the plane-strain two-dimensional initial-boundary value problem (IBVP) are illustrated in Fig. 2. In this, we consider a porous domain (gray) with an initial temperature of 55 C and a preexisting cracked zone along the porous domain with an initial fluid temperature of h F0 ¼ 20 C. A fluid is injected into the inlet point with a flow rate of q F = 0.00013 m 3 /s and heat flux of q h F ¼ q F q FR 0S C F V ðh F À h F0 Þ, which reduces the fluid temperature at the inlet to 20 C.
For simplicity, the underlying treatment neglects the thermal expansion effects, i.e., we apply a S h ¼ a F h ¼ 0 in the set of governing equations in Box (1).The applied material parameters for the hydro-thermal processes in this example are summarized in Table 1 and based on several sources, such as [20,52,57,96].
The preexisting cracked region is generated via applying an initial high value of the crack driving force, i.e., Dc ) 0, in the crack domain, which allows realizing the diffusive crack edges as illustrated in Fig. 3.In this, the internal length scale is chosen l c ¼ 0:14 m in Eq. (66).
Focusing on the fluid temperature at the outlet, which is affected directly by the heat exchange between the injection fluid and the surrounding ambient, we study in this example the following factors: between two plates, while the analytical solution assumes an axisymmetric geometry.Another reason is the diffuse nature of the crack edges in the model, while they are sharp in the analytical solution.
In the FE solution, the fluid temperature at the outlet (initially 20 C) is affected by the initial temperature of the ambient (55 C) and, thus, increases to a value of approximately 43 C before starting to decrease under the effect of the injected cold fluid.In the analytical solution, no initial temperature is applied (see, e.g., [20]), which explains the differences between the two solutions at the beginning.The effect of the initial temperature decreases over time and the two solutions become comparable with a slight deviation afterward.The decrease in the fluid outlet temperature is accompanied by decreasing in the ambient temperature due to heat flow from the surrounding hot regions to the colder fracture region.This is illustrated in Fig. 5, where the solid temperature at a point B (1 m above the fractured zone) decreases gradually, while the fluid outlet temperature keeps decreasing due to lower temperature gradients with the surrounding.
In connection with the fluid temperature change and heat exchange, Fig. 6 shows an example of the spatial and temporal change of h F at four subsequent time steps.It is clear the temperature changes much faster in the crack as the permeability is much higher than the surrounding domain.
In the numerical examples, an artificial thermal diffusivity term D is assigned to guarantee the numerical stability in the convection-dominated heat transport.This is particularly important in the case of high fluid velocity, where heat is transported mainly by the flow of the fluid as explained in Sect.3.1.The approach followed in the numerical examples is to apply as small values of D as possible so that the numerical stability is recovered while the accuracy is minimally impaired.In this connection,    In the underlying treatment with a non-equilibrium thermal condition (cf.Sect.2.5.4), the heat exchange through the interface between the solid and fluid phases is realized through the direct energy production terms êF , which is proportional to the specific heat transfer surface area x d and the heat transfer coefficient k I h .Having a higher value of the heat transfer coefficient allows for good heat exchange between the solid and the fluid phases and, thus, the fluid temperature at the outlet decreases much slower than the cases with lower k I h as illustrated in Fig 8 .In real fracture problems, e.g., in rocks, the presence of fragmented solid particles in the crack and the absence of flat crack edges cause the approximation of the flow in the crack as flow between two plates and the calculation of permeability based on the cubic law (29) to be inaccurate, i.e., it leads to an overestimation of the crack permeability.Therefore, we test the effect of the permeability in the crack on the outlet fluid temperature.The results presented in Fig. 9 show that the lower permeability in the crack allows for a slower flow and, thus, better heat exchange with the surroundings than the case of higher crack permeability.

Hydraulic fracturing in saturated granite rocks with temperature-dependent parameters
In the following, a 2D plane strain numerical example of phase-field hydraulic fracturing saturated porous media using the FE package FlexPDE 7.16 is presented.By considering temperature-dependent material parameters, the aim is to test the influence of initial domain temperature and injection-fluid temperature on the fracturing process, such as the pressure needed to initiate the crack.While the hydraulic fracturing process happens within a relatively short time, the time needed to have a significant change in the surrounding temperature and, thus, change in the values of temperature-dependent material parameters, as we saw in the previous example, is much longer.Therefore, for simplicity, we consider in this example merely the initial temperatures and the related values of material parameters of the solid and fluid constituents at that particular temperature.In analogy to the example presented in [52], we consider in the following hydraulic fracturing of granite rocks with water as the injection fluid.For the temperature-dependent material parameters, the focus is laid on Young's modulus of the solid phase E S , the tensile strength r c , and the dynamic fluid viscosity l FR .The choice of these values is based on several references, such as [35] and [89].In particular, the deterioration of the parameter values of granite in the range of 100-400 °C and that of water in the range of 20-100 °C is summarized in Table 2 for the four cases of study.
The other material parameters of the considered granite rock samples with presumably linear elastic material response until the onset of the fracture are summarized in Table 3, see [36,52] for more details.
The geometry and the boundary conditions of the considered problem are illustrated in Fig 10 .In this, a rectangular domain of dimensions 450 Â 300 mm 2 and a notch at the center are considered.The fluid injection is applied at the center of the domain with a constant injected flow rate of q F = 0.0635 mm 3 /s until the onset and propagation of the crack.For the boundary conditions, the pressure confining stresses r 1 and r 2 are applied to the lateral and upper  In the finite element implementation, triangular elements with quadratic shape functions are used.The mesh is refined in advance at the expected crack path boundaries, respectively.Additionally, the boundaries of the sample are considered drained, i.e., p 0 ¼ 0. The other boundaries are chosen in a way that no rigid-body motion occurs.
Under the assumption of constant temperature during the cracking process, the governing balance equations with primary variables u S ; p; d S are the overall momentum balance (62), the fluid mass balance (63), and the phasefield evolution Eq. ( 66).The finite element solution of the governing equations is accomplished in the FE package FlexPDE 7.16 with triangular quadratic elements, whereas the second-order implicit Backward Difference Formula (BDF2) is applied for the stepping.As a solution strategy, a staggered scheme is applied, where at each time step the phase-field equation is solved first and the coupled porous media equations are solved next.The stability of the solution is enhanced with the application of a pressure projection stabilization procedure discussed in Sect.3.1.In this, a stabilization term rÁ ða st r _ pÞ with a st ¼ 0:01 is added to the fluid mass balances.
As exemplary results of the finite element simulations, Fig. 11 presents the contour plots of the phase-field variable d S , the fluid pressure p, and the vertical displacement u y for case (4) of Table 2.In this, the crack propagates as expected normal to the maximum confining stress r 1 .The fluid pressure has its maximum value in the crack and in the nearby area, which corresponds to a clear high-pressure gradient in that region.Moreover, the crack opening is associated with a displacement jump as can be interpreted from u y in Fig. 11.
In discussing the effect of initial temperatures on the cracking process, Fig. 12 shows the time histories of the pressure at the injection source for the four cases of temperature-dependent material parameters presented in Table 2.
Based on these results, it is clear that the solid phase (granite) temperature governs the maximum pressure needed to initiate the crack.In this, for cases (1) and ( 2) with a lower solid temperature and, thus, higher values of Young's modulus E S and the tensile strength r c , higher values of injection pressure are needed to start the cracking.The effect of fluid phase (water) temperature is obvious in the crack propagation phase (post-crack initiation phase).In this, the higher temperature [cases (2) and ( 4)] is associated with lower fluid viscosity and, therefore, one has a more toughness-dominated hydraulic fracturing response with more physical oscillations in the fluid pressure during crack propagation, see, e.g., [48].The oscillations in the fluid pressure field during crack propagation prove the capability of the model in capturing the stepwise crack advancement.With a lower fluid temperature [cases (1) and (3)] and, thus, higher fluid viscosity, a more stable crack propagation (viscosity-dominated) can be observed, see, e.g., [52] for experimental evidence.In this regard, the topic of stepwise crack propagation in porous media and the related fluid pressure oscillations has been studied numerically and in many research works as reviewed in, e.g.[18].In this, experiments from the petroleum industry [10] proved the stepwise behavior of hydraulic fracture, where these results are reported in [38].
The stepwise crack propagation and the related pressure oscillation are affected by several factors, including the mechanical boundary conditions, the injection rate, the fluid viscosity, the heterogeneity the material, and preexisting cracks.Studying the different factors is, however, beyond the scope of the underlying work.

Conclusions and future work
Within hydraulic fracturing under non-isothermal conditions, this work presented a continuum mechanical modeling framework that utilizes the theory of porous media (TPM) together with the phase-field method (PFM).In this, a two-phase porous material consisting of thermally nonequilibrated and mechanically interacting solid-fluid constituents is considered.The continuum mechanical approach with assumptions of small deformations and nomass exchange between the constituents resulted in the most general case in a coupled five partial differential equations.These are the overall momentum balance, the fluid mass balance, the fluid energy balance, the solid energy balance, and the phase-field evolution equation.The numerical implementation of the strongly coupled equations succeeded within a finite element framework, where a number of numerical stability challenges had to be overcome.This includes stabilizing the solution algorithm in the limit of very small permeability or materially incompressible fluid by adding a Laplace pressure term to the fluid mass balance equation within a forcing quasi-compressibility approach.The other stability challenge that had to be overcome is related to the fluid energy balance equation in the convection-dominated heat transport mechanism with high fluid velocity and heat transport mainly by the flow of the fluid.In this, a stabilization technique based on adding an isotropic artificial thermal conductivity term to the fluid energy balance equation was used.
To show the capability of the proposed approach and the stability and robustness of the applied numerical algorithm, two initial-boundary-value problems were presented.In the first one, the simulation of heat transfer in a fractured medium was compared with a reference analytical solution.In addition to demonstrating good agreement with the reference analytical solution, this example discussed the various factors that affect output accuracy, such as permeability in the crack, the conductivity stabilization parameter, and the heat transfer coefficient.In the second numerical example, the influence of temperature-dependent material properties on crack initiation and propagation as well as the influence of temperature on fluid flow in the crack were discussed.This example showed that rocks with lower temperatures require higher pressure for crack initiation and propagation.Moreover, a higher temperature of the fluid lead to a more toughness-dominated hydraulic fracturing response, whereas a viscosity-dominated response is obtained for lower-temperature fluids.
To conclude, the novelty of this contribution lies in deepening the understanding of non-isothermal hydraulic fracturing and its associated coupled processes.It addressed the effects of temperature variations on the fluid and rock during the fracturing process and provided methods to account for these effects and interactions in numerical simulations.This work is of critical importance to industries that use non-isothermal hydraulic fracturing, such as enhanced oil and gas recovery, geothermal power generation, and underground carbon dioxide storage.In our future work on modeling hydraulic fracturing under non-isothermal conditions, a number of important issues have been identified to be addressed.This includes validation of the model by comparing the numerical results with experimental data from laboratories or from in-field enhanced geothermal systems.In this, it is also important to carry out efficient simulations within three-dimensional configurations in order to get more accurate numerical results.Another topic to be addressed is related to extending the model to capture the possible water-vapor phase change under high temperatures in the cracking process.Finally, the issue of convection-dominated heat transport in case of turbulent flow will be studied in future works.In this, the heat transfer coefficients have to be carefully adjusted as they can be much greater than that of laminar flow and the Pe ´clet number is anyhow much greater.
Laminar fluid flow in the crack and in the intact porous regions.• The saturated porous material includes incompressible solid and compressible fluid.• No mass production from the interaction between the solid and the fluid phases, however, a fluid source term qF is considered in the mass balance equation to represent the external fluid supply.• The model is under local thermal non-equilibrium conditions, where the different temperature variables are chosen for the solid and the fluid phases.• The natural convection due to buoyant forces is neglected.• The fracturing takes place under quasi-static conditions.

Fig. 3 Fig. 4
Fig. 3 Contour plot of the phase-field variable d S , which illustrates the cracked region in the FE model

Fig. 5 Fig. 6
Fig. 5 The solid temperature history at point A (20,3.14) in the domain and the fluid temperature history at point B (20,2.14) at the outlet for the case of k I h = 10 W/(m 2 K), K S c = 10 À8 m 2 , D = 200

10 Fig. 7
Fig. 7 The effect of the stabilization parameter D on the solution of the fluid temperature at the outlet for the case of k I h = 100 W/(m 2 K) and K S c = 10 À8 m 2

Fig. 8 −9 m 2 Fig. 9
Fig.8The effect of the heat transfer coefficient k I h on the fluid temperature at the outlet for the case of K S c = 10 À8 m 2 and D = 10

Fig. 7
Fig.7illustrates the effect of D on the accuracy, whereas choosing higher values of stabilization parameter than needed might lead to over-estimate or under-estimate the temperature of the fluid in the fractured region.In the underlying treatment with a non-equilibrium thermal condition (cf.Sect.2.5.4), the heat exchange through the interface between the solid and fluid phases is realized through the direct energy production terms êF , which is proportional to the specific heat transfer surface area x d and the heat transfer coefficient k I h .Having a higher value of the heat transfer coefficient allows for good heat exchange between the solid and the fluid phases and, thus, the fluid temperature at the outlet decreases much slower than the cases with lower k I h as illustrated in Fig 8.In real fracture problems, e.g., in rocks, the presence of fragmented solid particles in the crack and the absence of flat crack edges cause the approximation of the flow in the crack as flow between two plates and the calculation of permeability based on the cubic law (29) to be inaccurate, i.e., it leads to an overestimation of the crack permeability.Therefore, we test the effect of the permeability in the crack on the outlet fluid temperature.The results presented in Fig.9show that the lower permeability in the crack allows for a slower flow and, thus, better heat exchange with the surroundings than the case of higher crack permeability.

Fig. 10
Fig. 10 Illustration of the geometry and boundary conditions of the two-dimensional model with a rectangular notch of 18 Â 1:75 mm 2 .In the finite element implementation, triangular elements with quadratic shape functions are used.The mesh is refined in advance at the expected crack path

Fig. 11 FluidFig. 12
Fig.11Exemplary contour plots of the phase-field variable d S , the fluid pressure p, and the vertical displacement u y at t ¼ 4000 s for case(4) and characterized by five primary variables nðx; tÞ :¼ ½u S p h F h S d S T .We start by discussing the boundary conditions in connection with each of the variables.As such, denoting B as the domain of interest with oB as its boundary, this can be split for each variable into Dirichlet (oB D ) and Neumann (oB N ) boundaries that fulfill the conditions oB D [ oB N ¼ oB and oB D \ oB N ¼ ;.In particular, the Dirichlet and Neumann boundary conditions are given, receptively, for each variable as the discretized trial spaces and dn % dn h over the discretized spaces of weighting functions.
S as the space-discrete phase-field equation, the sequential steps at a time t ðkþ1Þ ¼ t ðkÞ þ Dt with Dt as the time step size are summarized in Box 2.

Table 2
[35,89]r cases of temperature-dependent material parameters with Young's modulus E S , the tensile strength r c , and the dynamic fluid viscosity l FR[35,89]