Diffuse interface modeling and variationally consistent homogenization of fluid transport in fractured porous media

We critically assess diffuse interface models for fluid transport in fractured porous media. Such models, often called fracture phase field models, are commonly used to simulate hydraulic stimulation or hydraulic fracturing of fluid-saturated porous rock. In this paper, we focus on the less complex case of fluid transport in stationary fracture networks that is triggered by a hydro-mechanical interaction of the fluid in the fractures with a surrounding poroelastic matrix material. In other words, fracture propagation is not taken into account. This allows us to validate the diffuse interface model quantitatively and to benchmark it against solutions obtained from sharp interface formulations and analytical solutions. We introduce the relevant equations for the sharp and diffuse, i.e. fracture phase field, interface formulations. Moreover, we derive the scale-transition rules for upscaling the fluid-transport problem towards a viscoelastic substitute model via Variationally Consistent Computational Homogenization. This allows us to measure the attenuation associated with fluid transport on the sub-scale. From the numerical investigations we conclude that the conventional diffuse interface formulation fails in predicting the fluid-transport behavior appropriately. The results even tend to be non-physical under certain conditions. We, therefore, propose a modification of the interpolation functions used in the diffuse interface model that leads to reasonable results and to a good approximation of the reference solutions.


Introduction
In recent decades, hydraulic stimulation, often called hydraulic fracturing in the context of shale gas production, has become a promising technology for the production of deep geothermal energy in a carbon neutral economy. In the past, however, stimulation of deep geothermal reservoirs induced earthquakes, as observed in e.g. the Rhine valley (Charl� ety et al., 2007). For the societal acceptance of this technology it is, therefore, of utmost importance to assess the risk of hydraulic stimulation in a reliable fashion, in particular for its application in densely populated regions.
Recently, manifold approaches have been launched that aim at modeling and simulating hydraulic stimulation using computational mechanics. Among others, the family of fracture phase field models has become most popular. Seminal early contributions are those of Miehe et al. (Miehe and Mauthe, 2016;Miehe et al., 2010Miehe et al., , 2015 and Mikelic et al. (Mikeli� c et al., 2015a, 2015b who proposed that the damage variable interpolates smoothly between the fully disintegrated state at the fracture and an undamaged poroelastic background. In practice, the fractures contribute to a Darcy-type fluid transport by an additional permeability term whilst, in turn, the fluid pressure at the fracture tip triggers the fracture propagation. These concepts were extended towards fracture propagation in strongly heterogeneous materials (Xia et al., 2017) and partially saturated porous media (Cajuhi and Lorenzis, 2018). Moreover, Stokes-type transport was included (Ehlers and Luo, 2017;Wilson and Landis, 2016), and the fracture phase field formulation was embedded in the Theory of Porous Media (Ehlers and Luo, 2017;Heider and Markert, 2017).
In this paper, we aim at investigating how well fluid transport in fracture networks embedded in fluid-saturated porous media can be predicted using diffuse interface formulations such as the fracture phase field model. However, quantitative benchmarks for the problem of hydraulic stimulation are difficult to assess, if existent at all. We, therefore, focus on the sub-problem of hydro-mechanically triggered fluidtransport in stationary fracture networks. Hence, fracture propagation is not taken into account. This allows us to validate hydro-mechanically coupled fluid transport in fracture networks, obtained using a diffuse interface formulation, against suitable numerical and analytical benchmark experiments.
This scenario of fluid transport in stationary fracture networks is closely related to seismic attenuation which can be observed if seismic waves travel through a fluid-saturated fractured rock. If mechanical waves, usually at seismic (i.e. very low) frequencies, travel through such a medium, the heterogeneity of the fracture network leads to strong fluid pressure gradients. The pressure gradients are equilibrated via redistribution of fluid in the fracture network. This process is dissipative and governed by pressure diffusion along the fractures (O'Connell and Budiansky, 1977;Vinci et al., 2014). Hence, part of the wave energy is lost and the seismic wave is attenuated. Moreover, it was shown that the effect of seismic attenuation correlates with the degree of fracture connectivity in the network, i.e. with the amount of fractures that are connected and which communicate hydraulically (Quintal et al., 2014(Quintal et al., , 2016Rubino et al., 2013). On the much larger length scale of seismic wave propagation, the local fluid redistribution in the fracture network can be interpreted as viscoelastic damping. The correlation between pressure diffusion in fractures and the properties of an overall viscoelastic medium were investigated in (J€ anicke et al., 2019a, 2019b) in terms of Variationally consistent Computational Homogenization (VCH) together with Numerical Model Reduction (NMR). The motivation to focus on the sub-problem of seismic attenuation associated with fluid transport in fracture networks is that phenomenon is of high relevance in the context of hydraulic stimulation: Fracture networks increase the hydraulic conductivity of a rock tremendously (Odling, 1992;Snow, 1969). Therefore, the ability to detect, to simulate, and to understand seismic attenuation enables practitioners to decide whether or not a rock formation is suitable for hydraulic stimulation. Alternatively, during hydraulic stimulation, a tool is provided for monitoring the fracturing process, i.e. the generation of induced fracture networks. Another reason to focus on that particular sub-problem is that there are well-established formulations available which treat the fractures as sharp interfaces (e.g. using interface elements) and which can be used to benchmark the phase field approach quantitatively. Hence, the material response associated with diffusive fluid transport in sub-scale fracture networks in sharp or diffuse formulation can be directly used for FE 2 type simulations of the seismic wave propagation experiment.
The paper is organized as follows: We first recall the sharp interface formulation of hydro-mechanically coupled fluid transport in fracture networks in Section 2. Next, we introduce the upscaling procedure for selective homogenization from a poroelastic Representative Volume Element (RVE) with fluid-filled fractures, represented as sharp interfaces, towards an apparently viscoelastic substitute model via VCH in Section 3. In Sections 4 and 5, we derive the simplified phase field model for stationary fractures with poroelasticity and the corresponding selective VCH approach. Finally, we investigate the performance of the diffuse interface model compared to the sharp interface reference and an analytical solution in Section 6 and discuss our findings in Section 7.
Throughout this manuscript, vector and tensor quantities are written as bold types. Taking into account Einstein's sum convention we write out simple and double contractions in terms of the Cartesian components as a⋅b ¼ a i b i and A : B ¼ A ij B ij . Time derivatives are denoted _ ⋄.

Preliminaries
In this section, we recall a sharp interface formulation for pressure diffusion in a fractured poroelastic medium that was initially presented in (J€ anicke et al., 2019a). We refer e.g. to (Coussy, 2004;Ehlers et al., 2002;Verruijt, 2013) for a more detailed overview on the modeling of (undamaged) poroelastic media. The fracture opening is assumed to be on the length scale of pores and grains building up the poroelastic background material. In other words, the typical fracture length is orders of magnitude larger than the fracture opening. Thus, the fractures can be approximated as lower-dimensional interfaces with the fracture opening τ. Moreover, we assume the fractures to be mechanically and hydraulically open. Hence, fracture surfaces are not in contact, such that pressure diffusion, associated with fluid transport, takes place along the fracture. The fracture network is stationary, and we neglect fracture propagation. The primary global fields for the domain of the poroelastic background medium are the displacement of the solid skeleton u and the pore fluid pressure p M . The fluid pressure in the fracture domain is denoted p F .

Strong and weak form of the fine-scale problem
The general formulation of the fracture interface model can be stated as follows: Restricting to quasi-statics, we seek the displacement field uðx; tÞ : Ω M � R þ →R 3 , the pore pressure field in the poroelastic domain p M ðx; tÞ : Ω M � R þ →R, and the fluid pressure in the fracture domain p F ðx; tÞ : Γ F � R þ →R that solve the system (1) together with the initial conditions Here, the fracture network Γ F is assumed to consist of i ¼ 1; 2; …; n F possibly intersecting fractures such that Γ F ¼ [ nF i¼1 Γ F i . Γ M is the exterior boundary of Ω M whilst L F i is the exterior boundary associated with Γ F i , see Fig. 1. We introduce for simplicity the gradient operator r F ¼ ½I À n F � n F �⋅rthat is tangential to the fracture, where n F represents the vector normal to the fracture, contrary to the normal vector of the domain n. Further, Φ M and Φ F are storage functions that are defined subsequently. ⋄ p denotes a prescribed quantity (Dirichlet/Neumann boundary condition). (1a) represents static equilibrium whilst (1e) and (1i) express mass conservation of the fluid saturating the poroelastic material and the fractures. The constitutive relations for the equilibrium stress σ and the seepage velocities w M and w F become where ε is the linear symmetric strain, and square brackets �½⋄� refer to operational dependency of � on ⋄. E is the fourth order elasticity tensor of the undamaged material which is, in the simplest case of isotropy, defined by the shear modulus G and the bulk modulus K of the dry solid background skeleton. K M and K F are the second order permeability tensors of the poroelastic material and the fractures and α is the Biot-Willis coefficient. In the simplest case of isotropy, the permeability tensors can be expressed as with the permeability k of the poroelastic material and the viscosity η of the saturating fluid. The permeability of the fracture zone is computed from the Poiseuille flow rule, i.e. we assume a quadratic velocity profile with non-slip conditions in the fractures, which results in the apparent fracture permeability τ 2 0 =12 where τ 0 ¼ τðt ¼ 0Þ is the given initial fracture opening. Here, we used the pressure gradient operator ζ½p M � for the background medium and the tangential gradient operator ζ F ½p F � that is evaluated in the plane tangential to the fracture domain with the normal vector n F . Furthermore, q L denotes the leak-off of fluid from the fracture into the surrounding poroelastic background material whereas b q i denotes the fluid exchange between intersecting fractures. Moreover, we introduce the storage functions Φ M and Φ F that represent the current volume fraction of fluid in the background material and the fracture as The Biot coefficients for the poroelastic material are defined as where φ is the porosity of the poroelastic material, K f the bulk modulus of the fluid saturating pore space and fractures, and K s the bulk modulus of the solid grains constituting the poroelastic background material. The volume change of the fracture due to deformation of the background material is denoted e F . We specify the hydro-mechanical coupling between the poroelastic background material in Ω M and the fluid-filled fractures Γ Fi . We, therefore, define the volumetric strain in the fracture due to the deformation of the surrounding poroelastic material and the leak-off of fluid from the fracture into the poroelastic material as where k � k F is the jump of a quantity across the fracture. Finally, we assume the pressure field to be continuous such that Henceforth, we simplify notation and write p instead of p M and p F . With these definitions at hand, we introduce the standard spacevariational format corresponding to (1) as follows: Find uðx; tÞ; pðx; tÞ in the appropriately defined spaces U � P that solve Z ΩM ½ε½u� : E : ε½δu� À αpI : where U 0 ; P 0 are appropriately defined test spaces.

VCH for the sharp interface model
In this section, we aim at establishing a scheme for the upscaling from the fluid-filled fractured rock towards a viscoelastic substitute medium. For the sharp interface model, a similar procedure is described in (J€ anicke et al., 2019a). Therefore, we first introduce the general concept of first-order homogenization. Afterwards, we derive scale transition rules that constitute the macro-scale problem and establish the RVE problem in a variationally consistent way.

First-order homogenization in the spatial domain
We replace the single-scale problem in Section 2 by a two-scale problem upon introducing (i) running averages in the weak format and (ii) scale separation via first-order homogenization. The first step is to replace the space-variational problem (9) by that of finding ðuðx; tÞ; pðx; tÞÞ 2 U FE 2 � P FE 2 that solves Z ΩM h a ðuÞ □;M ðu; δuÞ À b □;M ðp; δuÞ À b □;F ðp; δuÞ where the pertinent space-variational forms in (10) are given as The volume averaging operators are defined as These forms represent running averages on square (2D) or cubic (3D) domains Ω □ with the side length L □ which are located at each macroscale spatial point x 2 Ω. We note that the introduced averaging results in a change of the original problem. The weak format (10) thus serves as the basis for VCH as discussed below.
Next, we introduce scale separation and first-order homogenization. We decompose the sub-scale field u into a macro-scale part u M and a micro-scale or fluctuation part u s within each RVE such that Furthermore, we introduce the simplification that pressure diffusion associated with fluid transport in the pores and fractures is a local process within the RVE. Hence, we presume the absence of macroscopic pressure gradients ζ½p� and consider the RVE to be undrained, i.e. the amount of fluid stored in the RVE remains constant. We thus choose This upscaling procedure with undrained conditions for poroelastic RVEs, sometimes called "selective homogenization", has been discussed in more detail in (J€ anicke et al., 2016, 2019a).
We define the two-scale trial and test spaces as where U and U s □;j are the finite element discretized trial spaces for the macro-and microscopic problem on RVE Ω □;j , respectively. Similarly, U 0 is the trial space for the macroscopic problem and P s □;j is the trial and test space for RVE Ω □;j .
Hence, we introduce the trial and test spaces for the macro-scale problem U and U 0 and for the micro-scale problem U s □;j ; P s □;j , where we use the notation In other words, there exist components ðu s ; p s Þ and corresponding spaces ðU s □;j ; P s □;j Þ for each individual RVE Ω □;j . We are now in the position to restate the two-scale problem as follows: Find ðu; fu s

Macro-scale problem and stress averaging
We obtain the homogenized macro-scale problem from (17) upon considering only macroscopic contributions to the test functions. Hence, δp s ¼ 0 and δu s ¼ 0. We find that the macro-scale problem reduces to Z where we assume that t p is the suitably homogenized traction vector on the Neumann boundary. Moreover, we find the relation With the standard argumentation, the stress averaging rule (19) can be expressed as the surface integral

Micro-scale problem on an RVE
In order to generate a uniquely solvable RVE problem we adopt the additional model assumption that u s and p s are micro-periodic, i.e.
We, therefore, use the difference operator k � k □ ðxÞ ¼ �ðx þ Þ À �ðx À Þ. We call x 2 Γ þ □ image points and x 2 Γ À □ the corresponding mirror points. We apply micro-periodicity in a variational form such that d with the space L 2 of square integrable functions. The space-variational RVE format can now be stated as follows: For a given history εðtÞ, find uðx; tÞ 2 U □ , pðx; tÞ 2 P □ ; λ ðuÞ 2 T ðuÞ □ , and λ ðpÞ ¼ ðλ Finally, we skipped the index ð�Þ M for the RVE formats in (25a) and (25c) for simplicity reasons.
Remark 1. Setting δλ ðuÞ ¼ δσ⋅n with σ 2 R 3�3 sym into (25) and considering the balance of momentum of the fluid phase in the fracture network, one can show that 〈ε〉 □ ¼ ε: (27) Hence, the volume averaging rule for the strain tensor ε is built-in the problem formulation. (25a) and (25b) and taking into account (20), we are able to derive the Hill-Mandel macrohomogeneity condition as 4. Diffuse interface model for fractured poroelastic media

Preliminaries
In this section, we develop a model for pressure diffusion in a damaged poroelastic medium under geometrically linear and quasistatic conditions. The damage field is constructed in a way that it represents networks of fluid-saturated fractures in a diffuse fashion. Hence, the fractures are approximated in the spirit of a fracture phase field formulation. As stated earlier, the purpose of our study is to investigate how well fluid pressure diffusion induced by hydro-mechanical coupling in a fractured poroelastic medium can be approximated by a diffuse fracture formulation. Therefore, our starting point is an existing fracture network whilst fracture propagation is neglected. This simplification results in a partly decoupled equation system. The primary global fields are the damage variable d, the (continuous) displacement of the solid skeleton u, and the pore fluid pressure p.

Strong and weak format of the fine-scale problem
The general formulation of a fracture phase field in poroelastic media can be stated as follows: Restricting to quasi-statics, we seek the damage field dðx; tÞ : Ω � R þ →½0; 1�, the displacement field uðx; tÞ : Ω � R þ →R 3 , and the pore pressure pðx; tÞ : Ω � R þ →R that solve the system (29) Here, (29e) and (29h) express static equilibrium and the continuity equation of the two-phase solid-fluid mixture. In order to study stationary fracture networks, we decouple the damage conservation law (29a) from static equilibrium and the continuity equation and introduce simplified constitutive relations for the damage flux f, the damage source term Θ and the viscous regularization D such that Here, γ is a length scale parameter representing the smoothing length of the diffuse fracture representation. It is important to remark that the smoothing length γ has no physical meaning but is introduced to regularize the problem. 1 ζ½d� is the damage gradient. Obviously, the particular choice in (30) leads to a decoupling of (29a) from (29e) and (29h) and to a stationary damage field. This stationary solution can be found independently of static equilibrium and continuity equation in a pre-computation step. The resulting damage distribution dðxÞ can be treated as an input parameter to (29e) and (29h).
The constitutive relations for the equilibrium stress σ and the seepage velocity w become σðε½u�; pÞ : ¼ aðdÞE : ε À αðdÞpI; ε½u� : ¼ ðu � rÞ sym ; (31a) wðζ½p�Þ : ¼ À KðdÞ ⋅ ζ½p�; ζ½p� : ¼ rp: Note that d is now a parameter, not a solution field. The decay of the elastic properties associated with the damage field d is built-in via the degradation function aðdÞ which will be specified further below. In the case of isotropy, the second order permeability tensor KðdÞ can be expressed in standard fashion as with the permeability kðdÞ and the viscosity η of the saturating fluid.
Finally, the Biot-Willis coefficient is defined in accordance with (31a) using the degradation function aðdÞ such that It remains to introduce the storage function ΦðdÞ that represents the current volume fraction of fluid in the damaged poroelastic medium as Φðd; ε½u�; pÞ ¼ φðdÞ þ αðdÞI : ε½u� þ βðdÞp; where φðdÞ is the porosity, and βðdÞ is the intrinsic compression compliance of the pore fluid which is, in the case of the diffuse fracture model, computed as Altogether, it is important to remark that the material parameters of the fractured poroelastic medium are scaled between the undamaged poroelastic state and the fully disintegrated state at the fracture. The scaling is built-in via (i) the degradation function aðdÞ applied to the overall elastic properties of the solid skeleton ðE; K; GÞ, (ii) the porosity φðdÞ, and (iii) the permeability kðdÞ. Those functions will be defined further below.
With these definitions at hand, we introduce the standard spacevariational format corresponding to (29a) -(29d) as follows: Find dðxÞ in the appropriately defined space D that solves Z And, from (29e) -(29j), find uðx; tÞ; pðx; tÞ in the appropriately defined spaces U � P that solve where D 0 , U 0 and P 0 are the appropriately defined test spaces.

Degradation function
It remains to specify the functions aðdÞ, φðdÞ, and kðdÞ that interpolate between the undamaged poroelastic medium and the fracture zone. We follow (Miehe et al., 2010;Mikeli� c et al., 2015a) and introduce the degradation function as where the small parameter κ is used for numerical stability. Hence, the degradation function interpolates between the undamaged state aðd ¼ 0Þ ¼ 1 and the fully disintegrated stated aðd ¼ 1Þ ¼ κ. This choice for aðdðxÞÞ guarantees C 1 -continuity in x.
In our numerical studies, we will investigate the degradation function on the basis of an alternative damage parametrization d * . For this purpose, we refer to the analytical solution of one single fracture with d ¼ 1 at the interface Γ F in an infinite domain Ω ∞ . Under these circumstances, the solution of the simplified Euler equation (29a) can be computed as see Fig. 2. We use this analytical result to define a threshold value and the modified damage parametrization As visualized in Fig. 2 b), the alternative damage function d * is constant within a band of �γ * along the fracture. This feature will be discussed in more detail below. For now, we note that the degradation function can be reformulated in terms of d * such that

Interpolation of porosity and permeability
The porosity φðdÞ represents the volume fraction of fluid stored in the pore space of the poroelastic medium as well as in the damage zone. In other words, we aim for defining φðdÞ such that, in an integral sense, the fractured poroelastic medium in a diffuse interface representation contains the same amount of fluid as in a sharp interface representation of the fracture network. It is of utmost importance to evaluate the amount of fluid stored in the fracture network consistently since, otherwise, the storage capacity of the fractures would be substantially overestimated by a diffuse interface formulation. Based on the analytical solution (39), we compute the apparent fracture porosity b from the relation 1 Note that physical transition zones exist at the interface between fluid domains and porous media, cf. (Beavers and Joseph, 1967). Such effects are neglected in the present approach assuming that the length of the "true" transition zone is much smaller than γ.
Here, we use the (true) fracture porosity φ D ¼ 1, i.e. the fluid volume fraction of the resolved fracture, and the crack opening τ of the corresponding sharp interface model. Combining (39) and (43) Hence, the apparent porosity of the fracture zone depends on the smoothing length ðγ; γ * Þ and represents nothing else than the ratio of the true thickness of the fracture and the width of the diffuse damage zone. In other words, the larger the smoothing length is chosen the smaller the apparent fracture porosity becomes.
The exponent ν in (45) allows us to interpolate between the arithmetic mean ν ¼ 1 and the harmonic mean ν ¼ À 1. The interpolation function is visualized in Fig. 3 for two different porosity ratios φ s b . It can be seen that the exponent ν allows to tune the properties of the transition zone.
Finally, we use the same mixture rule to interpolate the permeability kðdÞ between the fully disintegrated fracture zone, where the permeability is computed as k F ¼ ðk s Þ aðdÞ À k F � 1À aðdÞ ; ν ¼ 0: (46b)

VCH for the diffuse interface model
Subsequently, we establish the two-scale formulation that homogenizes the fractured rock problem towards a viscoelastic substitute medium in analogy to the procedure presented in Section 3. Here, the

First-order homogenization in the spatial domain
We replace the single-scale problem in Section 4 by a two-scale problem upon introducing (i) running averages in the weak format and (ii) scale separation via first-order homogenization. As presented in Section 3, the first step is to replace the space-variational problems (36) and (37) where the pertinent space-variational forms in (47) Note again, that (47a) is decoupled from fu; pg whilst d serves as parameter in (47b) and (47c) instead of a solution field.
The forms introduced in (48) represent running averages on square (2D) or cubic (3D) domains Ω □ with the side length L □ which are located at each macro-scale spatial point x 2 Ω. The weak format (47) serves as the basis for variationally consistent homogenization as discussed below. Furthermore, we note that the representation in (47) still implies that the stationary damage problem (47a) can be treated independently from that of solving for ðuðx; tÞ; pðx; tÞÞ.
Next, we introduce scale separation and first-order homogenization. We decompose the sub-scale field u into a macro-scale part u M and a fluctuation part u s within each RVE such that u ¼ u M ½u� þ u s ; with u M ðx; tÞ : ¼ εðtÞ ⋅ ðx À xÞ; ε : ¼ ε½u�: (49) Moreover, we introduce the simplification that the stationary damage variable d representing the fracture network is a local property of the RVE. Hence, the damage field is not controlled by any macroscopic damage d or its gradient ζ½d�, and we impose where the fluctuation field d s will be specified later. Similarly, we consider pressure diffusion to be a local process within the RVE as discussed in Section 3. Hence, we presume the absence of macroscopic pressure gradients ζ½p� and consider the RVE to be undrained, i.e. the amount of fluid stored in the RVE remains constant. We write We define the two-scale trial and test spaces as Here, we introduced the trial and test spaces for the macro-scale problem U and U 0 and for the micro-scale problem D s □;j ; U s □;j ; P s □;j and C s □;j where we used the notation In other words, there exist components ðd s ; u s ; p s Þ and corresponding spaces ðD s □;j ; U s □;j ; P s □;j Þ for each singe RVE Ω □;i . Moreover, D s □;j and C s □;j are constructed in a way that they contain the conditions that d ¼ 1and We are now in the position to restate the two-scale problem (47) as follows: Find ðfd s j g; u; fu s j g; fp s j gÞ 2 ½� j D s 8ðfδd s g; δu; fδu s g; fδp s gÞ 2 ½� j C s □;j � � U 0 � ½� j U s □;j � � ½� j P s □;j � where we implicitly used the composition of (d; u; p) from (53).

Macro-scale problem and stress averaging
We obtain the homogenized macro-scale problem from (54) upon considering only macroscopic contributions to the test functions. Hence, δd s ¼ 0, δu s ¼ 0, and δp s ¼ 0. We find that the macro-scale problem reduces to Z where we assume that t p is the suitably homogenized traction vector on the Neumann boundary. We find the well-established relation Alternatively, the stress averaging rule (56) can be expressed as the surface integral σ ¼ 1 jΩ □ j Z Γ□ ðt � ½x À x�Þ sym dΓ: (57)

Micro-scale problem on an RVE
In order to generate a uniquely solvable RVE problem we adopt the additional model assumption that d s , u s , and p s are periodic, i.e.
Here, we use the difference operator k � k □ ðxÞ ¼ �ðx þ Þ À �ðx À Þ as introduced in Section 3. We apply micro-periodicity in a variational form such that The Lagrange multiplier λ ðuÞ represents a field of boundary tractions on Γ þ □ whereas λ ðdÞ and λ ðpÞ represent boundary fluxes on Γ þ □ . Moreover, we introduce the test spaces for boundary tractions and fluxes as with the space L 2 of square integrable functions. The space-variational RVE problem can now be stated as follows: Find d 2 D □ and λ ðdÞ 2 T For a given history εðtÞ, use the damage distribution dðxÞ from (62) and find uð �; tÞ 2 U □ , pð �; tÞ 2 P □ , λ ðuÞ ð �; tÞ 2 T  (63) we obtain the condition 〈ε〉 □ ¼ ε: (64) Hence, the volume averaging rule for the strain tensor ε is built-in the problem formulation.

Remark 2.
Setting δu ¼ _ u, δp ¼ p, summing up (63a) and (63b) and taking into account (57), we are able to derive the Hill-Mandel macrohomogeneity condition as a property built-in the problem formulation for any damage distribution dðxÞ computed from (62).

Numerical examples
After having completed the fine-scale formulations as well as the upscaling principles via VCH, we are now in the position to investigate the properties of the diffuse interface model and compare it with the sharp interface model as well as with an analytical solution. We investigate two benchmark experiments: (i) Outflux of fluid from a single drained fracture under compression and the resulting stress. (ii) Seismic attenuation in a prototype RVE problem consisting of two perpendicular and intersecting fractures.
In both cases, we examine how the solution of the hydromechanically coupled problem depends on the FE discretization of the diffuse fracture representation. In a first step, we choose the smoothing length γ and keep it constant whilst we vary the ratio γ=h e . Here, h e is the size of the finite elements in the vicinity of the interface. In other words, we investigate how many finite elements are at least needed to resolve the fracture zone in a mesh-independent fashion. In a second step, we study the convergence of the diffuse interface approximation towards the reference solution by reducing the smoothing length γ whilst the FE resolution in terms of the previously identified ratio γ=h e is kept constant. To simplify notation, the initial fracture aperture is denoted τ instead of τ 0 throughout the remainder of the paper.
For all numerical examples, we employ triangular Finite Elements with Lagrange shape functions by means of a Taylor-Hood formulation, i.e. quadratic in u, linear in p, quadratic in the (uncoupled) d.

Convergence study: single fracture in an impermeable rock matrix
We use the setup shown in Fig. 4 with L □ ¼ 2 m and τ ¼ 1e-5 m. The material properties are specified in Table 1. In this study, we aim at investigating the transport properties of the fracture under hydromechanical stimulation and measure the outflux via the left and right boundary and, therefore, suppress the leak-off of fluid, i.e. the drainage of fluid from the fracture into the surrounding material. This is achieved by choosing the permeability of the poroelastic material much smaller than the permeability associated with the fracture.
We apply displacement boundary conditions via (13) and (49) with ε 22 ¼ -1e-6 whilst all other strain components are ε ij ¼ 0. Note that the deformation is small enough that the fracture remains open throughout the time-dependent simulation. The fluctuation field at the surface is suppressed, u s ¼ 0. Moreover, we apply drained conditions (p ¼ 0) at the left and the right boundary and undrained boundary conditions (w ¼ 0) at the bottom and the top boundary of the sample. For illustration purposes, we introduce an apparent fracture pressure p F :¼ p ð1 À aÞ by incorporating the degradation function. The temporal evolution of the apparent fracture pressure in the computational domain is shown in Fig. 5 where one observes the successive decay of pore pressure. After having solved the transient problem, we evaluate the response of the structure via the overall outflux of pore fluid via the drained boundary defined as w ¼ w F ⋅n; and w ¼ w⋅n for the sharp and the diffuse interface model, respectively. Moreover, we use the definitions (19) and (56) to compute the overall stress σ 22 from the sub-scale fields. Finally, we execute all investigations using the interpolation functions based on the conventional damage field d as well as based on the modified damage field d * . FE discretization of the diffuse interface at constant smoothing length: We choose γ ¼ γ * ¼ 2e À 2 m and vary the size of the finite element in the vicinity of the interface. The element size is quantified by the element-specific scalar h e which is defined as the length of the longest edge of the used triangular elements. The outflux via the drained boundary and the homogenized stress response over time are plotted in Figs. 6 and 7.
We observe that all solutions tend towards a converged state if h e is chosen sufficiently small. Henceforth, we consider solutions computed with h e =γ ¼ 0:1 (conventional model) and h e =γ ¼ 1 (modified model) as converged. However, we also observe differences in the material response of the conventional and the modified model: The convergence rate seems to be higher for the modified model.

Size of the smoothing length:
We choose now the element size such that h e =γ ¼ 0:1 for the conventional and h e =γ ¼ 1 for the modified model and vary the smoothing length γ. For convenience, we set γ * ¼ γ in the model using the modified damage field d * . The overall outflux of fluid and the homogenized stress response are plotted in Figs. 8 and 9 and compared with the corresponding sharp interface solution, denoted as reference.
Surprisingly, we observe that the conventional formulation does not converge towards the correct solution (σ 22 →0 for t→∞). 2 The reason for this somewhat intriguing observation is that, e.g. under compression, the volumetric strain in the damage zone, which is directly related to the fluid pressure in the fracture, is concentrated in one single element layer along the fracture zone. In practice that means that the smoothing length for the fluid pressure in the fracture zone corresponds to the element size h e whilst the smoothing length γ of the interpolation function is chosen as a material parameter. In other words, the smoothing length of the   However, κ is chosen very small and is not responsible for the persisting stress e. g. observed in Fig. 7 a) and 9 a).  pressure and the solid stress are not matching. By contrast, the modified formulation results in a very good approximation of the reference solution if γ is chosen sufficiently small. This behavior can be explained by the fact that d * ¼ 1 ¼ const: over a band of 2 γ * . Hence, the problem is not mesh dependent as long as h e is chosen small enough. Moreover, the smoothing length of the pressure and the solid stress are matching. Nevertheless, it is important to note that even the modified formulation requires a very high resolution at the interface for a good agreement with the reference. This high resolution, however, might limit the applicability of the method for more complex problems due to numerical costs.

Validation against an analytical solution:
Finally, we compare the numerical results with the analytical solution of the so-called Terzaghi problem presented in (Verruijt, 2013). For that purpose, we modify the setting in Fig. 4 such that, instead of applying a strain ε 22 , we prescribe a constant initial pressure pðt ¼ 0Þ ¼ p 0 in the fracture with drained conditions, cf. Fig. 10 a). Hence, the setup is converted into an 1D problem of a thin poroelastic layer under homogeneous confining pressure. We can compute the strain rate of the poroelastic layer as with λ ¼ K À 2 3 G. Taking into account that, as a loading condition, the total stress is kept constant throughout the process, i.e. _ σ ¼ 0, we can insert the strain rate into the 1D continuity equation of the porous layer which yields the uncoupled diffusion equation In our case, the thin layer of a poroelastic medium is representing a fluid conduit, i.e. the solid phase is absent and the diffusion equation simplifies towards where the consolidation coefficient is defined as The analytical solution of the problem yields The pressure profiles, observed at different time steps, are plotted in Fig. 10 b) for the modified damage field d * (γ ¼ γ * ¼ 4e-3 m, h e =γ ¼ 1) in comparison with the analytical results. We find an excellent agreement between the modified diffuse interface model and the analytical solution.

