Sanghyun Lee

Andro Mikeli Ć †

Mary F Wheeler

AND Thomas Wick

PHASE-FIELD MODELING OF TWO PHASE FLUID FILLED FRACTURES IN A POROELASTIC MEDIUM

Keywords: Finite elements, Phase-field fracture, Biot system, Fixed-stress iterative coupling, Multiphase flow, Fracture propagation

We propose an immiscible two phase flow fracture model, based on phase-field for treating crack propagation in porous media. This multifluid model is an extension of classical flow models and we take into account non-zero capillary pressure. Using lubrication theory, we provide details of the determination of effective parameters: absolute and relative permeabilities. The phase-field formulation is a generalization of previous works by the authors and extends the single phase model to the two phase case. Here the resulting flow system has four unknowns: resident and injected pressures and saturations, respectively. The solid contribution consists of displacements and a phasefield variable. Both systems are coupled employing a fixed-stress splitting. Therein, the flow problem is treated with an iterative scheme and the solid problem fully implicitly. Modeling and algorithms are substantiated with several numerical tests.

1. Introduction. Coupled multiphase fluid flow and deformation in fractured porous rock are important in energy and environmental applications. Hydraulic fracturing and horizontal drilling are considered to be the major reasons U.S. has been having an energy revolution, one that has changed the energy picture from scarcity to abundance. Fracking, the forcing opening of fissures in subsurface rocks by introducing liquid at high pressure, is letting the U.S. tap vast oil and natural gas reserves that previously were locked away in shale and other tight-rock formations. Up to 95 percent of natural gas wells drilled in the next decade will utilize hydraulic fracturing [START_REF] Energyfromshale | What is Fracking?[END_REF] and by 2040, over 80 percent of US natural gas production is projected to come from formations that need to be hydraulically fractured to be accessible [1]. Additional important examples involving subsurface fluid structure interactions include geothermal energy production, deep underground injection of hazardous wastes, and storage of carbon sequestration into saline aquifers and depleted oil and gas reservoirs. The latter has been considered to be a major technology for reducing the emission of greenhouse gases to the atmosphere. Numerical simulations are essential in understanding the performance of the above coupled systems in a highly complex heterogeneous porous media with multiple temporal and spatial scales. For example modeling hydraulic fracturing requires three dimensional compositional fluid models, multiphase, thermal, proppant transport with gravitational settling, water solute transport and includes quasi-Newtonian fluid rheology with modified power law coupled with a geomechanics elasticity model that can treat fracture propagation. These simulator capabilities enable comparison of proposed fracturing designs for predicted production and refracturing near depleted wells and optimization designs for slick water injection and possible usage of other fracturing fluids such as CO 2 and foam. Strategies for carbon sequestration and geothermal production utilize similar mathematical models and numerical simulators. These numerical studies are important in addressing environmentally prudent procedures.

There is a long list of environmental issues that need to be considered such as groundwater contamination and surface water pollution. In addition to water quality issues, fracking wells release compounds into the air, such as benzene, ethylbenzene, toluene and n-hexane; longterm exposure to these has been linked to birth defects, neurological problems, and blood disorders. In addition many areas of the United States, not considered earthquake-prone, such as Texas and Oklahoma, are now experiencing relatively strong seismic activity. In 2015, an area of North Texas had nine confirmed earthquakes in a 24-hour period. Fracking is being investigated as the cause; hence, a quantitative understanding of fracking is crucial.

There has been considerable research in modeling fluid flow in a fractured porous media going back to the 1960s with the development of the PKN model [START_REF] Nordgren | Propagation of a vertical hydraulic fracture[END_REF][START_REF] Perkins | Widths of hydraulic fractures[END_REF] and the KGD model [START_REF] Zheltov | On hydraulic fracturing of an oil-bearing stratum[END_REF]. These two dimensional models assumed fractures with constant height and differed primarily because of the assumption of plane strain in the vertical and horizontal directions, respectively. In recent years, a common practice has been the linking of simulators, each specialized to flow or mechanics, as in the case of TOUGH2 and FLAC3D [START_REF] Rutqvist | A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, and deformation in fractured porous rock[END_REF]. For more details regarding the coupling between fracture geometry and fluid flow problems, the reader is referred to [START_REF] Adachi | Computer simulation of hydraulic fractures[END_REF][START_REF] Economides | Reservoir stimulation[END_REF]. Recently, several other approaches have been discussed in [27, 30, 37-39, 48, 59, 71, 73, 78, 86].

