Modeling Fractures in a Poro-Elastic Medium

— Un mode`le de fracture dans un milieu poro-e´lastique — Nous pre´sentons un mode`le de fracture dans un milieu poro-e´lastique. Il de´crit la fracture comme une courbe ou une surface, selon la dimension, l’e´paisseur e´tant incorpore´e dans l’e´quation de l’e´coulement dans la fracture. La discre´tisation utilise des e´le´ments ﬁnis mixtes pour le ﬂuide et des e´le´ments ﬁnis continus pour le de´placement du milieu. Le sche´ma est re´solu par un algorithme qui de´couple, d’une part, le calcul du de´placement de celui de l’e´coulement, et d’autre part, le calcul de l’e´coulement dans le re´servoir de celui dans la fracture. Le mode`le est illustre´ par un essai nume´rique. Abstract — Modeling Fractures in a Poro-Elastic Medium — We present a fracture model in a poro-elastic medium. The model describes the fracture as a curve or surface according to the dimension, the width of the crack being included into the equation of ﬂow in the fracture. The discretization uses mixed ﬁnite elements for the ﬂuid and continuous ﬁnite elements for the porous medium’s displacement. The numerical scheme is solved by an algorithm that decouples, on one hand, the computation of the mechanics from that of the ﬂuid, and on the other hand, the computation of the ﬂow in the reservoir from that in the fracture. The model is illustrated by a numerical experiment.


INTRODUCTION
In petroleum and environmental engineering, studies of multiscale and multiphysics phenomena such as reservoir deformation, surface subsidence, well stability, sand production, waste deposition, hydraulic fracturing, and CO 2 sequestration [1] require a clear understanding of both the fluid flow and solid phase mechanical response.
Traditionally, fluid flow problems were focused on reservoir flow whereas the influence of porous media deformation on pore pressure was usually simplified and approximated as a constant rock compressibility that could not account for rock response especially in naturally fractured and/or stress-sensitive reservoirs [2]. In this setting, models for narrow fractures can be derived by treating the width as a small parameter e and letting this parameter tend to zero. One approach in this direction can be found in the theoretical analysis of Morales and Showalter [3,4], where only flow is considered, the crack has a flat basis and vertical height of the order of e, and the pressure is continuous at the interfaces. The former reference deals with the divergence form of Darcy's law and the latter with the mixed form. Another approach consists in treating the fracture as a thin domain in the framework of domain decomposition. We refer the reader to the extensive work of Alboin et al. To overcome the limitations of decoupled flow simulation, solving the coupled fluid flow and mechanics model in porous media setting has become more feasible with emerging computational machine power.
There are three approaches that are currently employed in the numerical coupling of fluid flow and the mechanical response of the reservoir solid structure [9,12]. They have been referred to in the literature as fully implicit, loose or explicit coupling, and iterative coupling. The fully implicit involves solving all of the governing equations simultaneously and involves careful implementation with substantial local memory and complex solvers. The loosely or explicitly coupled is less accurate and requires estimates of when to update the mechanic response. Iterative coupling is a sequential procedure where either the flow or the mechanics is solved first, followed by incorporating the latest solution information. At each step, the procedure is iterated until the solution converges within an acceptable tolerance. This last approach appears to be more scalable on emerging exascale computer platforms; it can also be used as a preconditioner for the fully implicit coupling. But of course, iterative coupling schemes must be carefully designed. For instance, two iterative coupling schemes, the undrained split and the fixed stress split, are of particular interest; they converge for slightly compressible linearized flow [13]. In contrast, a Von Neumann stability analysis shows that the drained split and the fixed strain split iterative methods are not stable for the same model problem [9].
Hydraulic fracking is one of the primary ways manufacturers retrieve natural gas. Hence, modeling fractures in porous media, that are either naturally present or are created by stimulation processes, is a high profile topic, but it is also a source of additional and challenging complexities. This work focuses on the simulation of the time-dependent flow of a fluid in a deformable porous medium that contains a crack. The medium in which the crack is embedded is governed by the standard equations of linear poro-elasticity and the flow of the fluid within the crack is governed by a specific channel flow relation. The two key assumptions of this channel flow equation are: -the width of the crack is small compared to its other relevant dimensions, and is such that the crack can be represented as a single surface. The relevant kinematical data is the jump in the displacement of the medium across the crack, i.e. the crack's width; -the permeability in the crack is much larger than in the reservoir. One advantage of treating the crack as a single surface is that it alleviates the need for meshing a thin region, thus avoiding all the complexities related to anisotropic elements.
In this paper, we formulate a discretization that is suitable for irregular and rough grids and discontinuous full tensor permeabilities that are often encountered in modeling subsurface flows. To this end, we develop a multiphysics algorithm that couples Multipoint Flux Mixed Finite Element (MFMFE) methods for the fluid, both in the medium and in the crack, with continuous Galerkin finite element methods (CG) for the elastic displacement. The MFMFE method was developed for Darcy flow in [14][15][16]. It is locally conservative with continuous fluxes and can be viewed within a variational framework as a mixed finite element method with special approximating spaces and quadrature rules. It allows for an accurate and efficient treatment of irregular geometries and heterogeneities such as faults, layers, and pinchouts that require highly distorted grids and discontinuous coefficients. The resulting discretizations are cell-centered with convergent pressures and velocities on general hexahedral and simplicial grids. The MFMFE method is extended to poro-elasticity in [17]; the analysis is based on the technique developed by Phillips and Wheeler in [18].
We solve the numerical scheme by an extension to discrete fractures of the fixed stress splitting algorithm. Its salient features at each time step are: -we decouple the computation of the displacement from that of the fluid flow until convergence; -we decouple the computation of the fluid flow in the reservoir from that of the flow in the crack until convergence.
The model described here does not include crack propagation, i.e. the crack front is stationary. However, the work presented here can be viewed as a starting point toward hydraulic fracture modeling in which non planar crack propagation is simulated. Our future plans include coupling the present software with the HYFRAC boundary element method developed by Mear and discussed in [19], as well as treating multiple cracks or multiphase flow in the reservoir and non Newtonian flow in the fracture.
This article is organized as follows. The modeling equations are defined in Section 1. Section 2 summarizes theoretical results on the linear model. The numerical scheme and algorithm are described in Section 3. A numerical experiment illustrating the model is presented in Section 4. We end with some conclusions.

PROBLEM FORMULATION
Let the reservoir X be a bounded domain of R d d ¼ 2 or 3, with piecewise smooth Lipschitz boundary @X and exterior normal n. Let the fracture C b X be a simple piecewise smooth curve with endpoints a and b when d ¼ 2 or a simple piecewise smooth surface with piecewise smooth Lipschitz boundary @C when d ¼ 3. To simplify, we denote partial derivatives with respect to time by the index t.

Equations in X\C
The displacement of the solid is modeled in XnC by the quasi-static Biot equations for a linear elastic, homogeneous, isotropic, porous solid saturated with a slightly compressible viscous fluid. The constitutive equation for the Cauchy stress tensor r por is: r por ðu; pÞ ¼ rðuÞ À a p I ð1:1Þ where I is the identity tensor, u is the solid's displacement, p is the fluid pressure, r is the effective linear elastic stress tensor: Here, k > 0 and G > 0 are the Lame´constants, a > 0 is the dimensionless Biot coefficient, and eðuÞ ¼ ðru þ ru T Þ=2 is the strain tensor. Then the balance of linear momentum in the solid reads: where f is a body force. For the fluid, we use a linearized slightly compressible single-phase model. Let p r be a reference pressure, q f > 0 the fluid phase density, q f ;r > 0 a constant reference density relative to p r , and c f the fluid compressibility. We consider the simplified case when q f is a linear function of pressure: Next, let u Ã denote the fluid content of the medium; it is related to the displacement and pressure by: where u 0 is the initial porosity, and M a Biot constant.
The velocity of the fluid v D in XnC obeys Darcy's Law: where K is the absolute permeability tensor, assumed to be symmetric, bounded, uniformly positive definite in space and constant in time, l f > 0 is the constant fluid viscosity, g is the gravitation constant, and g is the distance in the vertical direction, variable in space, but constant in time. The fluid mass balance in XnC reads: where q is a mass source or sink term taking into account injection into or out of the reservoir. The equation obtained by substituting (1.4) and (1.5) into (1.7) is linearized through the following considerations. The fluid compressibility c f is of the order of 10 À5 or 10 À6 , i.e. it is small. The fraction p=M is small and the divergence of u is also small. Therefore these terms can be neglected. With these approximations, when dividing (1.7) by q f ;r , substituting (1.6), and settingq ¼ q=q f ;r , we obtain: Kðr p À q f ;r gr gÞ Thus the poro-elastic system we are considering for modeling the displacement u and pressure p in XnC is governed by (1.1), (1.3) and (1.8).

Equation in C
The trace of p on C is denoted by p c and r is the surface gradient operator on C, it is the tangential trace of the gradient. These quantities are well defined for our problem. The width of the fracture is represented by a nonnegative function w defined on C; it is the jump of the displacement u in the normal direction. Since the medium is elastic and the energy is finite, w must be bounded and must vanish on the boundary of the fracture. We adopt a channel flow relation for the crack in which the volumetric flow rate Q on C satisfies [1]: and the conservation of mass in the fracture satisfies: Here, q I is a known injection term into the fracture, and q L is an unknown leakoff term from the fracture into the reservoir that guarantees the conservation of mass in the system. Approximating q f by q f ;r in the time derivative, linearizing the diffusion term as in the previous section, dividing by q f ;r , and settingq In order to specify the relation between the displacement u of the medium and the width w of the fracture, let us distinguish the two sides (or faces) of C by the superscripts þ and À; a specific choice must be selected but is arbitrary. To simplify the discussion, we use a superscript H to denote either þ or À. Let X H denote the part of X adjacent to C H and let n H denote the unit normal vector to C exterior to X H , H ¼ þ; À. As the fracture is represented geometrically by a figure with no width, then n À ¼ Àn þ . For any function g defined in XnC that has a trace, let g H denote the trace of g on C H , H ¼ þ; À. Then we define the jump of g on C in the direction of n þ by: The width w is the jump of u Á n À on C: Therefore the only unknown in (1.9) is the leakoff term q L .
Summarizing, the equations in XnC are (1.3) and (1.8), and the equation in C is (1.9); the corresponding unknowns are u, p andq L . These equations are complemented in the next section by interface, boundary and initial conditions.

Interface, Boundary, and Initial Conditions
The balance of the normal traction vector and the conservation of mass yield the interface conditions on each side (or face) of C: Then the continuity of p c through C yields: The conservation of mass at the interface gives: General conditions on the exterior boundary oX of X can be prescribed for the poro-elastic system, but to simplify, we assume that the displacement u vanishes as well as the flux Kðr p À q f ;r gr gÞ Á n. According to the above hypotheses on the energy and medium, we assume that w is bounded in C and vanishes on oC. Finally, the only initial data that we need are the initial pressure and initial porosity, p 0 and u 0 . Therefore the complete problem statement, called Problem (Q), is: Find u, p, andq L satisfying (1.1), (1.3), (1.8) in XnC and (1.9) in C, for all time t 20; T ½, with the interface conditions (1.10) and (1.11) on C: where p c ¼ p jC and: with the boundary conditions: and the initial condition at time t ¼ 0: Problem (Q) is highly non linear and its analysis is outside the scope of this work. Therefore we present here some theoretical results on a linearized problem where (1.12) is substituted into the first term of (1.9) while the factor w 3 in the second term is assumed to be known. Knowing w 3 amounts to linearizing the nonlinear term with respect to w, such as could be encountered in a time-stepping algorithm.

Variational Formulation
It is convenient (but not fundamental) to generalize the notation of Section 1.2 by introducing an auxiliary partition of X into two non-overlapping subdomains X þ and X À with Lipschitz interface C containing C, X H being adjacent to C H , H ¼ þ; À. The precise shape of C is not important as long as X þ and X À are both Lipschitz. Let C H ¼ oX H nC; for any funtion g defined in X, normed by the graph norm: The space for the displacement is: with the norm of W d : The space for the pressure is more subtle: it is H 1 in X, but in C, it is: equipped with the norm: where H 1=2 ðCÞ is the space of traces on C of H 1 ðXÞ functions. Thus p belongs to the space: equipped with the graph norm: Finally, the space for the leakoff variableq L is H 1 w ðCÞ 0 , the dual space of H 1 w ðCÞ, equipped with the dual norm. Then, by multiplying (1.3), (1.8) and (1.9) with adequate test functions, applying Green's formula, using the interface conditions (1.10) and (1.11), the boundary conditions (1.13) and (1.14), substituting (1.12) into the first term of (1.9), and assuming sufficiently smooth data, we obtain the variational formulation: Find p 2 H 1 ð0; T ; L 2 ðXÞÞ \ L 1 ð0; T ; QÞ, u 2 H 1 ð0; T ; V Þ, and q L 2 L 2 ð0; T ; H 1 w ðCÞ 0 Þ solution a.e. in 0; T ½ of: 2GðeðuÞ; eðvÞÞ þ kðr Á u; r Á vÞ À aðp; r Á vÞ þðp c ; ½v C Á n þ Þ C ¼ ðf ; vÞ; 8v 2 V ð2:7Þ where the scalar products are taken in XnC, with the initial condition at time t ¼ 0

Existence and Uniqueness
The analysis of Problem (2.7-2.10) relies on suitable properties of H 1 w ðCÞ, that in turn depend on the regularity of w. For the practical applications we have in mind, w has the following properties when d ¼ 3; the statement easily extends to d ¼ 2. Hypothesis 1. The non negative function w is H 1 in time and is smooth in space away from the crack front, i.e. the boundary @C. It vanishes on @C and in a neighborhood of any point of @C, w is asymptotically of the form: where y is locally parallel to the crack front, x is the distance to the crack, and f is smooth.
With this hypothesis, we can validate the integrations by parts required to pass from Problem (Q) to its variational formulation (2.7-2.9) and we recover the leakoff unknownq L from (2.11), thus showing equivalence between all three formulations. This equivalence allows to ''construct" a solution of Problem (Q) by discretizing (2.14) with a Galerkin method, more precisely by approximating p in a truncated basis, say of dimension m, of the space Q. This leads to a square system of linear Ordinary Differential Equations of order one and of dimension m, with an initial condition. A close examination of this system's matrices shows that it has a unique solution. Then reasonable regularity assumptions on the data, Hypothesis 1, and the theory developed by [18], yield a set of basic a priori estimates on the discrete pressure, say p m , and displacement uðp m Þ. These are sufficient to pass to the limit in the discrete version of (2.14), but, as is the case of a poro-elastic system, they are not sufficient to recover the initial condition on p. This initial condition can be recovered by means of a suitable a priori estimate on the time derivative p 0 m , but it requires a slightly stronger regularity assumption on the data (still fairly reasonable) and an assumption on the time growth of w. It allows the time evolution of a crack, but it cannot be used in the formation of a crack, since it forbids an opening of a region that is originally closed. More precisely, we make the following hypothesis. Hypothesis 2. There exists a constant C such that: a:e: in CÂ0; T ½; w 0 C w ð2:16Þ With this additional assumption, we can prove existence and uniqueness of a solution of Problem (Q). Theorem 1. If f is in H 2 ð0; T ; L 2 ðXÞ d Þ,q in L 2 ðXÂ0; T ½Þ, q I in H 1 ð0; T ; L 2 ðCÞÞ, p 0 in Q, and w satisfies Hypotheseses 1 and 2, then the linearized problem (2.7-2.10) has one and only one solution p in H 1 ð0; T ; L 2 ðXÞÞ \ L 2 ð0; T ; QÞ, u in H 1 ð0; T ; V Þ, andq L in L 2 ð0; T ; H 1 w ðCÞ 0 Þ. This solution depends continuously on the data.

Heuristic Mixed Formulation
A mixed formulation for the flow is useful in view of discretization because it leads to locally conservative schemes. The mixed formulation we present in this short section is derived heuristically in the sense that it is written in spaces of sufficiently smooth functions in the fracture. As the discrete functions are always locally smooth, all terms in the discrete formulation will make sense. For the pressure in XnC, we define the auxiliary velocity z by: and for the pressure in C, we define a surface velocity f by: f ¼ Àr p c À q f ;r gg À Á ð2:18Þ Let: In view of (1.13), (1.11), and the regularity in Theorem 1, we have z 2 Z and we see that it must satisfy the essential jump condition: We use the same space V for u, but we reduce the regularity of p and set: As the trace of functions in Q are no longer meaningful, we introduce an additional variable p c in the space: and we assume that the leakoff term is in H C . Then the mixed variational formulation reads (to simplify, we do not specify the dependence of the spaces on time): Find u 2 V , p 2 Q, p c 2 H C ,q L 2 H C , z 2 Z, and f 2 Z C such that (2.7), (2.10), and (2.19) hold: 2GðeðuÞ; eðvÞÞ þ kðr Á u; r Á vÞ À aðp; r Á vÞ þðp c ; ½v C Á n þ Þ C ¼ ðf ; vÞ; 8v 2 V and for all h in Q: :21Þ 8q 2 Z; K À1 z; q À Á ¼ ðp; r Á qÞ À ðp c ; ½q C Á n þ Þ C þ r q f ;r gg À Á ; q À Á ð2:22Þ 8l 2 Z C ; ðw 3 f; lÞ C ¼ ðp c ; r Á ðw 3 lÞÞ C þ ðrðq f ;r ggÞ; w 3 lÞ C ð2:23Þ It follows from (2.22) that p belongs to H 1 ðXÞ and p c ¼ pj C . Hence p c acts like a Lagrange multiplier that enforces the continuity of the traces on C. Moreover, when the solutions are sufficiently smooth, the mixed and primitive formulations are equivalent. In addition, q L can be immediately eliminated by substituting (2.19) into (2.21); this yields: . In order to avoid handling curved elements or approximating curved surfaces, we assume that both oX and the crack C are polygonal or polyhedral surfaces. We shall refer to this scheme as the MFMFE-CG scheme. A convergence analysis of the MFE-CG scheme for the poro-elasticity system without a crack is done in [17], where two versions of the method: one with a symmetric quadrature rule on simplicial grids or on smooth quadrilateral and hexahedral grids, and one with a non-symmetric quadrature rule on rough quadrilateral and hexahedral grids are treated. Theoretical and numerical results demonstrated first order convergence in time and space for the fluid pressure and velocity as well as for solid displacement. In the formulation below, we restrict our attention to the symmetric quadrature rule for the flow.

MFMFE Spaces
We first describe the interior discretization. Let X be exactly partitioned into a conforming union of finite elements of characteristic size h, i.e. without hanging nodes. For convenience we assume that the elements are quadrilaterals in 2D and hexahedra in 3D. Let us denote the partition by T h and assume that it is shape-regular [20]. For convenience, we also assume that T h triangulates exactly X þ and X À ; in particular, no element crosses C. The displacement, velocity and pressure finite element spaces on any physical element E are defined, respectively, via the vector transformation: and via the scalar transformation: where F E denotes a mapping from the reference element E to the physical element E, DF E is the Jacobian of F E , and J E is its determinant.
whereV ðÊÞ,ẐðÊÞ andQðÊÞ are finite element spaces on the reference elementÊ. Note that since C is an open set and V h is contained in V , the functions of V h are continuous on oC, i.e. they have no jump on oC.
Let P k denote the space of polynomials of degree k in two or three variables, and let Q 1 denote the bilinear or trilinear polynomial spaces in two or three variables, according to the dimension.

Convex Quadrilaterals
In the case of convex quadrilaterals,Ê is the unit square with verticesr 1 ¼ ð0; 0Þ T ,r 2 ¼ ð1; 0Þ T ,r 3 ¼ ð1; 1Þ T , and r 4 ¼ ð0; 1Þ T . Denote by r i , i ¼ 1; . . . ; 4 the corresponding vertices of E. In this case, F E is the bilinear mapping given as: þ r 4 1 Àx ð Þŷ the space for the displacement is: V ðÊÞ ¼ Q 1 ðÊÞ 2 and the space for the flow is the lowest order BDM 1 [23] space:Ẑ where r and s are real constants.

A Quadrature Rule
The integration on a physical element is performed by mapping to the reference element and choosing a quadrature rule onÊ. Using the Piola transformation, we write ðK À1 Á; ÁÞ in (2.22) as: where: Then, we denote the trapezoidal rule onÊ by TrapðÁ; ÁÞÊ: where fr i g k i¼1 are the vertices ofÊ. The quadrature rule on an element E is defined as: and we define ðK À1 z; qÞ Q by summing (3.8) on all elements E of T h . As the trapezoidal rule induces a scalar product on the above finite element spaces that yields a norm uniformly equivalent to the L 2 norm [14,15], and the tensor K is symmetric, bounded, and uniformly positive definite, this quadrature rule also yields a norm on these spaces uniformly equivalent to the L 2 norm.

Discretization in the Fracture
Since C is assumed to be polygonal or polyhedral, we can map each line segment or plane face of C onto a segment in the x 1 line (when d ¼ 2) or a polygon in the x 1 À x 2 plane (when d ¼ 3) by a rigid-body motion that preserves both surface gradient and divergence, maps the normal n þ into a unit vector along x 3 , for exemple Àe 3 , and whose Jacobian is one. After this change in variable, all operations on this line segment or plane face can be treated as the same operations on the x 1 axis or x 1 À x 2 plane. To simplify, we do not use a particular notation for this change in variable, and work as if the line segments or plane faces of C lie on the x 1 line or x 1 À x 2 plane. Let S i , 1 i I, denote the line segments or plane faces of C; to simplify, we drop the index i. Let T S;h be a shape regular partition of S, for instance the trace of T h on S (but it can be something else) let e denote a generic element of T S;h , with reference element e, and let the scalar and Piola transforms be defined by the same formula as above, but with respect to e instead of E. Then we define the finite element spaces on C by: with:

The Scheme
Let N ! 1 be a fixed integer, Át ¼ T =N the time step, and t i ¼ iÁt, 0 i N , the discrete time points. We define the approximations f n of f ðt n Þ andq n ofqðt n Þ for almost every x 2 X þ [ X À or X by: x; t ð Þd t

ð3:12Þ
Similarly, we define the approximationq n I ofq I ðt n Þ for almost every s 2 C by: On the other hand, we use point values for w in time, i.e. define: 8s 2 C; w n s ð Þ ¼ w s; t n ð Þ Then we propose an iterative algorithm for solving the fully discrete version of (2.7), (2.10), (2.20), (2.22), (2.23), (2.24) that uses the quadrature rule of Sections 3.2 and 3.3. A flowchart of the scheme is given in Figure 1. It uses the Fixed Stress split iterative algorithm of Mikelic´and Wheeler [25]. The mean stress is given by where r por is defined by (1.1). Thus the discrete mean stress reads: For a given displacement and mean stress, the algorithm decouples sequentially the computation of the flow in the reservoir from that in the fracture.
Recall that the given initial pressure p 0 is sufficiently smooth and has a trace, say p 0 c , on C; this is the initial pressure in C.
P h are suitable approximation operators from Q into Q h and H C into H C;h respectively. Then the initial displacement is approximated by discretizing the elasticity equation (2.13) in XnC: This gives the initial mean stress: . Computing the solution relies on the work of Liu in [26]. The system (3. 18-3.19) is solved in Section 4 below by means of a mimetic formulation [27], which is closely related to a mixed finite element method when using quadrilateral elements.

NUMERICAL EXPERIMENT
We provide the following numerical example as a simple illustration of the model and its predictive capability. Following the flowchart in Figure 1, we have implemented the fixed stress iterative coupling scheme with the inner iteration between reservoir flow and fracture flow into the IPARS reservoir simulation code.
The domain is taken to be the cube X ¼ ð0; 10Þ 3 m, with a square fracture C centered on the plane y ¼ 5 f gm of size 7:5 Â 7:5 m. The domain is discretized into 5 Â 6 Â 5 structured hexahedral elements. There are uniform mesh widths of 2 m in the x and z directions, and mesh widths of 2:475; 2:475; 0:05; 0:05; 2:475; f 2:475g m in the y direction.   The initial fluid pressure in both the reservoir and fracture is taken to be p 0 ¼ 3:5 Â 10 6 Pa. The external fluid boundary conditions are set to be p ¼ p 0 on @X, which will allow fluid to escape as pressure rises. The external mechanical boundary conditions are such that the bottom face is completely pinned: u ¼ 0 m; the four lateral faces are traction free: r por n ¼ 0 psi; and the top face has an overburden traction of r por n ¼ À3:8; 0; 0 ð ÞÂ10 6 psi. The stabilization parameter is taken to be b ¼ 1:0 Â 10 À4 . The remaining input parameters are summarized in Table 1. The Lame´constants are given by k ¼ Em=½ð1 þ mÞð1 À 2mÞ and G ¼ E=½2ð1 þ mÞ. Fluid is injected into the center of the fracture with a constant rate of q I ¼ 1 000 kg Á s À1 , which will induce a pressure gradient in the fracture and leak off into the domain. Both the leak off rate q L and the fracture pressure p c are unknowns determined by the coupling of the reservoir and fracture flow models. The fluid pressure computed in the fracture plane is shown in Figure 2. The pressurized fracture also induces a traction on C, which creates a discontinuity in the displacement field, allowing the calculation of a dynamic width, which is shown in Figure 3.

CONCLUSIONS
We have studied the numerical approximation of a fracture model in a poroelastic medium, where the fracture is represented as a curve or surface and its width is incorporated into the flow equation in the fracture. We have used a MFMFE method for the flow in the reservoir and a Mimetic Finite Difference (MFD) method in the fracture, and a Continuous Galerkin (CG) method for the geomechanics. The scheme was solved by a fixed stress split iterative coupling for the flow and mechanics and an iterative algorithm coupling the flow in the reservoir and in the crack. A numerical experiment illustrated the relevance of the mechanical model and the efficiency of the numerical method and algorithm.