RVE with two intersecting fractures
In the second study, we address the case of a simple fracture network consisting of two perpendicular intersecting fractures. As discussed in (Quintal et al., 2014;Rubino et al., 2013), this setup can be seen as a prototype to investigate seismic attenuation in fluid-saturated fracture networks. In practice, this setup is prone to seismic attenuation associated with viscous fluid transport from one fracture into the other under uniaxial loading where the fracture normal to the loading direction is compressed and fluid is squeezed out into the connected second fracture.
We adopt this setup as shown in Fig. 11 and apply an overall loading of the RVE by ε 22 ¼ À 1e À 6 in a sense of a stress relaxation experiment whilst the remaining strain components are zero. We use periodic boundary conditions for u and p as introduced in Sections 3 and 5. Again, the sharp interface formulation serves as reference solution. The apparent fracture fluid pressure distribution in the RVE is plotted in Fig. 12 for three example time steps. FE discretization of the diffuse interface at constant smoothing length: The results are shown in Fig. 13 in terms of the homogenized stress over time for the conventional and the modified diffuse interface model. We keep γ ¼ γ * ¼ 2e À 2 m constant and vary the ratio h e = γ as executed in the previous section. We find that (i) the homogenized stress shows the behavior of a stress relaxation curve with a transition from instantaneous response to the relaxed state t→∞ around t ¼ 30 s.
We find that (ii) the problem based on the modified formulation is less sensitive for the finite element discretization and that a ratio h e = γ ¼ f0:1; 1g gives a reasonably well converged solution for the conventional and modified model, respectively.
Size of the smoothing length: In a second step, we fix the ratio h e =τ and vary the smoothing length γ as shown in Fig. 14. We observe that the conventional model does not converge towards the correct result. By contrast, we find for the modified model that the smaller γ ¼ γ * the better the approximation of the reference solution. Whilst γ ¼ γ * ¼ 4e À 3 m matches the equilibrated state t→∞, the instantaneous plateau is slightly underestimated, and the transition zone is shifted to faster times by a factor � 0:5: Interpolation of the permeability: The transition zone in the stress relaxation curve between instantaneous response for t→0 and equilibrium for t→∞ is mostly dominated by how fast the fluid is being transported through the fracture zone. Hence, this property depends on the permeability that we assume for the fracture zone. We, therefore, use the scaling parameter ν as introduced in (45) to interpolate the permeability between the fracture zone and the poroelastic domain. The results are shown in Fig. 15, where only the modified damage model was used.
As shown in Fig. 3, the higher permeability of the fracture becomes more dominant in the transition zone between fracture and poroelastic matrix for the arithmetic mean ν ¼ 1 whilst the transition zone is mostly dominated by the lower poroelastic permeability for the harmonic mean ν ¼ À 1. The consequence can be found in Fig. 15: The arithmetic mean, associated with a higher permeability in the transition zone between fracture and poroelastic background, leads to a faster stress relaxation whilst the harmonic mean ν ¼ À 1 leads to an overestimation of the relaxation time by orders of magnitude. In our case, the choice ν ¼ 0 seems to give the best approximation of the reference results. In analogous fashion, this observation can also be made for the evaluation in frequency domain shown in Fig. 16. These curves are computed via a Fast Fourier Transformation of the stress relaxation curves plotted in  is used for a quantification of the attenuation, where Q : ¼ À R= IðC 2222 Þ.