The authors have been involved in the development of high fidelity phase-field models for fracture propagation in a poroelastic medium. Recent advances and numerical studies for treating multiphysics phase-field fractures include the following; thermal shocks and thermoelastic-plastic solids [START_REF] Bourdin | Morphogenesis and propagation of complex cracks induced by thermal shocks[END_REF][START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. Part ii. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic plastic solids[END_REF][START_REF] Miehe | Phase field modeling of fracture in multiphysics problems. Part i. Balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids[END_REF], elastic gelatin for wing crack formation [START_REF] Lee | Investigation of wing crack formation with a combined phase-field and experimental approach[END_REF], pressurized fractures [START_REF] Bourdin | A variational approach to the numerical simulation of hydraulic fracturing[END_REF][START_REF] Mikelić | Phase-field modeling of pressurized fractures in a poroelastic medium[END_REF][START_REF] Mikelić | A quasi-static phase-field approach to pressurized fractures[END_REF][START_REF] Wheeler | An augmented-Lagangrian method for the phase-field approach for pressurized fractures[END_REF][START_REF] Wick | Coupling fluid-structure interaction with phase-field fracture[END_REF][START_REF] Wick | 3D phase-field for pressurized fracture propagation in heterogeneous media[END_REF], fluid filled (i.e., hydraulic) fractures [START_REF] Heider | A phase-field modeling approach of hydraulic fracture in saturated porous media[END_REF][START_REF] Lee | Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model[END_REF][START_REF] Lee | Optimal design of hydraulic fracturing in porous media using the phase field fracture model coupled with genetic algorithm[END_REF][START_REF] Markert | Coupled multi-field continuum methods for porous media fracture[END_REF][START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part iii. crack driving forces in hydro-poro-elasticity and hydraulic fracturing of fluidsaturated porous media[END_REF][START_REF] Miehe | Minimization principles for the coupled problem of Darcy Biot-type fluid transport in porous media linked to phase field modeling of fracture[END_REF][START_REF] Mikelić | A phase-field method for propagating fluidfilled fractures coupled to a surrounding porous medium[END_REF][START_REF] Wick | Fluid-filled fracture propagation using a phasefield approach and coupling to a reservoir simulator[END_REF], proppant-filled fractures [START_REF] Lee | Phase-field modeling of proppantfilled fractures in a poroelastic medium[END_REF], variably saturated porous media [START_REF] Cajuhi | Phase-field modeling of fracture in variably saturated porous media[END_REF], and crack initiations with microseismic probability maps [START_REF] Lee | Initialization of phase-field fracture propagation in porous media using probability maps of fracture networks[END_REF]. These computational examples demonstrate the potential of phase-field for fracture propagation.

Phase-field is a regularized approach that clearly has advantages and shortcomings. A list from the authors' experience is as follows. The first advantage is a continuum description based on first physical principles to determine the unknown crack path [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and the computation of curvilinear and complex crack patterns. The model allows for nucleation, branching and merging and postprocessing of certain quantities such as stress intensity factors becomes redundant. Therefore, easy handling of fracture networks in possibly highly heterogeneous media can be treated. The formulation being described in a variational framework allows finite element discretizations and corresponding analyses. Phase-field is a fixed-mesh approach in which no remeshing or update of basis functions to resolve the crack path is needed. The mathematical model permits any dimension d = 2, 3, and in case the software allows as well, phase-field fracture applies conveniently to three-dimensional simulations. On the energy level, the formulation is non-convex constituting a challenge for both theory and design of numerical algorithms. A second challenge is the computational cost since either additional iterations in an alternating approach [START_REF] Bourdin | Numerical implementation of the variational formulation for quasi-static brittle fracture[END_REF] are required or a fully nonlinear problem has to be solved. However, our quasi-monolithic approach [START_REF] Heister | A primal-dual active set method and predictorcorrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF][START_REF] Lee | Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model[END_REF] has low Newton iteration numbers and therein the major cost goes into solving the vector-valued elasticity problem that also applies using most other crack models (such as extended/generalized finite elements for instance). Two bigger downsides are the accurate crack width computations, despite that some recent ideas have been proposed [START_REF] Lee | Iterative coupling of flow, geomechanics and adaptive phase-field fracture including level-set crack width approaches[END_REF][START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF], and also the smeared fracture transition zone when additional physics (a given pressure, proppant, or saturations as in the current paper) shall be described in, or around, the fracture or on the fracture boundary. Here, local mesh adaptivity helps to a great extend in order to achieve sufficient accuracy. As demonstrated in our own work and by others (see the references in the previous paragraph), ideas have been proposed how to cope with these challenges. In view of the increasing popularity, it however seems that the advantages outweight the shortcomings for today's applications.

In this paper, we are specifically interested in extending our porous media phase-field fracture model to two phase flow inside the fracture. While two phase flow for non-fractured media is well-established [START_REF] Aziz | Petroleum reservoir simulation[END_REF][START_REF] Chavent | Mathematical models and finite elements for reservoir simulation: single phase, multiphase and multicomponent flows through porous media[END_REF][START_REF] Peaceman | Fundamentals of Numerical Reservoir Simulation[END_REF], the extension to fractures that are propagating is a challenge. The novelties in the development of our two phase flow model for a propagating fracture is based on two steps. First, using lubrication theory, we compute the absolute and relative permeabilities. In the second step, these quantities are substituted into the full three dimensional equations for two phase flow.

With regard to numerical modeling, we finally need to solve four problems: a nonlinear displacement/phase-field system, a pressure system, a saturation equation, and crack width problem. Here, the two important aspects in question are the coupling algorithm and discretization schemes with appropriate conservation properties. For the displacement/phasefield fracture subproblem, we employ a quasi-monolithic approach [START_REF] Heister | A primal-dual active set method and predictorcorrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF][START_REF] Lee | Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model[END_REF][START_REF] Wick | An error-oriented Newton/inexact augmented Lagrangian approach for fully monolithic phase-field fracture propagation[END_REF]. To couple this fractured-geomechanics problem to flow, we use a splitting method. The numerical discretization in space is based on continuous Galerkin finite elements for displacements and phase-field, and an enriched Galerkin finite element approach [START_REF] Lee | Adaptive enriched Galerkin methods for miscible displacement problems with entropy residual stabilization[END_REF][START_REF] Lee | A locally conservative enriched Galerkin approximation and efficient solver for elliptic and parabolic problems[END_REF][START_REF] Sun | A locally conservative finite element method based on piecewise constant enrichment of the continuous galerkin method[END_REF] is applied for both the pressure and saturation equations in order to obtain local mass conservation with an efficient solver.

The outline of this paper is as follows: In Section 2, we briefly recapitulate phase-field modeling of hydraulic fractures in porous media and a width computation. Section 3 presents the derivation of a two phase flow model for propagating fractures; we discuss the single phase flow model for a reservoir in Section 4. In Section 5, we present the global solution algorithm and explain variational formulations and finite element discretization including decomposition of the fracture and reservoir domain as a diffraction system. In the short Section 6, we briefly describe the numerical solvers and emphasize the importance of local mesh refinement. Several numerical examples are provided in Section 7. A summary of our results is given in the conclusions in Section 8.

2.

Modeling phase-field fracture in a porous medium. We now recapitulate the essentials for phase-field modeling for pressurized and fluid filled fractures in porous media, as described in [START_REF] Lee | Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model[END_REF][START_REF] Mikelić | A quasi-static phase-field approach to pressurized fractures[END_REF][START_REF] Mikelić | Phase-field modeling of a fluid-driven fracture in a poroelastic medium[END_REF].

We seek a vector-valued displacement u(•, t) : Λ × (0; T ) → R d (d ∈ {2, 3}), a scalar-valued pressure p(•, t) : Λ×(0; T ) → R and a scalar-valued phase-field function (namely a smoothed indicator function) ϕ(•, t) : Λ × (0; T ) → [0, 1], which will satisfy the equations of the model in the weak sense. In order to satisfy the initial conditions, the above functions are considered as Hilbert space valued functions, defined on [0, T ]. They take values in appropriate Sobolev spaces with respect to the spatial variable. 2.1. Constitutive laws in the reservoir. The solid displacements and diffusive flow in a non-fractured porous medium are modeled in Ω R by the classical quasi-static ellipticparabolic Biot system for a porous solid saturated with a slightly compressible viscous fluid. Let p R = p| Ω R be the pressure in the reservoir and p F = p| Ω F be the pressure in the fracture.

The constitutive equation for the Cauchy stress tensor is given as

σ por (u, p R ) -σ 0 = σ R (u) -α(p R -p 0 )I in Ω R × (0, T ), (2.1)
where u is the displacement, p 0 is the initial reservoir pressure, I is the identity tensor, σ 0 is the initial stress value and α ∈ [0, 1] Biot's coefficient. The effective linear elastic stress tensor is

σ R := σ| Ω R = σ R (u) := λ(∇ • u)I + 2Ge(u), (2.2)
where λ, G > 0 are the Lamé coefficients. The linear elastic strain tensor is given as e(u) := 1 2 (∇u + ∇u T ). Then the balance of linear momentum in the solid reads

-∇ • σ por (u, p R ) = ρ s g in Ω R × (0, T ), (2.3)
where ρ s is the density of the reservoir solid and g is the gravity. For simplicity in this paper, we assume homogeneous Dirichlet boundary conditions on ∂Λ for the displacement u.

2.2.

A quasi-static phase-field model for pressurized fractures. As in the variational setting for linear elasticity in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], our phase-field formulation utilizes the elliptic (Ambrosio-Tortorelli) functional [START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functionals via γ-convergence[END_REF][START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF]. In the following, we recapitulate our phase-field fracture approach for pressurized and fluid filled fractures in porous media as established in [START_REF] Mikelić | Phase-field modeling of pressurized fractures in a poroelastic medium[END_REF][START_REF] Mikelić | A phase-field method for propagating fluidfilled fractures coupled to a surrounding porous medium[END_REF][START_REF] Mikelić | A quasi-static phase-field approach to pressurized fractures[END_REF].

Here, we assume that the fracture is a three dimensional thin object where the width is much smaller than the length, and where the leading order of the fracture fluid stress is derived from lubrication theory.

The effective fluid stress is given by σ F := σ| Ω F = -p F I. As transmission conditions, we assume continuity of normal stresses:

σ por • n F = σ F • n F on Γ F × (0, T ), or (σ R (u) -αp R I) • n F = -p F n F on Γ F × (0, T ), (2.4)
where Γ F := Ω R ∩ Ω F is the fracture interface (crack surfaace) and n F is the outward pointing normal vector on Γ F . In the following we use the second transmission condition; namely, continuity of pressures

p R = p F on Γ F × (0, T ). (2.5)
Consequently we identify p R = p F = p and obtain:

(σ R (u) -αpI) • n F = -pn F on Γ F × (0, T ). (2.6)
Next, we follow [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF] in which the stress tensor σ is additively decomposed as

σ = σ R = σ + (u) + σ -(u)
into a tensile part σ + (u) and a compressive part σ -(u) by:

σ + (u) := ( 2 d G + λ)tr + (e(u))I + 2G(e(u) - 1 d tr(e(u))I), (2.7)
σ -(u) := ( 2 d G + λ)tr -(e(u))I, (2.8)
and tr + (e(u)) = max(tr(e(u)), 0), tr -(e(u)) = tr(e(u)) -tr + (e(u)).

(2.9)

We notice that however the influence of different energy/stress splitting laws is still a subject of discussion in the fracture mechanics community (see [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF] and Section 2.2 in [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF]). Moreover, pressurized fractures are basically subject to tension forces in which the splitting into σ + and σ -becomes redundant; see Section 5.2 in [START_REF] Mikelić | A quasi-static phase-field approach to pressurized fractures[END_REF].

By applying (2.6) and Gauss' divergence theorem (see Section 3.2 in [START_REF] Mikelić | Phase-field modeling of a fluid-driven fracture in a poroelastic medium[END_REF]), the strong Euler-Lagrange Formulation 1 for a quasi-static model of a pressurized phase-field fracture is described below. By utilizing the definition of ϕ and the interface conditions; namely (2.5) and (2.6), we do not distinguish between σ R , p R and σ F , p F in the global dissipation form.

FORMULATION 1 (Pressurized displacement phase-field fracture). Let the phase-field regularization parameters are ε > 0 and κ > 0 with κ ε. Here ε is the length of a smooth phase-field transition zone where ϕ ∈ (0, 1); we refer to Figure 5.2 for details. Given a pressure p and a phase-field ϕ, find displacements u such that

-∇ • (1 -κ)ϕ 2 + κ σ + -∇ • σ -+ ϕ 2 ∇p + (α -1)∇ ϕ 2 p = 0 in (0, T ) × Λ, u = 0 on (0, T ) × ∂ D Λ.
Similarly, given displacements u and given a pressure p, find a phase-field variable ϕ such that

-G c ε∆ϕ - G c ε (1 -ϕ) + (1 -κ)σ + : e(u)ϕ +2(1 -α)ϕp(∇ • u) + 2ϕ ∇p • u ∂ t ϕ = 0 in (0, T ) × Λ, ∂ϕ ∂n = 0 on (0, T ) × ∂Λ, ϕ(0) = ϕ 0 in {t = 0} × Λ,
where n is the outward pointing normal vector on ∂Λ and G c > 0 is the energy release rate.

In addition, the phase-field variable is subject to an irreversibility constraint and yields an inequality system including a compatibility condition

-G c ε∆ϕ - G c ε (1 -ϕ) + (1 -κ)σ + : e(u)ϕ + 2(1 -α)ϕp∇ • u + 2ϕ∇p • u ≤ 0 in (0, T ) × Λ, ∂ t ϕ ≤ 0 in (0, T ) × Λ,
which is frequently referred as a Rice condition. Moreover, in the discretized system, we also require h K ε, where h K is the local discretization parameter.

REMARK 2.1. Because the displacement equations are quasi-static, without explicit timedepenence, the displacement value at t = 0 can be computed with ϕ 0 and p| t=0 .

2.3. Computation of the fracture width. In the displacement phase-field system, the width of the fracture is not explicitly defined. For accurate width computations of non-planar fractures an additional algorithm is required. Formally the definition of the width reads:

2w := 2w(u) = 2h = -[u • n F ],
which is the jump of the normal displacements, where n F is the unit normal on the crack surface (Γ F ). For the meaning of h, we refer the reader to Figure 3.1.

Using explicitly the phase-field variable, this definition can be rewritten as:

2w = -[u • n F ] = -2u • ∇ϕ LS ∇ϕ LS ,
by approximating the normal vector on the crack surface (Γ F ) with the level set values ϕ LS of a specific isoline of the fracture. To compute ϕ LS , we proceed as follows: Let Γ F be the fracture interface. We now define Γ F as the zero level-set of a function ϕ LS such that

ϕ LS > 0, x ∈ Ω R , ϕ LS < 0, x ∈ Ω F , ϕ LS = 0, x ∈ Γ F ,
where the isoline ϕ LS is obtained from the phase-field variable ϕ by

ϕ LS := ϕ -C LS .
(2.10)

Here, C LS ∈ (0, 1) is a constant selected for approximating the fracture boundary Γ F . Throughout this paper for simplicity, we set C LS = 0.1. Thus,

Γ F := {x ∈ Λ | ϕ(x, t) = C LS }, Ω R := {x ∈ Λ | ϕ(x, t) > C LS } and Ω F := {x ∈ Λ | ϕ(x, t) < C LS }.
Details of the entire procedure and underlying motivation are provided in [START_REF] Lee | Iterative coupling of flow, geomechanics and adaptive phase-field fracture including level-set crack width approaches[END_REF][START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF].

REMARK 2.2. In order to improve the fracture width representation inside the fracture, we finally solve a crack width interpolation problem; see [START_REF] Lee | Iterative coupling of flow, geomechanics and adaptive phase-field fracture including level-set crack width approaches[END_REF] (therein Formulation 6 and Figure 4). We omit the continuous problem statement here, but we directly provide the Galerkin approximation in Formulation 6.

3. Modeling two phase flow in a fracture. We concentrate now on a detailed derivation of a two phase flow model inside a propagating fracture. A characteristic situation is illustrated in Figure 3.1. An important aspect is to formulate a novel model with the surrounding reservoir. The goal is to derive a two dimensional pressure-saturation system from which we extract absolute and relative permeabilities and capillary pressure. These quantities are then used in a standard three-dimensional two phase flow system. To this end, we develop a lubrication model that not only computes p F but also the saturation and provides an explicit computation of the relative permeabilities. Our model is derived from first principle laws and results via upscaling into a three dimensional lubrication model. Finally, we use a physically based phenomenological constitutive law to account for the capillary pressure.

3.1.

A model for a two phase flow in a flat symmetric 3D penny shape fracture. In the following, we derive a model for two phase flow in a fracture that includes explicitly computed absolute and relative permeabilities. To this end, we start from first principles by working with the multifluid Stokes equations. For the fracture geometry we assume a flat three dimensional fracture occupying the domain Ω F and which is characterized by the ratio ε between its characteristic width H in the direction x 3 and its characteristic horizontal length L in directions x 1 and x 2 ; a sketch is provided in Figure 3.1.

3.1.1. Coupling conditions and a two phase flow Stokes model. Our first assumption is that the crack boundary propagates in a deformable porous medium. This boundary can be described as x 3 = h(x 1 , x 2 , t) for x 3 ≥ 0 and x 3 = -h(x 1 , x 2 , t) for x 3 < 0.

(3.1) Fig. 3.1: Sketch of the near-tip fracture region for a two dimensional setting. Fluid res is the residing fluid in the domain and Fluid inj is injected from a source in the fracture. p F,res and p F,inj are the corresponding pressure values. The interface A is defined between the two fluids.

Fracture boundary x 3 = ±h(x 1 , x 2 , t) moves in time. We assume that the characteristic width H of the fracture is much bigger than the pore size of the porous medium. We recall that for H = c 0 l 2/3 where l is the pore size, a coupled model must be taken into account.

In three dimensions, the x 2 coordinate goes into the plane and the one-dimensional curve h(x 1 , x 2 , t) becomes a two dimensional surface. We finally notice that H = c 0 l (2/3) is a kind of threshold value. If H is smaller than that value than the effective filtration does not see the fracture. If H is bigger than the threshold value then one has to couple two adjacent domains; see [START_REF] Bourgeat | Effective fluid flow in a porous medium containing a thin fissure[END_REF].

The so-called tip region is defined where h and -h meet at x 3 = 0. We further assume that the fluid flow in the fracture is given by the incompressible Navier-Stokes equations for two immiscible viscous fluids (see e.g., [START_REF] Temam | Navier-Stokes Equations: Theory and Numerical Analysis[END_REF]). The second fluid is being injected through an injection well in the middle of the crack and pushes the residing fluid. Between these two immiscible fluids, we introduce a free boundary x 3 = A(x 1 , x 2 , t) (i.e., the interface). On the interface A we deal with the usual transmission conditions:

• continuity of the velocities;

• continuity of normal stresses. REMARK 3.1 (Surface tension). Initially, we neglect the surface tension. These effects will be taken into account later via energy interface considerations.

Since we assume two immiscible fluids, we have a third condition on the interface A. Let v be the fluid velocity, then:

DA Dt -v 3 | x3=A = ∂A ∂t -v 3 | x3=A + (v 1 , v 2 )| x3=A • ∇ x1,x2 A = 0. (3.2)
To derive a two phase flow model, classically one would distinguish the two flow equations. We utilize the work of Demay, Nouri and Poupaud in [START_REF] Nouri | An existence theorem for the multi-fluid Navier-Stokes problem[END_REF][START_REF] Nouri | An existence theorem for the multi-fluid stokes problem[END_REF], who formulated a model, which assumes that the transmission and non-miscibility conditions at the interface follow from the global formulation of a linear transport equation for the viscosity. Furthermore, the momentum and mass conservation conditions are described by the standard (linear) Stokes equations. Consequently, our two fluid system has the form: Find v, p, η such that

∇p -∇ • 2ηe(v) = 0 in Ω F × (0, T ), (3.3) ∇ • v = q in Ω F × (0, T ), (3.4) ∂η ∂t + ∇ • ηv = qη in Ω F × (0, T ), (3.5
)

η| t=0 = η 0 (x) in Ω F × (0, T ), (3.6)
where q is a source term, v is the velocity, p is the fluid pressure, e(v) is the rate-of-strain tensor and η is the viscosity. Specifically, the viscosity can only take two values, i.e., η ∈ {η F,res , η F,inj }, where η F,res and η F,inj are the residing and injecting viscosities, respectively. We notice that the flow equations are given in the entire domain Ω F and to determine the specific phase is achieved with the help of (3.6).

A rescaled

Stokes system in a penny shaped geometry. Next, we derive a dimensionless form of the lubrication equation corresponding to a displacement of the fracture fluid by a radial source injection following the classical lubrication theory references (see e.g. [START_REF] Szeri | Fluid film lubrication[END_REF]). For the channel and tube flow such equations were derived in [START_REF] Fourar | A viscous coupling model for relative permeabilities in fractures spe 49006[END_REF][START_REF] Mikelić | On the derivation of the Buckley-Leverett model from the two fluid Navier-Stokes equations in a thin domain computational geosciences[END_REF]. Here we consider a penny-shaped fracture. We will follow the two scaled expansions introduced in [START_REF] Mikelić | On the derivation of the Buckley-Leverett model from the two fluid Navier-Stokes equations in a thin domain computational geosciences[END_REF].

Let U be the characteristic velocity in the directions x 1 and x 2 and Re= U L/η F,res is the Reynolds number. We assume that Re is small, which justifies using the Stokes equation ((3.3) and (3.4)) instead of the full Navier-Stokes equations.

Next, we set the characteristic quantities:

t = U L t; x = x 1 L , y = x 2 L , z = x 3 H , v x = v 1 U , v y = v 2 U , v z = v 3 L U H , a = A H , p = H 2 η F,res LU p, w(x, y, t) = 1 H h(x 1 , x 2 , t), q = qL U , ε = H L , ω = Ω F L , η ε = η η F,res , M = η F,inj η F,res .
We remark that the ε in this section is not linked to the phase-field regularization parameter.

For simplicity, we skip the bars in the following and denote the rescaled

Ω F = ω × (-h, h)
with the same symbol. In addition, we add the index ε to the unknowns, to indicate that the system under consideration contains a small parameter.

With these preliminaries, the system (3.3)-(3.6) can be re-written as

-∂ x p ε + ∂ z η ε ∂ z v ε x + O(ε 2 ) = 0 in Ω F , (3.7)
-∂ y p ε + ∂ z η ε ∂ z v ε y + O(ε 2 ) = 0 in Ω F , (3.8)
-∂ z p ε + O(ε 2 ) = 0 in Ω F , (3.9) ∂ x v ε x + ∂ y v ε y + ∂ z v ε z = q in Ω F , (3.10)
∂ t η ε + ∂ x η ε v ε x ) + ∂ y η ε v ε y ) + ∂ z η ε v ε z ) = qη ε in Ω F . (3.11)
Our main assumption is that Fluid res is already present in the fracture (and also in the reservoir) and that the injected fluid does not reach the fracture boundary:

η ε | t=0 = M χ {-a0(x1,x2)≤z≤a0(x1,x2)} + χ {-w≤z<-a0(x1,x2) and a0(x1,x2)<z≤w} , supp q ⊂ {-a 0 (x 1 , x 2 ) ≤ z ≤ a 0 (x 1 , x 2 )}, (3.12)
where a 0 ≥ 0 is a smooth surface with compact support in Ω F .

3.1.3. Normal and tangential vectors. We recall that the unit normal n and tangential vectors τ (j) , j = 1, 2 can be expressed in the form

n = -ε∂ x we 1 -ε∂ y we 2 + e 3 1 + ε 2 |∇ x,y w| 2 , τ (1) = e 1 + ε∂ x we 3 1 + ε 2 |∂ x w| 2 , τ (2) = e 2 + ε∂ y we 3 1 + ε 2 |∂ y w| 2 . Then for v = U (v ε x , v ε y , v ε z ) we have v • n = -U ε∂ x wv ε x -U ε∂ y wv ε y + U εv ε z 1 + ε 2 |∇ x,y w| 2 = U ε(v ε z -∇ x,y w • (v ε x , v ε y )) + O(ε 3 ); v • τ (1) = U v ε x + U ε∂ x wv ε z 1 + ε 2 |∂ x w| 2 = U v ε x + U ε(∂ x wv ε z -ε v ε x 2 (∂ x w) 2 ) + O(ε 3 ),
and analogously for v • τ (2) . The velocity in the fracture is larger than in the surrounding porous medium. Hence we have on the boundary z = w(x, y, t):

0 = v • τ (1) = U v ε x + U ε ∂ x wv ε z -ε v ε x 2 (∂ x w) 2 + O(ε 3 ) (3.13)
and an analogous expression for v • τ (2) .

3.1.4.

A kinematic condition for the fracture boundary. For later purposes, we recall the kinematic boundary condition for the fracture boundary (cf. [START_REF] Mikelić | A phase-field method for propagating fluidfilled fractures coupled to a surrounding porous medium[END_REF]). In the normal direction we impose

Dh Dt -v out = ∂ t h -v out + (v 1 , v 2 ) • ∇ x1,x2 h -v 3 = 0 on x 3 = h(x 1 , x 2 , t). (3.14)
The term v out ensures mass conservation inside the fracture and is well-known as the leak-off q L . In the dimensionless form, the condition reads:

∂ t w + (v ε x , v ε y ) • ∇ x,y w -v ε z = -v LO on z = w(x, y, t), (3.15)
where v LO = v out U ε .

3.1.5.

A second re-scaled system by using the lubrication assumptions. As in lubrication theory, we can now expand the velocities, pressure and viscosities as follows

v ε x = v 0 x + ε 2 v 1 x + O(ε 4 ); v ε y = v 0 y + ε 2 v 1 y + O(ε 4 ); v ε z = v 0 z + ε 2 v 1 z + O(ε 4 ); (3.16) p ε = p 0 + ε 2 p 1 + O(ε 4 ), (3.17)
η ε = η 0 + ε 2 η 1 + O(ε 4 ). (3.18)
Here, the upper indices 0 and 1 denote the respective order of the underlying expansions. We recall that the fracture boundary given by (3.1) is given in dimensionless form by: +w(x, y, t) for z ≥ 0, -w(x, y, t) for z < 0.

(

-∂ x p 0 + ∂ z η 0 ∂ z v 0 x = 0, -∂ y p 0 + ∂ z η 0 ∂ z v 0 y = 0, -∂ z p 0 = 0, (3.20)
∂ x v 0 x + ∂ y v 0 y + ∂ z v 0 z = q, (3.21)
∂ t η 0 + ∂ x η 0 v 0 x ) + ∂ y η 0 v 0 y ) + ∂ z η 0 v 0 z ) = qη 0 . (3.22)
Specifically, on the fracture boundary w, we have

v 0 x (x, y, w, t) = v 0 y (x, y, w, t) = 0, (3.23)
which can be seen with the help of (3.13) using the leading order term.

3.1.6. Derivation of a pressure lubrication equation. Following [START_REF] Szeri | Fluid film lubrication[END_REF], we use (3.20)-(3.23) to calculate v 0

x , v 0 y and obtain a PDE for p 0 . We first compute the permeability, which depends non-locally on the viscosity η:

K(η 0 , x, y, z, t) = z -w ξ + w η 0 (t, x, y, ξ) dξ - w -w ξ + w η 0 (t, x, y, ξ) dξ z -w 1 η 0 (t, x, y, ξ) dξ w -w 1 η 0 (t, x, y, ξ) dξ , (3.24)
v 0 x = K∂ x p 0 , v 0 y = K∂ y p 0 , ∂ t w = v 0 z (x, y, w) -v LO . (3.25)
Next we average the continuity equation (3.21) with respect to z, which yields

∂ x w -w v 0 x dz + ∂ y w -w v 0 y dz + v 0 z (x, y, w(x, y, t), t) -v 0 z (x, y, -w(x, y, t), t) = w -w q dz.
v 0 z (x, y, w(x, y, t), t) = ∂ t w + v LO . Then, equation (3.26) yields ∂ x w -w v 0 x dz + ∂ y w -w v 0 y dz + 2∂ t w -[v LO ] w -w = w -w q dz. (3.27)
The equation (3.27) is the Reynolds lubrication equation (see [START_REF] Szeri | Fluid film lubrication[END_REF]). It is usually written in terms of the pressure and we introduce K abs (which has no relationship with the absolute permeability) such that

w -w v 0 x dz = -K abs (x, y, t)∂ x p 0 and w -w v 0 y dz = -K abs (x, y, t)∂ y p 0 .
Using (3.24) and (3.25) yields

-K abs (x, y, t) = w -w (ξ + w)(w -ξ) η 0 (t, x, y, ξ) dξ - w -w ξ + w η 0 (t, x, y, ξ) dξ w -w w -ξ η 0 (t, x, y, ξ) dξ w -w 1 η 0 (t, x, y, ξ) dξ .
(3.28) In (3.28) K abs depends on η 0 non-locally and equation (3.27), and is written in the form of Reynolds' equation

2∂ t w -∂ x K abs ∂ x p 0 ) -∂ y K abs ∂ y p 0 ) = w -w q dz + [v LO ] w -w . (3.29)
3.1.7. Derivation of a PDE lubrication system for pressure and saturation. So far, the derivation in the previous section 3.1.6 has been performed without any additional assumptions on the two fluids. We now focus our attention on a specific case when η 0 can be downscaled. More precisely (see [START_REF] Mikelić | On the derivation of the Buckley-Leverett model from the two fluid Navier-Stokes equations in a thin domain computational geosciences[END_REF]), we assume

η ε = M χ {-a(x1,x2)≤z≤a(x1,x2)} + χ {-w≤z<-a(x1,x2) and a(x1,x2)<z≤w} , (3.30)
which means that the injecting fluid with viscosity η F,inj does not touch the fracture boundary w. Moreover, we note that it is compatible with the initial condition (3.12). The integrals involving η 0 can be evaluated as a function of the unknown free boundary a yielding

K abs = 2 3 a 3 M + w 3 -a 3 . (3.31) REMARK 3.2.
Permeability laws for single-phase fluids have been derived in [START_REF] Mikelić | A phase-field method for propagating fluidfilled fractures coupled to a surrounding porous medium[END_REF] for Newtonian and in [START_REF] Lee | Phase-field modeling of proppantfilled fractures in a poroelastic medium[END_REF] for non-Newtonian flows. For instance, M = 1 yields

K abs = 2w 3 3 . (3.32)
At this stage, it is convenient to introduce the saturation

s = a w , (3.33) and (3.31) becomes
K abs = 2w 3 3

s 3 M + 1 -s 3 . (3.34)
We denote that s F,res is the saturation for the residing fluid and s F,inj is the saturation for the injection fluid into the fracture. Also, p F,res is the pressure for residing fluid and p F,inj is the pressure for injection fluid in a fracture. For simplicity, in this section, we abuse the notation by defining

s := s F,inj , 1 -s := s F,res and
p 0 := (p 0 ) F,inj .
We recall that the 0 in the superscript denotes the order of our previous expansions. Here, we only work with the leading order 0 and neglect all further orders 1, 2, 3, . . .. Moreover, recall that