Discussion
In this paper, we applied a fracture phase field method to model diffusive fluid transport in simple stationary fracture networks. Hence, fracture propagation was ignored. This scenario can be seen as a subproblem of hydraulic fracturing. The research issue was to explore the predictive performance of a fracture phase field method in such a simplified scenario, i.e. how well can the method handle the hydromechanical coupling between fractures and a poroelastic background material. There are two major arguments to investigate this particular sub-problem. (i) Diffusive fluid-transport in fracture networks associated with attenuation of seismic waves is well understood and a relevant mechanism in seismic exploration in the context of hydraulic stimulation. (ii) There are sharp interface formulations at hand that can serve as reference models at benchmarking the diffuse interface method.
We found that the method based on the conventional interpolation functions established in the literature does not converge towards the reference solution if the smoothing length γ→0. Moreover, the method might even lead to nonphysical results. The reason is that, e.g. under compression, the volumetric strain in the damage zone, which is directly related to the fluid pressure in the fracture, was found to be concentrated in one single element. In practice that means that the smoothing length for the fluid pressure in the fracture zone corresponds to the element size h e whilst the smoothing length γ of the interpolation function is chosen as a material parameter.
To overcome this inherent mesh dependence of the hydromechanical coupling, we modified the interpolation by introducing a modified damage variable d * . The latter is constructed in a way that d * ¼ 1 ¼ constant in a band of 2 γ * along the interface. We showed that, doing so, the diffuse interface model converges towards the reference solution. However, the smoothing length γ * had to be chosen pretty small, yet still orders of magnitude larger than the true interface thickness. Nevertheless, the high resolution of the transition zone needed for quantitatively reasonable solutions obviously leads to high numerical costs if applied to more complex setups.
Altogether, it is important to point out that the ability to predict and monitor hydraulic stimulation by numerical simulations is of utmost importance to assess the risks of the process, e.g. seismic events, and the societal acceptance of this technology. Fracture phase fields methods for fracturing of poroelastic media are well-established in literature, and it is accepted that such methods are able to predict hydraulic stimulation qualitatively. However, having our critical observations for an even simpler sub-problem in mind, there is a need (i) to further investigate fracture phase field methods in fluid-saturated porous media, and (ii) to develop quantitative benchmarks in realistic geophysical scenarios.