s F,res + s F,inj = 1 (3.35)
and

p c := p F,res -p F,inj , (3.36)
where p c is a capillary pressure which will be discussed in detail in Section 3.1.8.

We still miss one equation and we use the non-miscibility condition (3.2) and its zero order expansion approximation (3.22) to get

∂ t a + v 0 x (x, y, a, t)∂ x a + v 0 y (x, y, a, t)∂ y a -v z (x, y, a, t) = 0. (3.37)
Integrating (3.21), from a to w with respect to z, yields

-v z (x, y, a, t) = w a q dz -v 0 x (x, y, a, t)∂ x a -v 0 y (x, y, a, t)∂ y a -∂ x w a v 0 x dz -∂ y w a v 0 y dz -v LO -∂ t w. (3.38)
After plugging (3.38) into (3.37) and we obtain

∂ t a -∂ x w a v 0 x dz -∂ y w a v 0 y dz = v LO + ∂ t w. (3.39)
A direct calculation gives

w a v 0 x dz = - w 3 6 (1 -s) 2 (2 + s)∂ x p 0 = w a K dz ∂ x p 0 . ( 3
∂ t w(s -1) + ∇ x,y • w 3 6 (1 -s) 2 (2 + s)∇ x,y p 0 = v LO . (3.41) REMARK 3.3.
We note that the (3.41) is not valid in the fracture tip region and that the two fluid Stokes equations (3.3) -(3.6) do not simplify there. Computationally, we however do not solve the Stokes equations, but work with a cake region [START_REF] Kovalyshen | Fluid-driven fracture in poroelastic medium[END_REF] and appropriate interpolations that are motivated and explained in [START_REF] Mikelić | A phase-field method for propagating fluidfilled fractures coupled to a surrounding porous medium[END_REF](Remark 2.2) and [START_REF] Lee | Iterative coupling of flow, geomechanics and adaptive phase-field fracture including level-set crack width approaches[END_REF](Remark 3).

Our results can be summarized as follows: The lubrication approximation for the velocity and the pressure solutions to two phase Stokes system (3.3)-(3.6), (3.13)-(3.13), with the compatibility hypothesis (3.12), reads as follows:

FORMULATION 2. For s and 1 -s, we solve, respectively,

∂ t w(1 -s) -∇ x,y • w 3 3 (1 -s) 2 2 (2 + s)∇ x,y p 0 = -v LO + w a q in ω × (0, T ), (3.42)
∂ t (ws) -∇ x,y • w 3 3M (s 3 + 3s 2 (1 -s 2 )M )∇ x,y p 0 = a 0 q in ω × (0, T ) (3.43)
and for the leading order of the pressure p 0 , we solve

∂ t w -∇ x,y • w 3 3 (1 -s 3 + s 3 M ) ∇ x,y p 0 = w 0 q dz-v LO in ω × (0, T ). (3.44)
The initial conditions for the interface between the two fluids, the fracture width, and saturation are given by a| t=0 = a 0 , (3.45) w| t=0 = w 0 , (3.46) s| t=0 = a 0 /w 0 , in ω.

(3.47)

The relative permeabilities are given by

k F,res (s) := 1 2 (1 -s) 2 (2 + s), k F,inj (s) := s 3 + 3s(1 -s 2 ) 2 M, (3.48)
and the fractional flow curve is

F (s) := s 3 + 3s(1 -s 2 ) 2 M s 3 + (1 -s 3 )M . ( 3.49)
REMARK 3.4. We remark that the relative permeabilities are equivalent to the laws stated in [START_REF] Fourar | A viscous coupling model for relative permeabilities in fractures spe 49006[END_REF], which were compared to the experimental data from [START_REF] Fourar | Experimental study of air-water two-phase flow through a fracture (narrow channel)[END_REF]. The agreement was not very good and justified by an inaccuracy of the experimental data. On the other, the agreement was good when the saturation variable was replaced by the Lockhart-Martinelli factor χ LM = k F ;res k F ;inj . Dropping the surface pressure between fluids could also be a source of inaccuracy and our solution consists in adding the capillary pressure introduced in [START_REF] Babchin | On the capillary pressure function in porous media based on relative permeabilities of two immiscible fluids[END_REF] to the model. REMARK 3.5. Under assumptions that the BV-norm of η ε is uniformly bounded with respect to ε, a justification of the upscaled model was obtained by Paoli in [START_REF] Paoli | Asymptotic behavior of a two fluid flow in a thin domain: From stokes equations to Buckley-Leverett equation and Reynolds law[END_REF].

3.1.8. The capillary pressure function for a penny shape fracture. In the previous section, we have assumed that the capillary pressure is zero. Knowledge of the relative permeabilities is crucial for the numerical simulation of the multiphase flows through porous media. Nevertheless, the models where the dissipation is absence lead to the Buckley-Leverett equation, which is a hyperbolic conservation law. The solutions lose regularity in finite time and discontinuous solutions are not unique. In order to have uniqueness, the Oleinik entropy criterion (or equivalently the Welge method [START_REF] Welge | A simplified method for computing oil recovery by gas or water drive[END_REF] in the petroleum engineering) is employed; see [START_REF] Concus | Numerical solution of a nonlinear hyperbolic equation by the random choice method[END_REF]. Nevertheless, recent questions about entropies associated with two phase flow and the application of Oleinik entropy is arisen. For more details we refer to the analytical work by van Duijn, Peletier and Pop [START_REF] Van Duijn | A new class of entropy solutions of the Buckley-Leverett equation[END_REF] justifying experimental observations by DiCarlo in [START_REF] Dicarlo | Experimental measurements of saturation overshoot on infiltration[END_REF].

A natural remedy is to add the physical dissipation linked to the surface tension. Classically it is taken into account through the capillary pressure. A more recent approach is to use a Cahn-Hilliard type phase-field equations, as in [START_REF] Cueto-Felgueroso | A phase-field model of two-phase Hele-Shaw flow[END_REF], which involves a fourth order differential operator in the flow equations. This is non-standard for the complex flows simulations with industrial applications. We choose to calculate the capillary pressure function from our relative permeabilities, using the interface energy approach from [START_REF] Babchin | On the capillary pressure function in porous media based on relative permeabilities of two immiscible fluids[END_REF].

Following the approach from [START_REF] Babchin | On the capillary pressure function in porous media based on relative permeabilities of two immiscible fluids[END_REF], we use as capillary pressure:

p c (s) = γ res cos θ φ K c (k F,res (s)) -1/2 + (k F,inj (s)) -1/2 (1 -s) + C, (3.50)
where γ res is the phase surface energy between two different phases, θ is the injecting-phase contact angle, and

K c = 2H 2 w 2

3

. The constant C corresponds to initial displacement pressure For imbibition, we set C = 0 and k F,res and k F,inj are given by (3.48).

3.1.9. The final 2D lubrication system: an upscaled two phase model for fracture flow. After including the capillary pressure, our upscaled dimensional model for the two phase flow in the fracture Ω F reads FORMULATION 3. Find the saturations 1 -s and s such that:

∂ t h(1 -s) -∇ x,y • h 3 3η F,res k F,res ∇ x,y p F,res = -v out + h sh q dx 3 in ω × (0, T ), (3.51)
∂ t (hs) -∇ x,y • h 3 3η F,inj k F,inj ∇ x,y p F,inj = sh 0 q dx 3 in ω × (0, T ), (3.52)
and find the pressure p F,inj such that

∂ t h -∇ x,y • h 3 3 k F,res η F,res + k F,inj η F,inj ∇ x,y p F,inj + h 3 3 k F,res η F,res ∇ x,y p c (s) = 1 2 h -h q dx 3 -v out in ω × (0, T ), (3.53)
or find the pressure p F,res such that

∂ t h -∇ x,y • h 3 3 k F,res η F,res + k F,inj η F,inj ∇ x,y p F,res - h 3 3 k F,inj η F,inj ∇ x,y p c (s) = 1 2 h -h q dx 3 + 1 2 [v LO ] w -w in ω × (0, T ) (3.54)
where

k F,res η F,res + k F,inj η F,inj = s 3 η F,inj + 1 -s 3 η F,res ,
with the initial conditions

s| t=0 = s 0 = A 0 /h 0 , (3.55
) h| t=0 = h 0 .

(3.56)

Here the phase velocities are

1 2h -sh -h + h sh (v 0 x , v 0 y ) = 1 h h sh (v res x , v res y )dx 3 = - K F η F,res k F,res (s)∇ x,y p F,res , (3.57) 1 2h -sh -h + h sh (v 0 x , v 0 y ) = 1 h sh 0 (v inj x , v inj y )dx 3 = - K F η F,inj k F,inj (s)∇ x,y p F,inj , (3.58)
with the absolute permeability 3.1.10. A 3D two phase flow system. We now use (3.51) -(3.53) and substitute them into a classical three dimensional two phase flow model in a fracture. The background of this procedure is similar to the ideas presented in [START_REF] Lee | Phase-field modeling of proppantfilled fractures in a poroelastic medium[END_REF] and [START_REF] Mikelić | A phase-field method for propagating fluidfilled fractures coupled to a surrounding porous medium[END_REF] for single-phase flow and quasi-Newtonian flow, respectively. In our global systems, we recall that we have two variables for the saturation in a fracture denoted by s F,inj (the saturation for the residing fluid) and s F,res (saturation for the injection fluid), and three variables for the pressure denoted by p R , p F,inj , p F,res , where p R is the pressure for the reservoir, p F,res is the pressure for residing fluid in a fracture and p F,inj is the pressure for the injection fluid in a fracture. However after employing the relation s F,inj + s F,res = 1, and p c = p F,res -p F,inj ,

K F = h 2 3 .
where p c is referred as the capillary pressure in a fracture, we only solve for the saturation and the pressure of the injecting fluid, s F,inj , p F,inj in a fracture as described in Remark 3.6.

In the reservoir, we have p R , but a saturation variable is not defined since we assume a single phase fluid in the reservoir; see Figure 3.3 for more details. Thus, we obtain the following Formulation 4 from Formulation 3 by rewriting (3.52) and (3.53).

FORMULATION 4. Find s F,inj such that

∂ t (φ F s F,inj ) -∇ • K F k F,inj (s F,inj ) η F,inj ∇p F,inj = q F,inj in Ω F × (0, T ), (3.60a) K F k F,inj (s F,inj ) η F,inj ∇p F,inj • n F = 0 on Γ F (t) × (0, T ), (3.60b)
where φ F is the porosity of the fracture which is set to one, and η F,res ,η F,inj , are the viscosities for each fluid, respectively. Moreover, q F,inj := sF,inj qF,inj , where sF,inj , qF,inj are the saturation injection/production terms and flow injection/production, respectively. If qF,inj > 0, sq,inj is the injected saturation of the fluid and if qF,inj < 0, sq,inj is the produced saturation. We note that sq,inj + sq,res = 1.

In addition, find the pressure p F,inj such that

-∇ • K F λ tot ∇p F,inj + K F λ F,res ∇p c (s F,inj ) = (q F ) tot in Ω F (t) × (0, T ), (3.61)
Here the total mobility, the mobility for the injecting fluid, and the mobility for the residing fluid, are defined respectively:

λ tot := λ tot (s F,inj ) = λ F,inj (s F,inj ) + λ F,res (s F,inj ), (3.62a)
λ F,inj := λ F,inj (s F,inj ) = k F,inj (s F,inj ) η F,inj , (3.62b)
λ F,res := λ F,res (s

F,inj ) = k F,res (s F,inj ) η F,res , (3.62c)
and (q F ) tot := q F,inj + q F,res .

(3.62d)

The above formulation is similar to the two phase flow model in porous media as discussed in [START_REF] Arbogast | A discontinuous Galerkin method for two-phase flow in a porous medium enforcing H (div) velocity and continuous capillary pressure[END_REF][START_REF] Chen | Computational methods for multiphase flows in porous media[END_REF][START_REF] El-Amin | An iterative implicit scheme for nanoparticles transport with two-phase flow in porous media[END_REF][START_REF] Hoteit | Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures[END_REF][START_REF] Lee | Enriched Galerkin approximations for two phase flow in porous media with capillary pressure[END_REF]], but we have computed new relative permeabilities (3.48) and employ a capillary pressure (3.50) to consider two phase flow in a fracture.

REMARK 3.7. Similar to the literature [START_REF] Witherspoon | Validity of cubic law for fluid flow in a deformable rock fracture[END_REF][START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF], the original lubrication law has a cubic power in which the width parameter h in the time derivative replaces the porosity. By dividing the flux over the cross-sectional area (see, for example, page 5 in [START_REF] Zimmerman | Hydraulic conductivity of rock fractures[END_REF]), we 'lose' one order and arrive at a quadratic law.

4. Single-phase flow in the reservoir. We briefly recapitulate the standard single phase flow model that we employ in the porous medium. Continuity of mass is modeled as:

∂ t (ρ R φ R ) + ∇ • (ρ R v R ) = ρ R q R in Ω R × (0, T ). (4.1)
Here ρ R is the fluid density in reservoir and q R is a source/sink term. We assume the fluid in the reservoir is slightly compressible, thus we define the fluid density as

ρ R := ρ 0 R exp(c R (p R -p 0 R )) ≈ ρ 0 R [1 + c R (p R -p 0 R )], (4.2)
where p R : Ω R × [0, T ] → R is the pressure, p 0 is the initial pressure at t = 0, ρ 0 R is the reference density and c R is the fluid compressibility. In addition, φ R is the reservoir fluid fraction and

φ R := φ 0 + α∇ • u + 1 M B (p R -p 0 ). (4.3)
Here u : Ω×[0, T ] → R d is the solid displacement, α ∈ [0, 1] is the Biot coefficient, M B > 0 is a given Biot modulus, and φ 0 is the initial fluid fraction value. Next, we describe the flow given by Darcy's law:

v R = - K R η R ∇p R , (4.4)
where η R is the fluid viscosity, K R is the permeability for the reservoir.

Following the general reservoir approximation with the assumption that c R is small enough, we use ρ R = ρ 0 R and assume p 0 = 0, to rewrite the (4.1) by

∂ t ( 1 M B p R + α∇ • u) -∇ • K R η R ∇p R = q R in Ω R × (0, T ). (4.5)
The system is supplemented with initial and boundary conditions. The initial conditions are given by:

p R (x, 0) = p 0 R for all x ∈ Ω R (t = 0),
where p 0 R is a smooth given pressure. To determine the initial location of Γ F at t = 0, we also have φ(x, 0) = φ 0 for all x ∈ Λ(t = 0), where φ 0 is a given smooth initial fracture. We finally prescribe the boundary and interface conditions for pressure as

K R ∇p R • n = 0 on ∂Λ × (0, T ), (4.6a
) p R = p F,res on Γ F × (0, T ), (4.6b) K R η R ∇p R • n F = K F k F,res η F,res ∇p F,res • n F on Γ F × (0, T ), (4.6c)
where ρ F,res is the residing fluid density in fracture.

5.

Coupling algorithms, variational formulations, and discretization. This section presents frameworks for the numerical solution of the previously introduced equations. We introduce the following equation settings: a displacement/phase-field framework, a crackwidth problem, a pressure diffraction problem and a saturation formulation. For coupling flow (pressure and saturation) and geomechanics (displacement/phase-field), we employ iterative fixed-stress coupling. Moreover, for consistent modeling of the flow equations, we introduce diffraction systems for which we explain the decomposition of the domain Λ(t) into appropriate time-dependent sub-domains.

Preliminaries.

For the spatial discretization, we utilize Galerkin finite element methods. In particular, we employ continuous Galerkin finite elements (CG) for displacements and phase-field, and an enriched Galerkin finite elements (EG) for both the pressure and saturation equations in order to obtain local mass conservation with an efficient solver. In porous media flow it is well known that the velocities should be locally mass conservative to avoid unphysical spurious oscillations. This can either be achieved by post-processing of the velocities [START_REF] Sun | Projections of velocity data for the compatibility with transport[END_REF] or by a direct solution from the appropriate function space. We employ the latter method, the enriched Galerkin (EG), that was first proposed in [START_REF] Sun | A locally conservative finite element method based on piecewise constant enrichment of the continuous galerkin method[END_REF] for elliptic problems and was extended to parabolic and hyperbolic cases in [START_REF] Lee | Enriched Galerkin approximations for two phase flow in porous media with capillary pressure[END_REF][START_REF] Lee | Adaptive enriched Galerkin methods for miscible displacement problems with entropy residual stabilization[END_REF][START_REF] Lee | A locally conservative enriched Galerkin approximation and efficient solver for elliptic and parabolic problems[END_REF] with efficient solvers and dynamic mesh adaptivity. We consider a shape regular mesh family and we assume that each mesh T is a subdivision of Λ made of disjoint elements K, i.e., squares when d = 2 or cubes when d = 3. Each subdivision is assumed to exactly approximate the computational domain, thus Λ = ∪ K∈T K. The diameter of an element K ∈ T is denoted by h K and we denote h min and h max for the minimum and maximum diameters, respectively.

The enriched Galerkin method requires (similar to discontinuous Galerkin) additional notations due to face term contributions. We let E h be the set of all edges/faces and E o h and E ∂ h are the collection of all interior and boundary edges/faces, respectively. In the following, we use the notation for edges in two dimension but the results hold analogously for faces in three dimension. The boundary edges E ∂ h can be further decomposed into

E ∂ h = E D,∂ h ∪ E N,∂ h , where E D,∂ h
is the collection of edges where Dirichlet boundary conditions are imposed, while E N,∂ h is the collection of edges where Neumann boundary conditions are imposed. In addition, we denote that

E 1 h := E o h ∪ E D,∂ h .
For any e ∈ E o h , let K + and K -be two neighboring elements such that e = ∂K + ∩ ∂K -. We denote by h e the maximum length of the edges e. Let n + and n -be the outward normal unit vectors to ∂K + and ∂K -, respectively (n ± := n |K± ). For any given function, vector, or matrix ξ defined on the triangulation T , we denote ξ ± by the restrictions of ξ to K ± . The space H k (T ) (k ∈ IR) is the set of element-wise H k functions on T , and L 2 (E h ) refers to the set of functions whose traces on the elements of E h are square integrable. In addition, we define the average { {•} } as follows:

{ {w} } := 1 2 (w + + w -) for w ∈ [L 2 (T )] d , (5.1)
{ {w} } := 1 2 (w + + w -) for w ∈ [L 2 (T )] d×d , (5.2)
and { {w} } := w for e ∈ E ∂ h .

(

The jump across the interior edge will defined as For any integer k ≥ 1 and any K ∈ T , we denote by Q k (K) the space of scalar-valued multivariate polynomials over K of partial degree of at most k. The vector-valued counterpart of

Q k (K) is denoted Q Q Q k (K).
Temporal discretization is performed with an A-stable backward Euler finite difference scheme.

Here, we define a partition of the time interval 0 =: t 0 < t 1 < • • • < t N := T and denote the time step size by δt := t n -t n-1 .

Global algorithms and fixed stress iterative coupling.

To facilitate the understanding of the overall numerical solution, we first present in this section the global coupling algorithm based on a fixed-stress iteration. The fixed-stress coupling scheme is well-known in subsurface modeling, environmental and petroleum engineering problems [START_REF] Castelletto | Accuracy and convergence properties of the fixed-stress iterative solution of two-way coupled poromechanics[END_REF][START_REF] Kim | Stability, accuracy, and efficiency of sequentiel methods for flow and geomechanics[END_REF][START_REF] Kim | Stability and convergence of sequentiel methods for coupled flow and geomechanics: fixed-stress and fixed-strain splits[END_REF][START_REF] Mikelić | Convergence of iterative coupling for coupled flow and geomechanics[END_REF][START_REF] Mikelić | Numerical convergence study of iterative coupling for coupled flow and geomechanics[END_REF][START_REF] Settari | A coupled reservoir and geomechanical simulation system[END_REF][START_REF] Settari | Advances in coupled geomechanical and reservoir modeling with applications to reservoir compaction[END_REF]. The idea is to iterate between flow and geomechanics at each time step until sufficient accuracy is achieved. For detailed explanations and numerical demonstrations using phase-field fracture, we refer to [START_REF] Lee | Iterative coupling of flow, geomechanics and adaptive phase-field fracture including level-set crack width approaches[END_REF]. The algorithm used in this paper is provided in Figure 5.1. A theoretical justification of the fixed-stress phase-field fracture scheme in which the phase-field values are assumed to be given has been recently shown in [START_REF] Almani | Multirate coupling for flow and geomechanics applied to hydraulic fracturing using an adaptive phase-field technique[END_REF]. Computational evidence for several numerical tests, including a fluid filled Sneddon test [START_REF] Sneddon | Crack problems in the classical theory of elasticity[END_REF] have been performed in [START_REF] Lee | Iterative coupling of flow, geomechanics and adaptive phase-field fracture including level-set crack width approaches[END_REF]. 5.2.1. Notations for spatial, temporal and fixed-stress discretizations. For the following, we introduce the necessary notation. We denote by l the fixed-stress iteration index and by n the time step index. Thus, the spatial approximations of the functions u(x, t), ϕ(x, t), ϕ LS (x, t), w(x, t), p(x, t), s(x, t) are denoted by U, Φ, Φ LS , W, P, S.

The functions ϕ LS (x, t), w(x, t) and Φ LS , W , respectively, have not been introduced yet.

They stand for a level-set function and a width approximation problem introduced in [START_REF] Lee | Iterative coupling of flow, geomechanics and adaptive phase-field fracture including level-set crack width approaches[END_REF]. Details of the numerical solution are provided below in Section 5.4.

The approximation of each function in time

0 ≤ n ≤ N are denoted by U n , Φ n , Φ n LS , W n , P n , S n ,
where U n := U(t n ) at time t n and so forth. At each time step n, we have an inner loop indexed by l = 1, 2, 3, • • • and we denote the fixed stress iterative solutions for each time step by U l := U n,l , Φ l := Φ n,l , Φ l LS := Φ n,l LS , W l := W n,l , P l := P n,l , S l := S n,l .

Moreover, we set as initial fixed-stress step U n,0 := U n-1 and so forth.

The final algorithm.

The overall solution loop is presented in Algorithm 1.

Displacements and phase-field.

Based on Formulation 1, we derive a variational incremental formulation. Therein, the continuous irreversibility constraint ∂ t ϕ ≤ 0 is approximated by ϕ ≤ ϕ old .

The discrete solution variables are denoted by U and Φ (approximating u and φ), respectively. In more detail, we have

U ∈ C 0 ([0, T ]; V V V(T )) and Φ ∈ C 0 ([0, T ]; Z(T )), where V V V(T ) := {V ∈ C 0 ( Λ; R d ) | V = 0 on ∂Λ, V | K ∈ Q Q Q 1 (K), ∀K ∈ T }, (5.4
)

Z(T ) := {Z ∈ C 0 ( Λ; R)| Z n ≤ Z n-1 ≤ 1, Z| K ∈ Q 1 (K), ∀K ∈ T }. (5.5)
Moreover, we extrapolate Φ (denoted by E(Φ)) in some terms in order to avoid an indefinite Hessian matrix [START_REF] Heister | A primal-dual active set method and predictorcorrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF]. That is to say, we take the two previous time step solutions at t n-1 and t n-2 and evaluate at time t:

E(Φ n,l ) = Φ n-2 + (t -t n-1 -t n-2 ) (t -t n-1 ) -(t -t n-1 -t n-2 ) (Φ n-1 -Φ n-2 ).
Algorithm 1 Iterative coupling for two phase flow phase-field fractures including level-set and width computations Initial conditions : Let Φ n-1 , U n-1 , P n-1 , S n-1 , Φ n-1 LS , W n-1 be given. At time t n : repeat For l = 1, 2, 3 . . .:

-Solve the pressure diffraction Formulation 7 for the reservoir pressure and for the injecting phase P l = P l F,inj -Solve the saturation diffraction Formulation 8 for the injecting phase S l = S l F,inj -Solve the quasi-monolithic displacement/phase-field Formulation 5 for (U l , Φ l ) -Compute the level-set Φ l LS and solve the width Formulation 6 for W l until the stopping criterion for fixed-stress split is satisfied:

max{ P l -P l-1 , U l -U l-1 , Φ l -Φ l-1 } ≤ TOL FS , TOL FS > 0 Set: (P n , U n , Φ n , S n ) := (P l , U l , Φ l , S l ).
The other two variables Φ n LS and W n are obtained from Φ n and U n . Increment the time n → n + 1.

For the solution process, we define a semi-linear form in the following, yielding: FORMULATION 5. Let us assume that P l := P n,l is a given pressure at the time t n and at the fixed-stress iteration l. Given the initial conditions U l-1 and Φ l-1 we seek

{U l , Φ l } ∈ V V V(T ) × Z(T ) such that A(U l , Φ l )(w, ψ -Φ l ) ≥ 0, ∀{w, ψ} ∈ V V V(T ) × Z(T ), (5.6)
where

A(U l , Φ l )(w, ψ -Φ l ) := Λ (1 -k)(E(Φ l ) 2 + k)σ + (U l ) : e(w) dx + Λ σ -(U l ) : e(w) dx - Λ (α -1)E(Φ l ) 2 P l ∇ • w dx + Λ E(Φ l ) 2 ∇P l • w dx + (1 -k) Λ Φ l σ + (U l ) : e(U l ) • (ψ -Φ l ) dx -2(α -1) Λ Φ l P l ∇ • U l • (ψ -Φ l ) dx + Λ 2Φ l ∇P l • U l • (ψ -Φ l ) dx -G c Λ 1 ε (1 -Φ l ) • (ψ -Φ l ) dx + G c Λ ε∇Φ l • ∇(ψ -Φ l ) dx.
The solution of this nonlinear variational inequality is briefly explained in Section 6 with all details presented in [START_REF] Heister | A primal-dual active set method and predictorcorrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF].

Level set and interpolation width.

For the solution of the fracture pressure equation, we need the permeabilities, which depend on the fracture opening displacements. Based on the solution of Formulation 5, we are able to compute the fracture opening displacements (width). We briefly recapitulate the main ideas in this section.

The spatially discretized solution variables for the level-set and the width are denoted by Φ LS (x, t) and W (x, t), respectively. Those functions are approximated by using continuous piecewise polynomials given in their respective finite element spaces,

V LS (T ) := {V LS ∈ C 0 ( Λ; R) | V LS = 0 on Γ F , V LS | K ∈ Q 1 (K), ∀K ∈ T }
for the level-set and

V w (T ) := {V w ∈ C 0 ( Λ; R) | V w = 0 on ∂Λ, V w | K ∈ Q 1 (K), ∀K ∈ T }
for the width. After the computation of the displacement field U n,l and the phase-field Φ n,l at time t n and iteration number l, the level set function introduced in (2.10) is simply computed as follows:

Φ n,l LS = Φ n,l -C LS . (5.7)
Details of the entire procedure and different algorithms to consider the level set evaluation are provided in [START_REF] Lee | Iterative coupling of flow, geomechanics and adaptive phase-field fracture including level-set crack width approaches[END_REF][START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF]. As described in Section 2.3, we compute the width W l D on the fracture boundary Γ F by

2W l D = -[U l • n F ] = -2U l • ∇Φ l LS ∇Φ l LS .
(5.8) REMARK 5.1. We recall for the convenience of the reader that our definition of the width W l D concerns only one half of the fracture (see Figure 3.1). Thus, in order to get the full width, we need to multiply by the factor 2.

Finally, we interpolate the above computed width values W l D on Γ F to values inside the fracture (see Remark 2.2), and denote the solution by W l . Here, we briefly introduce the weak formulation for the interpolation process to solve for W l as described in [START_REF] Lee | Iterative coupling of flow, geomechanics and adaptive phase-field fracture including level-set crack width approaches[END_REF]:

FORMULATION 6. Find W l ∈ C 0 ([0, T ]; V w (T )) such that A W (W l , ψ) = F W (ψ) ∀ψ ∈ V w (T )
where

A W (W l , ψ) = (∇W l , ∇ψ) + θ Γ n,l F W l ψ ds, F W (ψ) = θ Γ n,l F W l D • ψ ds.
5.5. Decomposing the domain Λ into Ω R , Ω F , and Ω ε . Based on the previous two subsections, we have obtained U l , Φ l and W l . Then1 we increment l → l + 1; see Figure 5.1. We now address the coupling of flow and saturations between the fracture and the porous medium yielding P l+1 and S l+1 .

A conceptional challenge is that we consider both diffusive and sharp interfaces. We assume a sharp interface for the pressure coupling and modeling but the phase-field fracture framework yields a diffusive (smeared) interface (see Figure 5.2). To this end, we partition the transition zone (Ω ε ) into a fracture part (Ω ε,F ) a porous medium part (Ω ε,R ), which we describe in this section. This decomposition is the basis for so-called diffraction systems in which the same type of equations is imposed on the entire domain, but the different parts (fracture and porous media) are distinguished by time-dependent and spatially varying coefficients.

First, we define the fracture domain Ω F and the reservoir domain Ω R by introducing following two linear indicator functions χ F and χ R for the two different sub-domains as

χ R := χ R (Φ) = 1 in Ω R (t), and χ R = 0 in Ω F (t),
(5.9)

χ F := χ F (Φ) = 1 in Ω F (t), and χ F = 0 in Ω R (t).
(5.10)

The transition zone Ω ε is characterized as follows:

χ R + χ F = 1 in Ω ε ,
where χ R = 1 and χ F = 1 and have intermediate values between 0 and 1.

By employing the already computed phase-field values Φ n,l as an indicator function, we can specifically define

χ n,l R := 1 + H h (Φ n,l ) 2 , and χ n,l F := 1 -H h (Φ n,l ) 2 , (5.11)
where H h (.) is an approximation of the Heaviside function given as

H h (w) =     1 if w > 0.5 + D , -1 if w < 0.5 -D , w -0.5 D otherwise.
(5.12)

The transition zone Ω ε is characterized by the third case. Here, we set D ε = 0.1 as discussed in Section 5.4.

5.6. Discretization of flow/pressure using enriched Galerkin (EG). The EG space approximation P (x, t) of the pressure function is approximated by piecewise polynomials in the following finite element space:

V EG k := M k 0 (T ) + M 0 (T ), (5.13)
where

M k (T ) := M ∈ L 2 (Λ)| M | K ∈ Q k (K), ∀K ∈ T , (5.14)
and

M k 0 (T ) := M k (T ) ∩ C(Λ). (5.15)
The space M k 0 (T ) is the subspace of M k (T ) consisting of the globally continuous piecewise continuous functions. Here we choose linear finite elements with k = 1; see [START_REF] Becker | A reduced P1-discontinuous Galerkin method[END_REF][START_REF] Lee | Enriched Galerkin approximations for two phase flow in porous media with capillary pressure[END_REF][START_REF] Lee | Adaptive enriched Galerkin methods for miscible displacement problems with entropy residual stabilization[END_REF][START_REF] Lee | A locally conservative enriched Galerkin approximation and efficient solver for elliptic and parabolic problems[END_REF][START_REF] Sun | A locally conservative finite element method based on piecewise constant enrichment of the continuous galerkin method[END_REF] for more details. Using Section 5.5 and employing phase-field as an indicator function, we can couple the two phase flow equation (3.61) in the fracture with the single phase equation (4.5) in the reservoir in terms of a pressure diffraction system. The following discretization is obtained using a standard derivation by multiplying the pressure diffraction system [START_REF] Lee | Phase-field modeling of proppantfilled fractures in a poroelastic medium[END_REF][START_REF] Lee | Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model[END_REF][START_REF] Mikelić | A phase-field method for propagating fluidfilled fractures coupled to a surrounding porous medium[END_REF] by a discontinuous test function in V EG k , integration by parts and summing over elements . We recall our notation introduced in Section 5.2.1, namely that P n := P n F,inj and P l := P n,l . We focus on the specifics of the fixed-stress iteration between flow and geomechanics/fracture as describe in Section 5.2 Then, the discretized problem reads as follows: FORMULATION 7. At time t n and the fixed-stress level l + 1, let U l , Φ l , W l , S l be given. Furthermore, let the previous time step solutions U n-1 , P n-1 be given. We iterate for l = 1, 2, 3, . . . such that: Find P l+1 ∈ V EG 1 :

A P (P l+1 , ψ) = F P (ψ) ∀ ψ ∈ V EG 1 , (5.16)
where the Incomplete Interior Penalty Galerkin (IIPG) is employed to be compatible with the saturation system [START_REF] Dawson | Compatible algorithms for coupled flow and transport[END_REF][START_REF] Scovazzi | Analytical and variational numerical methods for unstable miscible displacement flows in porous media[END_REF]. Here the variational form is defined as

A P (P l+1 , ψ) := K∈T K θ l P l+1 -P n-1 δt ψ dx + K∈T K K l eff ∇P l+1 • ∇ψ dx - e∈E 1 h e { {K l eff ∇P l+1 } } [[ψ]] dγ + e∈E 1 h e α p h e [[P l+1 ]] [[ψ]] dγ, (5.17)
where α p > 0 is a penalty parameter and the right hand side is defined as

F P (ψ) := K∈T K q l ψ dx + K∈T K f l ψ dx - K∈T K θ l c ∇P c (S l ) • ∇ψdx + e∈E 1 h e { {θ l c ∇P c (S l ))} } [[ψ]] dγ. (5.18)
Here the diffraction coefficients defined as

θ l := χ R (Φ l ) 1 M B + 3α 2 3λ + 2µ + χ F (Φ l )(c F ), (5.19)
K l eff := χ R (Φ l ) K R η R + χ F (Φ l ) K F (W l )λ tot (S l ) , (5.20)
q l := χ R (Φ l ) (q n R ) + χ F (Φ l )(q F ) n tot , (5.21)
f l := χ R (Φ l ) 3α 2 3λ + 2µ P l -P n-1 δt -α∇ • U l -U n-1 δt (5.22) + χ F (Φ l ) • 0, θ l c := χ R (Φ l ) • 0 + χ F (Φ l ) K F (W l )λ F,res (S l ) . (5.23)
Here we assume no flow boundary conditions on ∂Λ.

REMARK 5.2. For the stability of the numerical solution, we have added (as it is very often done for such problems) the time derivative ∂ t χ F (Φ l )c F p in Formulation 7, which corresponds to assuming a slightly compressible fluid with the fracture fluid compressibility c F .

Finally, for the saturation approximation in the next section, we introduce flux variables.

Here, the conservative flux variables of the injected fluid in fractures, V l+1 (U l , W l , S l , P l+1 ), can be obtained as shown in [START_REF] Lee | Enriched Galerkin approximations for two phase flow in porous media with capillary pressure[END_REF][START_REF] Lee | A locally conservative enriched Galerkin approximation and efficient solver for elliptic and parabolic problems[END_REF][START_REF] Sun | A locally conservative finite element method based on piecewise constant enrichment of the continuous galerkin method[END_REF]:

V l+1 | T := -K l eff ∇P l+1 , ∀K ∈ T (5.24) V l+1 • n| e := -{ {K l eff ∇P l+1 } } • n + h -1 e [[P l+1 ]], ∀e ∈ E o h , (5.25)
where n is the unit normal vector of the boundary edge e of K.

Numerical discretization of saturation.

In this subsection, we provide the enriched Galerkin approximation of the saturation Formulation 4 described in Section 3.1.10. The space approximation of the injected saturation s F,inj (•, t) is denoted by S F,inj (x, t) using the same linear EG finite element space V EG 1 as the pressure system. As short hand notation, we denote the approximation of S F,inj (x, t n ) by S n := S n F,inj and as before S l := S n,l . As for the pressure system, we formulate the saturation system in the spirit of a diffraction problem, despite the fact that the saturation is only present in the fracture in this work. However, this formalism allows us an easy extension to two phase flow in the reservoir once such a model shall be employed. FORMULATION 8. Let the previous fixed-stress solutions U l , Φ l , W l , P l+1 and the previous time step solutions S n-1 be given. At time t n we iterate for l = 1, 2, 3, . . .:

Find S l+1 ∈ V EG 1 such that A S (S l+1 , ψ) = F S (ψ), ∀ ψ ∈ V EG 1 , (5.26)
where W is the IIPG IP method and the variational form defined as

A S (S l+1 , ψ) := K∈T K θ l s S l+1 -S n-1 δt ψdx + K∈T K µ l+1 stab ∇S l+1 • ∇ψdx - e∈E 1 h e { {µ l+1 stab ∇S l+1 } } [[ψ]] dγ + e∈E 1 h e α s h e { {µ stab } } |e [[S l+1 ]] [[ψ]] dγ, (5.27)
where α s is a penalty parameter. We note that terms weighted with µ l+1 stab are artificial diffusion terms for numerical stabilization to ensure s ∈ [0, 1] and which avoid spurious oscillations for the advection system. Then,

F S (ψ) := K∈T K V l+1 ∇ψdx - e∈E 1 h e V l+1 • n [[ψ]] dγ + K∈T K q l s ψ dx.
(5.28)

The artificial viscosity is defined as

µ l+1 stab := c s h K V l+1 L ∞ (K) ,
where c s is a positive constant. For more details, the readers are referred to [START_REF] Lee | Enriched Galerkin approximations for two phase flow in porous media with capillary pressure[END_REF]. Here the diffraction coefficients are defined as

θ l s := χ R (Φ l ) • 0 + χ F (Φ l ) (φ F ) , (5.29)
q l s := χ R (Φ l ) • 0 + χ F (Φ l ) q n F,inj . (5.30)
6. Parallel solvers and mesh adaptivity. We briefly discuss the numerical solvers to treat the discretized Formulations 5, 6, 7, 8. The overall problem is a large coupled framework in three dimensions which is treated by parallel computing and local mesh adaptivity in order to keep the computational costs reasonable.

We describe the solvers in the following with more details: both, the pressure and saturation diffraction problems, are solved with generalized minimal residual method (GMRES) solvers with diagonal block-preconditioning. The nonlinear quasi-monolithic displacement/phasefield system presented in Formulation 5 is solved with Newton's method and line search algorithms. The constraint minimization problem is treated with a semi-smooth Newton method (i.e., a primal-dual active set method). Both methods are combined in one single loop leading to a robust and efficient iteration scheme that is outlined in [START_REF] Heister | A primal-dual active set method and predictorcorrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF]. Within Newton's loop we solve the linear equation systems with GMRES solvers with diagonal block-preconditioning from Trilinos [START_REF] Heroux | An Overview of Trilinos[END_REF]. Finally, the linear-elliptic level-set and width problems are solved with a parallel conjugate gradient (CG) solver and symmetric successive over-relaxation (SSOR) preconditioning.

To allow for small phase-field regularization parameters ε, we employ a predictor-corrector scheme that chooses an initial ε > h K at the beginning of the computation [START_REF] Heister | A primal-dual active set method and predictorcorrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF]. Since this is a model parameter, we do not want to change the model during the computation and therefore keep ε fixed. However, the crack propagates and in coarse mesh regions ε > h K may be violated. Then, we take the first step as a predictor step, and refine the mesh such that ε > h K holds again and finally recompute the solution. It has been shown in [START_REF] Heister | A primal-dual active set method and predictorcorrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF] (two dimensional) and [START_REF] Lee | Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model[END_REF][START_REF] Wick | 3D phase-field for pressurized fracture propagation in heterogeneous media[END_REF] (three dimensional) that this procedure is efficient and robust. As refinement criterion we take simply a threshold value of the phase-field variable, i.e., once Φ < C A with 0 < C A < 1, e.g., C A ∼ 0.8, we flag a mesh element for refinement.

The implementation is done in our software IPACS (integrated parallel advanced crack simulator), which is an extension of [START_REF] Lee | Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model[END_REF][START_REF] Wick | 3D phase-field for pressurized fracture propagation in heterogeneous media[END_REF] based on the open-source MPI parallel phase-field fracture framework [START_REF] Heister | A primal-dual active set method and predictorcorrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF] using deal.II [START_REF] Bangerth | The deal.II library, version 8.4[END_REF], p4est [START_REF] Burstedde | P4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees[END_REF] and Trilinos [START_REF] Heroux | An Overview of Trilinos[END_REF].

Numerical examples.

In this section, we present several numerical examples to substantiate our proposed algorithms and models.

Example 1:

A fracture propagation in a homogeneous media. In this first example, we consider a single propagating fracture. A first fluid, the residing fluid, is already present. We then inject a second fluid. This example is split into two subsections. In the first test (Example 1a), the crack is aligned with the x-axis and thus finally on existing mesh faces. In the second example (Example 1b) the initial crack is shifted by 45 degree in order to show that phase-field fracture allows for arbitrary initial cracks, not necessarily aligned with the mesh, but leading to the same findings (same crack length, maximal pressure, etc.). This is a very important feature of phase-field fracture and shows that our entire derivation holds irrespectively of the initial geometric configuration. ample 1b. The initial mesh is 5 times uniformly refined. We then allow for three levels of predictor-corrector mesh refinement.

Boundary and initial conditions. For the displacements we prescribe homogeneous Dirichlet conditions on ∂Λ. For the phase-field and the pressure system, homogeneous Neumann conditions are employed. The boundary conditions for the crack-width problem have been explained in Section 5.4. The initial phase-field values are set to zero for the initial fracture described above and ϕ = 1 otherwise. Also the initial displacement and pressure values are set to zero.

Parameters. The time step size is chosen as δt = 0.01 s. The total number of time steps is 1000 with the final time 10 s. The phase-field value for the predictor-corrector refinement is chosen as 0.8; thus below 0.8 the mesh will be refined up to three times. The penalization parameter for the enriched Galerkin pressure system is set to 100. The penalization parameter for the saturation system is set to 1000. The capillary coefficient is γ res = 1. Furthermore, T OL F S = 10 -3 and the Newton tolerance for the phase-field displacement system is 10 -10 .

The regularization parameters are chosen as ε = 2h and κ = 10 -10 h. Biot's coefficient is α = 1. Furthermore the wellbore pressure is q F = 5 m 3 /s and M B = 1 × 10 8 Pa, c F = 1 × 10 -8 Pa. The viscosities are chosen as ν R = ν F = 1 × 10 -3 Ns/m 2 . The reservoir permeability is K R = 1 d. The densities are ρ 0 F = 1 kg/m 3 . The critical energy release rate is chosen as G c = 1 × 10 3 N m -1 , Young's modulus is E = 10 8 Pa, and Poisson ratio is set to ν = 0.2. The relationship to the Lamé coefficients µ and λ is given by:

µ = E 2(1 + ν) , λ = νE (1 + ν)(1 -2ν) .
Observation goals. We observe the length of the fracture in terms of the phase-field variable, the predictor-corrector mesh refinement, the crack-width distribution, and the pressure and saturation distributions.

Discussion of our findings. In Figure 7.1, we observe the crack pattern and the predictorcorrector mesh refinement at three different times. The fracture width solution is displayed in Figure 7.2. Here we observe a smooth width field, which is parabolic with respect to the injection point. The largest width can be identified in the injection point. The pressure and saturation distributions are shown in Figure 7.3. Pressure is similar to our previous results. Namely, we observe negative pressure at the fracture tips and pressure dropping when fracture starts propagating (see Figure 7.3). The saturation develops a slower than a priori expected.

On the other hand, the material is fully saturated and the injecting fluid needs to push the existing fluid away, which requires a certain time. We ran a second series of tests, referenced as Example 1b. A common claim in phase-field studies is that the fracture is independent of the mesh and does not to be aligned with mesh faces. In order to study this independence we shifted the initial fracture by 45 degrees. The results are displayed in the Figures 7.4 A quantitative comparison of the maximum pressure and the maximum crack width (crack opening displacement -COD) is provided in Figure 7.7. While there are slight differences in the pressures the width evolutions show excellent agreement for the straight and diagonal crack settings, respectively. We thus conclude that the model is independent with respect to the initial crack location; here a rotation. The crack width values also slightly differ, but much less. From these data, we conclude a very good agreement and thus an independence of the model with respect to the geometry (i.e., here a rotation of the initial fracture). We expect the data to converge closer on finer meshes, which we may further quantify in future studies.

With regard to the computational cost of the nonlinear solver for the displacement/phasefield system (details can be found in [START_REF] Heister | A primal-dual active set method and predictorcorrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF]) and fixed-stress iterations (see Figure 5.1 and Algorithm 1), we observe very similar iteration numbers in this example as recently published in [START_REF] Lee | Iterative coupling of flow, geomechanics and adaptive phase-field fracture including level-set crack width approaches[END_REF][Section 6.3]. Specifically we have in average about 5 Newton iterations per time step for solving Formulation 5. Fixed-stress iteration numbers are displayed in Figure 7.8, which are also very reasonable.

Configuration.

Here Ω := (0 m, 4 m) 2 and we prescribe three parallel initial cracks with half lengths l 0 = 0.2 m as illustrated in Figure 7.9. The initial mesh is locally refined an additional four times around the fracture than the reservoir (see Figure 7.9b).

Boundary and initial conditions. For the displacements we prescribe homogeneous Dirichlet conditions on ∂Λ. For the phase-field and the pressure system, homogeneous Neumann conditions are employed. The boundary conditions for the crack-width problem have been explained in Section 5.4. The initial fractures are positioned at (β -h min , β + h min ) × (2l 0 , 2 + l 0 ), where β = 0.75 m, 2 m, and 3.25 m with the phase value set to zero. The initial displacement and pressure values are also set to zero.

Parameters. The smallest mesh size is h min = 0.0055 m and the regularization parameters for phase-field are chosen as ε = 2h min and κ = 10 -10 h min . In the homogeneous media, the critical energy release rate is chosen as G c = 10 N m -1 , Young's modulus is E = 10 8 Pa, and Poisson ratio is set to ν = 0.2. For the fluid parameters, Biot's coefficient is set as α = 1, M B = 1 × 10 8 Pa, η R = η F,res = 3 × 10 -3 Ns/m 2 , η F,inj = 1 × 10 -3 Ns/m 2 , the reservoir permeability is K R = 1 D ≈ 1 × 10 -12 m 2 , q F = 10, and c F = 1 × 10 -12 Pa -1 . The capillary pressure coefficients are γ res = 0.01, θ = 0, φ = 1, for (3.50) and to avoid the singularity we use p c (S ε ) where S ε = 0.2. The numerical parameters are set to α p = 100, α s = 100, and c s = 1. The uniform time step size is δt = 0.01.

Observation goals. First, we observe the length of the fracture in terms of the phasefield variable and expect the stress shadowed effect [START_REF] Bai | Explanation for fracture spacing in layered materials[END_REF][START_REF] Lachenbruch | Depth and spacing of tension cracks[END_REF][START_REF] Olson | Predicting fracture swarms -the influence of subcritical crack growth and the crack-tip process zone on joint spacing in rock[END_REF] with close enough parallel fractures. Then, we plot the saturation distributions inside the fractures. Discussion of our findings. Figure 7.10 illustrates the fracture propagation for different time steps with corresponding phase-field values. Injected fluid saturation S F,inj is presented in Figure 7.11. Stress shadowed effect is observed where the left and right fractures are curving out. The optimal distances between fractures to maximize propagated fracture surfaces and to avoid these effects with the phase-field model are discussed in [START_REF] Lee | Optimal design of hydraulic fracturing in porous media using the phase field fracture model coupled with genetic algorithm[END_REF].

Example 3:

A propagating penny-shaped fracture in a three dimensional homogenous medium. In this final example, a penny shaped fracture is propagating in a three dimensional domain. The idea of this setting is similar to Example 1 and most parameters are the same. From a coding perspective the extension from 2D to 3D is a major step, but relatively easy to accomplish in our code.

Configuration. The domain is a cube Ω = (0 m, 4 m) 3 and a (prescribed) initial pennyshaped crack with radius l 0 = 0.25 m. The initial mesh is 2 times uniformly refined and we allow for additional 5 levels of predictor-corrector mesh refinement. The boundary conditions are the same as in Example 1, but now in 3D. The initial conditions similar to the 2D test cases: the initial phase-field value is 0 in the fracture and 1 outside. The remaining variables are initialized by zero values.

Parameters. The time step size is chosen as δt = 0.01 s. The total number of time steps is 250 with the final time 2.5 s. The phase-field value for the predictor-corrector refinement is chosen as 0.8; thus below 0.8 the mesh will be refined up to three times. The penalization parameter for the enriched Galerkin pressure system is set to 100. The penalization parameter for the saturation system is set to 1000. The capillary coefficient is γ res = 1. Furthermore, T OL F S = 10 -3 and the Newton tolerance for the phase-field displacement system is 10 -10 .

The regularization parameters are chosen as ε = 2h and κ = 10 -10 h. Biot's coefficient is α = 1. Furthermore the wellbore pressure is q F = 2 m 3 /s and M B = 1 × 10 8 Pa, c F = 1 × 10 -8 Pa. The viscosities are chosen as ν R = ν F = 1 × 10 -3 Ns/m 2 . The reservoir permeability is K R = 1 d. The densities are ρ 0 F = 1 kg/m 3 . The critical energy release rate is chosen as G c = 1 × 10 3 N m -1 , Young's modulus is E = 10 8 Pa, and Poisson ratio is set to ν = 0.2.

Observation goals. We observe the length of the fracture in terms of the phase-field variable, the predictor-corrector mesh refinement, the crack-width distribution, and the pressure and saturation distributions.

Discussion of our findings. Our findings are displayed in the Figures 7.12,7.13, and 7.14. We observe that the fracture propagates in a penny-shaped fashion and the mesh is locally refined with the predictor-corrector scheme. In particular in these 3D this procedure helps enormously to reduce the computational cost. In more detail, we start with a total number of 7792 at T = 0.01s, over 8184 at T = 1.25s, and ending with 23724 grid cells at T = 2.5s. The saturation (of the injecting fluid) propagates slower than the fracture itself, which we already observed for the 2D test cases. Our simulations confirm our modeling parts in which we assumed a radial injecting fluid.

Conclusion.

In this paper, we proposed a model for immiscible two phase flow in a propagating fracture using a variational phase-field model. A residing fluid can be considered as oil whereas the injecting fluid may be considered to be water.

The flow model in the fracture is derived from multifluid Stokes equations, employing upscaling. The calculation is a natural extension of our previous work [START_REF] Lee | Phase-field modeling of proppantfilled fractures in a poroelastic medium[END_REF] and [START_REF] Mikelić | A phase-field method for propagating fluidfilled fractures coupled to a surrounding porous medium[END_REF], where the lubrication approximation was used to calculate the permeability tensor for a single phase fluid flow through a penny-shaped fracture. In this paper, we have considered a thin pennyshaped fracture filled by two viscous incompressible fluids, separated by a free boundary. The wetting fluid was injected radially, while displacing the existing residing fluid. Using the approach from [START_REF] Mikelić | On the derivation of the Buckley-Leverett model from the two fluid Navier-Stokes equations in a thin domain computational geosciences[END_REF], the free boundary was described using a hyperbolic equation for the viscosity. Thus, the free boundary defines the saturation. We upscaled the multifluid Stokes equations to the Buckley-Leverett equations and calculated the absolute and relative permeabilities as functions of the saturations explicitly by using a two-scale expansion.

Moreover, we also account for capillary pressure. These quantities are then used in three dimensional equations for pressure and saturation in the fracture. Flow in the poroelastic medium is assumed to be of single-phase. In order to obtain an accurate fracture width (i.e., jump of displacements across the fracture boundaries), we adopt a level-set approach and a finite element fracture-width problem. The next challenging part is the coupling of the two phase flow system with single-phase flow on the fracture boundary. Here we employ diffraction systems with varying coefficients in space and time. These developments result in a four-field problem: pressure, saturation, phase-field-displacements, and fracture-width.

In order to ensure conservation properties of the continuous system, we develop specifically tailored discretization schemes. For the phase-field-displacement and crack-width systems, standard continuous finite element method are employed. The numerical discretization of pressure and saturation ask for local mass conservation. These are ensured by using an enriched Galerkin scheme. In order to keep the computational cost reasonable, we implemented local mesh adaptivity that allows for small phase-field regularization parameters in the fracture zone and an MPI-parallelized framework.

While the pressure, saturation and fracture-width systems are linear and can be solved directly with iterative solvers, the phase-field displacement system is nonlinear and respects a crack irreversibility condition. Here, we use a combined semi-smooth Newton solver based on a primal-dual active set strategy. In the final section, we present three prototype test cases that show the capabilities of the presented modeling and computational frameworks.

While this is a first study with detailed modeling and specifically-tailored construction of numerical discretization schemes, future work will be focusing on further treating more practical field problems.

  Furthermore, Ω F (t) is an open and connected set corresponding to the fracture and Ω R (t) is an open and connected set corresponding to the unbroken material, for every t ∈ [0, T ]. The phase-field is a smoothed indicator function of Ω R (t). Here, Λ = Ω F ∪ Ω R ∪ Γ F is the entire domain and Γ F = Ω F ∩ Ω R is the fracture boundary. The boundary ∂Λ splits into a Dirichlet part ∂ D Λ and a Neumann part ∂ N Λ.

(3. 26 )

 26 Employing the kinematic condition (3.15) with (3.23) yields:

  .40) Inserting (3.40) into (3.39) yields

REMARK 3 . 6 .

 36 Due to the relations (3.35) and (3.36), it is sufficient to solve only one of the saturation equations (3.51) -(3.52) and analogously only the pressure equation (3.53). To concentrate on one of the respective equations is done in Section 3.1.10.

Fig. 3 . 2 :

 32 Fig. 3.2: Two relative permeabilities and a capillary pressure are illustrated for the case M = 1. For simplicity, p c (s) = (k r1 (s)) -1/2 + (k r2 (s)) -1/2 (1 -s) is plotted by setting all the coefficients simply to one in (3.50).

Fig. 3 . 3 :

 33 Fig. 3.3: Illustration of different domains and variables. The middle void area is a fracture zone, Ω F (t). Here, p R (t) is the pressure value in the reservoir domain Ω R (t) and s F,inj and p F,inj are injecting fluid saturation and pressure in a fracture, respectively.

[

  [w]] = w + n + + w -n -on e ∈ E o h , and [[w]] = wn for e ∈ E ∂ h .

Fig. 5 . 1 :

 51 Fig. 5.1: The global algorithm flowchart. During the fixed-stress iteration, an iterative Implicit Pressure Explicit Saturation formulation (IMPES) is embedded.

  (a) Diffusive domains Ω F , Ω R , and Ω ε and ε. (b) Sharp domains Ω F , Ω R , and Γ F described by Ω F , Ω R , and Ω ε .

Fig. 5 . 2 :

 52 Fig. 5.2: Definitions for different domains are depicted with a phase-field variable. (a) Ω F , Ω R , and Ω ε , are considered as diffusive domains and defined where Φ = 0, Φ ∈ (0, 1) , and Φ = 1, respectively. Ω ε is the transition zone between Ω F and Ω R with the length ε. (b) Ω F , Ω R , and Γ F are considered as sharp domains and interface. Relations between the sharp and diffusive domains are illustrated, i.e Ω F := Ω F ∪ Ω ε,F and Ω R := Ω R ∪ Ω ε,R . The red color indicates the sharp fracture zone.

Fig. 7 .

 7 Fig. 7.1: Example 1a: fracture path displayed in terms of the phase-field variable on top. At the bottom the locally adapted mesh where the refinement criterion depends on the phase-field value at each cell. The visualizations are taken at T = 0.01s, 2s, 10s.

Fig. 7 . 3 :

 73 Fig. 7.3: Example 1a: the pressure and saturation distributions at T = 0.01s, 2s, 10s.

Fig. 7 . 4 :

 74 Fig. 7.4: Example 1b: fracture path displayed in terms of the phase-field variable on top. At the bottom the locally adapted mesh where the refinement criterion depends on the phase-field value at each cell. The visualizations are taken at T = 0.01s, 2s, 10s.

Fig. 7 . 5 :

 75 Fig. 7.5: Example 1b: the fracture width at T = 0.01s, 2s, 10s.

Fig. 7 . 6 :

 76 Fig. 7.6: Example 1b: the pressure and saturation distributions at T = 0.01s, 2s, 10s.

  ,7.5, and 7.6. Here we indeed observe the same range of values for the width, pressure and saturation as for the straight crack. The length of the fracture, displayed in Figure 7.4 corresponds to the straight crack length shown in Figure 7.1.

Fig. 7 . 8 :

 78 Fig. 7.8: Example 1: Fixed stress iteration numbers per time step for both fracture settings.

Fig. 7 . 9 :

 79 Fig. 7.9: Example 2: Setup

Fig. 7 . 10 :

 710 Fig. 7.10: Example 2: phase-field fracture values for different time step numbers n.

  (a) n = 10 (b) n = 50 (c) n = 100 (d) n = 150 (e) n = 175 (f) n = 225

Fig. 7 . 11 :

 711 Fig. 7.11: Example 2: Injected fluid saturation values for different time step numbers n.

Fig. 7 . 12 :

 712 Fig. 7.12: Example 3: the fracture in terms of the phase-field variable and the locally refined mesh at T = 0.01s, 1.25s, 2.5s.

Fig. 7 .

 7 Fig. 7.13: Example 3: the fracture width at T = 0.01s, 1.25s, 2.5s. The color scale is adjusted dynamically in each subfigure, i.e., rose/red colors display the highest width values at the current step. These maximal values are 3.185 × 10 -5 , 3.425 × 10 -4 and 4.16 × 10 -4 , respectively.

Fig. 7 . 14 :

 714 Fig. 7.14: Example 3: the pressure (first row) and saturation distributions (second row) at T = 0.01s, 1.25s, 2.5s. In the last row a zoom-in of the saturation values is provided.

  Quantitative comparison of the evolutions of the maximum pressure and the maximum crack width. At the highest pressures, the values differ by about 15% while the difference is about 3% at the end at T = 10s.

		700000				Ex. 1a			0.012
		600000				Ex. 1b			0.01
	Pressure [Pa]	200000 300000 400000 500000						Crack width [m]	0.004 0.006 0.008
		100000							0.002				Ex. 1a
		0							0				Ex. 1b
		0	2	4	6	8	10	12	0	2	4	6	8	10	12
					Time [s]							Time [s]
	Fig. 7.7: Example 1:

The reason why we ordered the subsections of Section 5 in that way, having the pressure/saturations last, is that for the pressure equation we need the width from the displacements. On the other hand, the classical idea of the fixed-stress algorithm is to solve first for flow (pressures) and afterwards for geomechanics.

Acknowledgments. The research by S. Lee and M.F. Wheeler was partially supported by Center for Frontiers of Subsurface Energy Security, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, and Office of Basic Energy Sciences, DOE Project #DE-SC0001114, and a DOE grant #DE FG02-04ER25617. M. F. Wheeler was also partially supported by Moncrief Grand Challenge Faculty Awards from The Institute for Computational Engineering and Sciences (ICES), the University of Texas at Austin. A. Mikelić and T. Wick were supported by The J. Tinsley Oden Faculty Fellowship Research Program. and Center for Subsurface Modeling at Institute for Computational Engineering and Science (ICES), the University of Texas at Austin. T. Wick thanks his former host institution Centre de Mathématiques Appliquées at École Polytechnique (France) for providing support and resources in preparing this paper.