Andro Mikelic
email: andro.mikelic@univ-lyon1.fr

Mary F Wheeler

Thomas Wick

Andro Mikelić

A quasistatic phase field approach to pressurized fractures *

Keywords: Pressurized cracks, phase-field, quasistatic model, incremental approximation, Biot's consolidation equations, MSC classcode 35Q86, 47J20, 49J40, 74R10

published or not. The documents may come

A quasistatic phase field approach to pressurized fractures

Introduction

In petroleum and environmental engineering, mathematical and numerical studies of multiscale and multiphysics phenomena such as reservoir deformation, surface subsidence, well stability, sand production, waste deposition, hydraulic fracturing, and CO 2 sequestration are receiving increasing attention. Hydraulic fracturing is fracturing of various rock layers by a pressurized liquid and is a technique used to release petroleum, natural gas (shale gas, tight gas, and coal seam gas) for extraction. The first use of hydraulic fracturing was in 1947 and in 1949 it was commercialized. As of 2010, 60% of all new oil and gas worldwide production was from fractured wells and as of 2012, 2.5M fracking jobs were performed.

Modeling is crucial for the understanding and the prediction of the physical behavior of fractured systems [START_REF] Dean | Hydraulic-Fracture Predictions with a fully coupled geomechanical reservoir simulator[END_REF]. Issues include evaluation of injection enhancement for various "frac-job" scenarios and modeling the interaction between hydraulic and discrete fractures. Specific environmental concerns involve possible contamination of groundwater, risks to air quality, and migration of gases and hydraulic fracturing chemicals to the surface and surface contamination from spills and flowback and health effects of these.

The computational modeling of the formation and growth of fluid filled fractures in a poroelastic media is difficult with complex crack topologies. Cracks can grow, form and interact. Tracking sharp cracks using classical methods can be computationally intensive, especially in realistic heterogeneous formations. These difficulties have recently been overcome by diffusive crack modeling based on the introduction of a crack phase field. We refer to the work of Bourdin, Francfort, and Marigo [START_REF] Bourdin | The Variational Approach to Fracture[END_REF] in which a regularized crack surface functional is shown to converge to a sharp crack surface functional. Our goal here is to generalize this approach to growing pressurized cracks in poroelastic media. Reports with our analysis and calculations with a fully implicit incremental formulation can be found in [START_REF] Mikelić | A phase field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF][START_REF] Mikelić | Phase-field modeling of pressurized fractures in a poroelastic medium[END_REF]. Our computational approach uses iterative coupling, in which mechanics and flow equations are solved sequentially. In this paper we analyze a quasistatic formulation of the phase field model that involves the introduction of a novel incremental approximation that differs from [START_REF] Mikelić | A phase field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF][START_REF] Mikelić | Phase-field modeling of pressurized fractures in a poroelastic medium[END_REF] in that the linear elasticity and nonlinear phase field problem are decoupled, a nonlinear term is treated using a discrete derivative. In [START_REF] Mikelić | A phase field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF][START_REF] Mikelić | Phase-field modeling of pressurized fractures in a poroelastic medium[END_REF] our consideration was limited to an incremental formulation. Specifically in [START_REF] Mikelić | A phase field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF], the Euler-Lagrange equations were considered; whereas [START_REF] Mikelić | Phase-field modeling of pressurized fractures in a poroelastic medium[END_REF] treats the problem as energy minimization.

The significance of our new quasistatic model is as follows. To the best of our knowledge, this is the first phase field model to treat pressurized crack propagation in a poroelastic medium. The characteristics of the quasi-static model include:

• Crack growth is strongly imposed as an entropy condition;

• The system is closed with the Rice condition interpreted as an energy inequality.

These results are established rigorously as the limit of a corresponding incremental model, as the time step goes to zero.

The outline of our paper is as follows: In Section 2, we formulate a two-field energy functional coupled with a crack phase field, involving a time derivative nonnegativity constraint. Our formulation follows Francfort and Marigo's variational approach from [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and [START_REF] Bourdin | The Variational Approach to Fracture[END_REF] and extends the latter to pressurized cracks in a poroelastic medium. Existence of a solution for this quasistatic problem is proven in Section 3 by introducing an incremental formulation. This is achieved by subsequently first treating a phase field step followed by an elasticity step and then establishing a quasistatic limit. A finite element approximation to our coupled nonlinear system is described in Section 4. The solution of several benchmark and prototype problems are presented in Section 5. These include comparisons with an approach based on the anisotropic energy storage function from [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF]. In the Appendix we prove a technical regularity proposition.

Model formulation

A pressurized crack is contained in (0, L) 3 and propagates into a poroelastic medium Ω ⊂ (0, L) 3 , as shown on Figure 1. We derive a phase-field model for crack propagation. In contrast to crack propagation in an elastic medium, the quasistatic Biot equations can not be formulated as an energy minimization problem. Therefore extending the variational phase-field approach of Francfort, Marigo and others is not straightforward.

As in [START_REF] Mikelić | A phase field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF], we approach the problem by applying the fixed stress splitting algorithm [START_REF] Settari | Advances in coupled geomechanical and reservoir modeling with applications to reservoir compaction[END_REF][START_REF] Mikelić | Convergence of iterative coupling for coupled flow and geomechanics[END_REF]. Specifically, p = p B denotes the effective fluid pressure in the poroelastic medium Ω and p = p f is the fracture fluid pressure. Next, the quasistatic Biot system in Ω

∂ t 1 M p B + div (αu) + div { K B f η (ρ f g -∇p B )} = 0; (1)
-div {Ge(u) -α(p B -p 0 )I + σ 0 } = 0, (2)
is solved in two steps. The first step consists of solving the pressure equation ( 1), with a given displacement u and an enhanced Biot's modulus M , 1/ M = 1/M + α 2 /K dr , where K dr is the drained bulk modulus and α is Biot's coefficient. In the second step, Navier's system (2) for the displacements is solved for a given pressure p = p B . Convergence of the fixed stress split iterations was established in [START_REF] Mikelić | Convergence of iterative coupling for coupled flow and geomechanics[END_REF].

We remark that in this paper we do not solve for the fixed-stress pressure equation ( 1) but assume that the pressures p f and p B are given a priori. Our goal here is to focus on how these pressure variables are modeled in the fixed-stress elasticity equation (2) (present Section 2), well-posedness (Section 3), numerical approximations and examples (Sections 4 and 5). The non-trivial extension, formulating both equations of the fixed stress algorithm within a phase-field framework is the purpose of another study [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF].

In order to avoid a possible confusion due to a large number of parameters in the Biot equations, we present them in Table 1. Let the boundary of (0, L) 3 be denoted by ∂(0, L) 3 . We assume homogeneous Dirichlet conditions for the displacements on ∂ D (0, L) 3 of ∂(0, L) 3 . On the remaining part ∂ N (0, L) 3 , Neumann conditions are defined. Furthermore, the crack domain is C, supposed to be smooth, and

∂ N Ω = ∂ N (0, L) 3 ∪ ∂C.
The unit exterior normal to Ω is n.

We follow Griffith's criterion [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF] and suppose that the crack propagation occurs in the domain Ω when the elastic energy restitution rate reaches its critical value G c . If τ is the traction force applied at part of the boundary ∂ N Ω, then we associate to the crack C the following total energy

E(u, C) = Ω 1 2
Ge(u) : e(u) dx -

∂ N Ω τ • u dS - Ω αp B div u dx + G c H 2 (C), (3)
where H 2 (C) is the Hausdorff measure of the crack. Note that τ = -p f n on ∂C. This energy functional is minimized with respect to the kinematically admissible displacements u and any crack set satisfying a crack growth condition. The computational modeling of this minimization problem requires approximation of the crack location and of its length. This is achieved by regularizing the sharp crack surface topology in the solid using diffusive crack zones described by a scalar auxiliary variable. This variable is a phase-field that interpolates between the unbroken and the broken states of the material.

A thermodynamically consistent framework for phase-field models has been proposed by Miehe et al. in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations[END_REF]. They developed models for quasistatic crack propagation in elastic solids, together with incremental variational principles.

We introduce the time-dependent crack phase field ϕ, defined on (0, L) 3 × (0, T ). The regularized crack functional reads

Γ ε (ϕ) = (0,L) 3 ( 1 2ε (1 -ϕ) 2 + ε 2 |∇ϕ| 2 ) dx. (4)
This regularization of H 2 (C), in the sense of the Γ-limit when ε → 0, was used in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF].

Our further considerations are based on the fact that the evolution of cracks is fully dissipative in nature. First, the crack phase field ϕ is intuitively a regularization of 1 -1 C and we impose its negative evolution as

∂ t ϕ ≤ 0. ( 5
)
Next we follow [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations[END_REF] and [START_REF] Bourdin | A Variational Approach to the Numerical Simulation of Hydraulic Fracturing[END_REF] and replace energy (3) by a global constitutive dissipation functional for a rate independent fracture process

E ε (u, ϕ) = (0,L) 3 1 2 (1 -k)ϕ 2 + k Ge(u) : e(u) dx - ∂ N Ω τ • u dS- (0,L) 3 αϕ 1+b p B div u dx + G c (0,L) 3 1 2ε (1 -ϕ) 2 + ε 2 |∇ϕ| 2 dx, b ≥ 0. ( 6
)
We note that k > 0 is a regularization parameter with k ε.

Remark 1. A straightforward use of the phase field would suggest the 3rd integral in [START_REF] Bourdin | A Variational Approach to the Numerical Simulation of Hydraulic Fracturing[END_REF] to be of the form (0,L) 3 αϕp B div u dx; that is b = 0. In the classical case of elastic cracks one has 0 ≤ ϕ ≤ 1. We establish this property for the continuous time problem. Nevertheless, for the time discretized problem there will be no invariant region nor maximum principle estimates; the phase field unknown ϕ may be negative and take values larger than 1. Thus, we replace ϕ by ϕ + in terms where negative ϕ could lead to erroneous conclusions. In addition the entropy condition (5) prevents ϕ being larger than 1.

To ensure smoothness, we replace ϕ by ϕ 1+b , b > 0. We note that such a choice does not modify the problem for ϕ = 1 or ϕ = 0; however, the intermediate values are slightly smoothed.

Modeling the interaction between a crack C and ∂ N (0, L) 3 is not considered in this research and is excluded. The most complicated situation is in the region where Dirichlet and Neumann boundary conditions meet and where it is not clear how to define the phase field. We suppose a priori that the crack does not reach this region, i.e. that ϕ = 1 in a neighborhood O of the contact surface between outer boundary conditions.

In a porous medium, due to the presence of the pore structure, cracks are tiny but three dimensional bodies. The energy E ε contains only implicitly the presence of the pressurized crack in the term ∂ N Ω τ • u dS and still has to be written in an acceptable form which does not include ∂C. Since Ω = (0, L) 3 \ C and ∂Ω = ∂(0, L) 3 ∪ ∂C, we have ∂C ⊂ ∂ N Ω. The stress in the crack C is -p f I and at the crack boundary we have the continuity of the contact force

σn = (Ge(u) -αp B I)n = -p f n. (7)
The interface and corresponding notation is described in Figure 1. The pressure continuity at ∂C allows us to work with the pressure field p with p = p f in C and p = p B in Ω = (0, L) 3 \ C.

Ω (0, L) 3 with p f Ω with p B C zoom-in C ∂C
Before introducing the phase field we eliminate the implicit dependence of E ε on the crack by transforming the traction crack surface integrals as follows

Ω αp B div w dx + ∂C σnw dS = Ω (α -1)p B div w dx- Ω ∇p B w dx + ∂ N (0,L) 3 p B w • n dS.
After replacing G c H 2 (C) by the phase field regularization (4), the Fréchet derivative of the functional (6) with respect to u, gives the elasticity equation

(0,L) 3 (1 -k)ϕ 2 + k Ge(u) : e(w) dx - (0,L) 3 αϕ 1+b p B div w dx - ∂ N (0,L) 3 τ • w dS = - ∂C p f w • n dS = - (0,L) 3 p B div w dx - (0,L) 3 ∇p B • w dx + ∂ N (0,L) 3 p B w • n dS, for admissible w. ( 8
)
Next we introduce the phase field ϕ in the terms at the right hand side of [START_REF] Dean | Hydraulic-Fracture Predictions with a fully coupled geomechanical reservoir simulator[END_REF]. It yields

- (0,L) 3 p B div w dx → - (0,L) 3 ϕ 1+b p B div w dx and (0,L) 3 ∇p B • w dx → (0,L) 3 ϕ 1+b ∇p B • w dx.
After inserting the above transformations into (8) and replacing ϕ by ϕ + , we obtain the following elasticity equation

(0,L) 3 (1 -k)ϕ 2 + + k Ge(u) : e(w) dx - (0,L) 3 (α -1)ϕ 1+b + p div w dx+ (0,L) 3 ϕ 1+b + ∇p • w dx - ∂ N (0,L) 3
(τ + pn) • w dS = 0, for admissible w.

We choose as functional space of admissible displacements V U = {z ∈ H 1 ((0, L) 3 ) 3 | z = 0 on ∂ D (0, L) 3 }. Then, Equation (9) becomes Formulation 1 (Weak form of elasticity including pressure force).

(0,L) 3 (1 -k)ϕ 2 + + k Ge(u) : e(w) dx - (0,L) 3 (α -1)ϕ 1+b + p div w dx+ (0,L) 3 ϕ 1+b + ∇p • w dx - ∂ N (0,L) 3 (τ + pn) • w dS = 0, ∀w ∈ V U , (10)
or in differential form Formulation 2 (Differential form of elasticity including pressure force).

-div

(1 -k)ϕ 2 + + k Ge(u) + ϕ 1+b + ∇p + (α -1)∇ ϕ 1+b + p = 0 in (0, L) 3 , (11)
u = 0 on ∂ D (0, L) 3 , (12)
(1

-k)ϕ 2 + + k Ge(u)n = τ + αpn on ∂ N (0, L) 3 . ( 13
)
It remains to write the phase field equation. In differential form, the phase field equation reads Formulation 3 (Differential form of phase-field including pressure force).

∂ t ϕ ≤ 0 on (0, T ) × (0, L) 3 and ∂ϕ ∂n = 0 on (0, T ) × ∂(0, L) 3 ; ( 14)

-G c ε∆ϕ - G c ε (1 -ϕ) + (1 -k)Ge(u) : e(u)ϕ + + (1 + b)(1 -α)ϕ b + p div u + (1 + b)ϕ b + ∇p • u ≤ 0 in (0, T ) × (0, L) 3 , (15)
-G c ε∆ϕ - G c ε (1 -ϕ) + (1 -k)Ge(u) : e(u)ϕ + + (1 + b)(1 -α)ϕ b + p div u + (1 + b)ϕ b + ∇p • u ∂ t ϕ = 0 in (0, T ) × (0, L) 3 . (16)
In the above formulation, the inequality ( 15) is motivated by the entropy condition ∂ t ϕ ≤ 0. Therefore, (15) can be viewed as an equation only when ∂ t ϕ < 0. The Rice condition ( 16) is a complementarity condition, which states that either ∂ t ϕ < 0 (i.e. the fracture grows and ( 15) is a phase field equation) or ∂ t ϕ = 0 (i.e. the geometry does not change from the previous time step and we satisfy only the inequality ( 15)). The Rice condition is well-known in fracture mechanics and for incremental formulations it leads to the complementarity condition for the classical obstacle problem.

In order to write the variational form, we rewrite Rice's condition [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF].

Lemma 1. Let {u, ϕ} be smooth functions satisfying ( 11)- [START_REF] Mikelić | A phase field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF]. Then equality ( 16) is equivalent to

∂ t (0,L) 3 1 2 (1 -k)ϕ 2 + + k Ge(u) : e(u) + G c 2 ε|∇ϕ| 2 + 1 ε (1 -ϕ) 2 + ϕ 1+b + (1 -α)p div u + ∇p • u dx - ∂ N (0,L) 3 (τ + pn) • u dS - (0,L) 3 ϕ 1+b + (1 -α)∂ t p div u + ∇∂ t p • u dx + ∂ N (0,L) 3 (∂ t τ + ∂ t pn) • u dS = 0 on (0, T ). ( 17
)
Proof. Let us suppose equality [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF]. We integrate it with respect to x over (0, L) 3 . Next we use ∂ t u as a test function in equation ( 10) and we add the two equalities. It yields directly [START_REF] Mikelić | Phase-field modeling of pressurized fractures in a poroelastic medium[END_REF].

In the opposite direction, we subtract from equality [START_REF] Mikelić | Phase-field modeling of pressurized fractures in a poroelastic medium[END_REF] the equality [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF], with the test function w = ∂ t u. This subtraction yields equality (17) integrated over (0, L) 3 . Using ( 14)- [START_REF] Mikelić | A phase field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF] we conclude that (16) holds true.

Remark 2. In fact, we note that in ( 16) being equal to zero can be replaced by ≤ 0, which is used in the definition of the weak solution.

Let p 0 = p(•, 0), τ 0 = τ (•, 0) and ϕ 0 the initial value of ϕ. We denote by u 0 the solution for equation [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF] with ϕ = ϕ 0 , p = p 0 and τ = τ 0 . Before introducing the variational problem, we recall the definitions of the appropriate functional spaces. A function ψ ∈ L 1 ((0, L) 3 × (0, T )) is a function of bounded variation if and only if

V (ψ, (0, T ) × (0, L) 3 ) := sup T 0 (0,L) 3 ψ (div ζ + ∂ t ζ) : ζ ∈ C ∞ 0 ((0, T ) × (0, L) 3 ), ζ C 0 ≤ 1 < +∞. ( 18
) Moreover, V (ψ, (0, T ) × (0, L) 3 ) = |Dψ|((0, T ) × (0, L) 3 ). The norm ψ BV := ψ L 1 + V (ψ, (0, T ) × (0, L) 3 ) endows BV ((0, T ) × (0, L) 3
) with a Banach space structure.

A sequence {ψ n } converges weakly in BV ((0, T ) × (0, L) 3 ) to ψ if ψ n → ψ strongly in L 1 ((0, T ) × (0, L) 3 ) and ∇ψ n and ∂ t ψ n converge weakly in the sense of measures to ∇ψ and ∂ t ψ. In fact a sequence converges weakly if and only if it converges in L 1 and it is bounded in the BV norm (see Section 3.1 in [START_REF] Ambrosio | Functions of bounded variations and free discontinuity problems[END_REF]). The Sobolev inequality gives the embeddings BV ((0, T ) × (0, L) 3 ) ⊂ L q ((0, T ) × (0, L) 3 ) for every q ∈ [1, 4 3 ]. The embedding is compact for q < 4 3 (see Corollary 3.49 of [START_REF] Ambrosio | Functions of bounded variations and free discontinuity problems[END_REF]). The variational formulation, corresponding to the equations ( 14)-( 16) is:

Formulation 4 (Weak form of phase-field including pressure). Find ϕ ∈ BV ((0, T ) × (0, L) 3 ) ∩ L 2 (0, T ; H 1 ((0, L) 3 )) such that ∂ t ϕ, ∂ t ϕ -and ∂ t ϕ + are nonpositive bounded measures; ( 19
) T 0 (0,L) 3 (1 -k)ϕ + ψGe(u) : e(u) + G c - 1 ε (1 -ϕ)ψ + ε∇ϕ • ∇ψ dxdt +(1 + b) T 0 (0,L) 3 ϕ b + ((1 -α)p div u + ∇p • u) ψ dxdt ≤ 0, ∀ ψ ∈ L ∞ ((0, T ) × (0, L) 3 ) ∩ L 2 (0, T ; H 1 ((0, L) 3 )), ψ ≥ 0 a.e. on (0, T ) × (0, L) 3 ; (20) ϕ(x, 0) = ϕ 0 (x) on (0, L) 3 , 0 ≤ ϕ 0 (x) ≤ 1; (21)
(0,L) 3 1 2 (1 -k)ϕ 2 + (t) + k Ge(u(t)) : e(u(t)) + G c 2 ε|∇ϕ(t)| 2 + 1 ε (1 -ϕ(t)) 2 +ϕ 1+b + (t) (1 -α)p(t) div u(t) + ∇p(t) • u(t) dx - ∂ N (0,L) 3 (τ (t) + p(t)n) • u dS - (0,L) 3 1 2 (1 -k)(ϕ 0 + ) 2 + k Ge(u 0 ) : e(u 0 ) + G c 2 ε|∇ϕ 0 | 2 + 1 ε (1 -ϕ 0 ) 2 +(ϕ 0 + ) 1+b (1 -α)p 0 div u 0 + ∇p 0 • u 0 dx - ∂ N (0,L) 3 (τ 0 + p 0 n) • u 0 dS - t 0 (0,L) 3 ϕ 1+b + (1 -α)∂ η p div u + ∇∂ η p • u dxdη+ t 0 ∂ N (0,L) 3 (∂ η τ + ∂ η pn) • u dSdη ≤ 0 a. e. on (0, T ). ( 22
)
Remark 3. For a smooth solution ϕ such that ∂ t ∇ϕ -∈ L 2 ((0, T ) × (0, L) 3 ), nonpositivity of ∂ t ϕ, ∂ t ϕ -and ∂ t ϕ + implies nonnegativity of ϕ. This follows from the following observations. Namely, nonnegativity of -∂ t ϕ -and -∂ t ϕ + allows using them as test functions in [START_REF] Sneddon | The distribution of stress in the neighbourhood of a crack in an elastic solid[END_REF]. From Rice's equality and [START_REF] Sneddon | The distribution of stress in the neighbourhood of a crack in an elastic solid[END_REF] with -∂ t ϕ -and -∂ t ϕ + as respective test functions, we have

T 0 (0,L) 3 (1 -k)ϕ + ∂ t ϕ -Ge(u) : e(u) dxdt - G c ε T 0 (0,L) 3 (1 -ϕ)∂ t ϕ -dxdt+ G c ε T 0 (0,L) 3 ∇∂ t ϕ -• ∇ϕ dxdt + (1 + b) T 0 (0,L) 3 ϕ b + ((1 -α)p div u+ ∇p • u) ∂ t ϕ -dxdt = 0. ( 23
)
Since the supports of ϕ + and ϕ -are orthogonal, ϕ + (x, t)∂ t ϕ -(x, t) = 0 on (0, T ) × (0, L) 3 and the first and the fourth integral cancel. For the integrand in the second term we have

-(1 -ϕ)∂ t ϕ -= -∂ t ϕ -+ 1 2 ∂ t ϕ 2 - and ε T 0 (0,L) 3 ∇∂ t ϕ -• ∇ϕ dxdt = ε 2 (0,L) 3 |∇ϕ -(T )| 2 dx ≥ 0.
Insertion of the above equalities in equality [START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF] and using the nonnegativity assumption ϕ -(0) = 0 of the initial condition yield

- G c 2ε T 0 (0,L) 3 ∂ t ϕ - ≤0 dxdt + G c 2 (0,L) 3 1 ε ϕ 2 -(T ) + ε|∇ϕ -(T )| 2 dx = 0 (24)
Equality ( 24) yields ∂ t ϕ -= 0. Unfortunately our constructed solution is not necessarily smooth and this argument does not apply to weak solutions. Note that even for smooth ϕ, ϕ -is only a Lipschitz function and ∂ t ∇ϕ -in general does not exist as a function.

Our goal in the remainder of this paper is to consider the weak settings; namely Formulation 1 and Formulation 4. We notice that the corresponding Galerkin finite element approximations are provided in Finite Element Formulation 2 and Finite Element Formulation 1 presented in Section 4.

Existence of the quasistatic problem

In this section we prove existence of a solution for the quasistatic problem by introducing an incremental formulation.

We start by making the following assumptions on the data:

• G is a positive definite constant rank-4 tensor,

• k is a positive constant,

• and the following regularity assumptions:

p ∈ W 1,1 (0, T ; W 1,2r ((0, L) 3 )), r > 3; τ ∈ W 1,1 (0, T ; W 1-1/r,r (∂ N (0, L) 3 )),
and

ϕ 0 ∈ W 2,r ((0, L) 3 ), 0 ≤ ϕ 0 ≤ 1. ( 25
) Let K p = {ψ ∈ H 1 ((0, L) 3 )) | ψ ≤ ϕ p a.e. on (0, L) 3 }, (26)
where ϕ p is the phase field and u p the displacement from the previous time step, respectively. We consider a different incremental problem than introduced in [START_REF] Mikelić | A phase field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF]; namely, a decoupling into a linear elasticity problem and a nonlinear phase field problem. The nonlinear term ϕ b + is treated using a discrete derivative. In addition, the friction term is eliminated.

Let O be a small neighborhood of the contact surface between our outer boundary conditions and let χ be the indicator function of its complement. The size of O will be defined later precisely by (45).

The system takes the form:

(0,L) 3 χ(1 -k)ϕ + (ψ -ϕ)Ge(u p ) : e(u p ) dx + G c (0,L) 3 - 1 ε (1 -ϕ)(ψ -ϕ)+ ε∇ϕ • ∇(ψ -ϕ) dx + (0,L) 3 ϕ 1+b + -(ϕ p ) 1+b + ϕ -ϕ p (1 -α)p p div u p + ∇p p • u p (ψ -ϕ) dx ≥ 0, ∀ ψ ∈ K p ∩ L ∞ ((0, L) 3 ), ( 27
) (0,L) 3 (1 -k)ϕ 2 + + k Ge(u) : e(w) dx - ∂ N (0,L) 3 (τ p + p p n) • w dS + (0,L) 3 ϕ 1+b + (1 -α)p p div w + ∇p p • w dx = 0, ∀w ∈ V U . ( 28
)
We briefly explain the rationale of ( 27) and in particular the replacing of (1 + b)ϕ b + by the finite difference

ϕ 1+b + -(ϕ p ) 1+b + ϕ -ϕ p . ( 29
)
We firste note that in the quasistatic limit ϕ is approaching ϕ p and the quotient from (29) behaves as (1 + b)ϕ b + . Next we take as test function ψ = ϕ p -ψ in (20) and we recall that the time derivative difference quotient only needs to approximatively satisfy the Rice condition. For a second test function ϕ p -ϕ we impose the equality in [START_REF] Sneddon | The distribution of stress in the neighbourhood of a crack in an elastic solid[END_REF]. Subtracting the inequality and equality yields (27). Regarding the finite difference (29), we remark that several estimates are only valid with this choice rather than the original term. Only this choice allows us to pass to the quasi-static limit, which was not possible for the implicit formulation as discussed in [START_REF] Mikelić | A phase field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF][START_REF] Mikelić | Phase-field modeling of pressurized fractures in a poroelastic medium[END_REF].

The phase field step

Let δ > 0 and θ δ be given by

θ δ (y) =    1, y ≤ 0; 1 -y/δ, 0 < y ≤ δ; 0, δ < y. (30)
Next let φ = inf{1, ϕ + }, ϕ ∈ L 1 ((0, L) 3 ), and let

k reg (ϕ p , g) = inf{0, -G c ε∆ϕ p - G c ε (1 -ϕ p ) + (1 -k)χGe(u p ) : e(u p )g+ inf{1, g 1+b + -(ϕ 1+b p ) + g -ϕ p } (1 -α)p p div u p + ∇p p • u p } on (0, L) 3 . ( 31
)
We consider the following penalized variant of the variational inequality (27

) (0,L) 3 (1 -k)χ φψGe(u p ) : e(u p ) dx + G c (0,L) 3 - 1 ε (1 -ϕ)ψ+ ε∇ϕ • ∇ψ dx + (0,L) 3 inf{1, ϕ 1+b + -(ϕ p ) 1+b + ϕ -ϕ p } (1 -α)p p div u p + ∇p p • u p ψ dx - (0,L) 3 k reg (ϕ p , ϕ)θ δ (ϕ p -ϕ)ψ = 0, ∀ ψ ∈ H 1 ((0, L) 3 ), (32)
In the following propositions we drop the domain notation for the function spaces.

Proposition 1. Let p p ∈ W 1,r , ϕ p ∈ W 2,r , ϕ p ≤ 1 and u p ∈ W 1,2r , r > 3. Then there exists a solution ϕ = ϕ δ ∈ H 1 for variational equation (32) satisfying

||ϕ δ || W 2,r ≤ C(||u p || 2 W 1,2r + ||p p || 2 W 1,r ), ( 33
)
where the constant C is independent of δ.

Proof. See Appendix A.

Proposition 2. Let p p ∈ W 1,r , ϕ p ∈ W 2,r , ϕ p ≤ 1 and u p ∈ W 1,2r , r > 3. Then ϕ δ ϕ weakly in W 2,r , as δ → 0, where ϕ ∈ W 2,r is a solution of the variational inequality (27) satisfying estimate (33).

Proof. We first prove that ϕ δ ∈ K p . Let ζ = ϕ δ -inf{ϕ δ , ϕ p } ≥ 0. We use ζ as a test function in equation ( 32), with the goal to prove that ζ = 0. Testing equation (32) with ζ yields

G c (0,L) 3 1 ε (ϕ δ -ϕ p )ζ + ε∇(ϕ δ -ϕ p ) • ∇ζ dx + (0,L) 3 (1 -k)χ φδ Ge(u p ) : e(u p ) + inf{1, (ϕ 1+b δ ) + -(ϕ p ) 1+b + ϕ δ -ϕ p } (1 -α)p p div u p + ∇p p • u p -G c ε∆ϕ p - G c ε (1 -ϕ p ) ζ dx - (0,L) 3 k reg (ϕ p , ϕ δ )θ δ (ϕ p -ϕ δ )ζ = 0. ( 34
) For ϕ δ < ϕ p , we have ζ = 0. Next, ϕ δ ≥ ϕ p yields θ δ = 1 and ζ = ϕ δ -ϕ p . Finally, inf{1, (ϕ δ ) 1+b + -(ϕ p ) 1+b + ϕ δ -ϕ p } (1 -α)p p div u p + ∇p p • u p + (1 -k)χ φδ Ge(u p ) : e(u p ) -G c ε∆ϕ p - G c ε (1 -ϕ p ) -k reg (ϕ p , ϕ δ ) ≥ 0 (a.e) on (0, L) 3
and we obtain from (34) ζ = 0. By estimate (33), the set {ϕ δ } δ>0 is bounded in W 2,r independently of δ. Hence, by the weak compactness, it contains a subsequence in K p which converge weakly in W 2,r and strongly in C 1,s , s < 1 -3/r, to an element ϕ of K p .

To show that ϕ satisfies equation ( 27), it is enough to follow [START_REF] Kinderlehrer | An introduction to variational inequalities and their applications[END_REF], page 109 and apply Minty's lemma to the monotone term defined by k reg θ δ .

The elasticity step

Proposition 3. Let p ∈ W 1,2r , τ ∈ W 1-1/r,r and ϕ ∈ W 2,r , r > 3.
Let O be a smooth neighborhood of the contact surface between Dirichlet and Neumann conditions. Then there exists a unique solution u ∈ W 2,r ((0, L) 3 \ O) 3 for variational equation (28).

Proof. Obviously, problem (28) has a unique solution u ∈ H 1 ((0, L) 3 ) 3 . We write in the following differential form

-div Ge(u) = Ge(u)∇ log (1 -k)ϕ 2 + + k - ϕ 1+b + ∇p p + (α -1)∇ ϕ 1+b + p p (1 -k)ϕ 2 + + k in (0, L) 3 , (35)
u = 0 on ∂ D (0, L) 3 , ( 36
) (1 -k)ϕ 2 + + k Ge(u)n = τ p + αp p n on ∂ N (0, L) 3 . ( 37
)
Then the regularity theory for the Navier equations of the linear elasticity (see e.g. [START_REF] Ciarlet | Three-dimensional elasticity[END_REF], Theorems 6.3-6) and a simple bootstrapping yield u ∈ W 2,r ((0, L) 3 \ O) 3 .

The quasistatic limit

In this subsection we suppose for simplicity b = 1. We start by stating auxiliary lemmas easily derived from elementary inequalities:

Lemma 2.

(

-k)ϕ + (ϕ -ϕ p )Ge(u p ) : e(u p ) + (1 -k)ϕ 2 + + k Ge(u) : e(u -u p ) ≥ 1 2 (1 -k)(ϕ 2 + -(ϕ p ) 2 + )Ge(u p ) : e(u p ) + (1 -k)ϕ 2 + + k (Ge(u) : e(u)- Ge(u p ) : e(u p ) = 1 2 (1 -k)ϕ 2 + + k Ge(u) : e(u)- (1 -k)(ϕ p ) 2 + + k Ge(u p ) : e(u p ) . ( 1
) Lemma 3. - 1 ε (1 -ϕ)(ϕ -ϕ p ) + ε∇ϕ • ∇(ϕ -ϕ p ) ≥ - 1 ε (ϕ -ϕ p ) + 1 2ε (ϕ 2 -ϕ 2 p ) + ε 2 (|∇ϕ| 2 -|∇ϕ p | 2 ). ( 38
) 39
Lemma 4.

(ϕ 2 + -(ϕ p ) 2 + )p p div u p + ϕ 2 + p p div (u -u p ) = p div u ϕ 2 + -p p div u p (ϕ p ) 2 + -ϕ 2 + div u(p -p p ), ( 40
) (ϕ 2 + -(ϕ p ) 2 + )∇p p • u p + ϕ 2 + ∇p p • (u -u p ) = ∇p • u ϕ 2 + -∇p p • u p (ϕ p ) 2 + -ϕ 2 + (∇p -∇p p ) • u. (41)
We suppose the quasistatic problem [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF], ( 20)-( 22) is discretized with a uniform time step ∆t. Given solutions at discrete times t j , j = 0, . . . , N , {ϕ ∆t (t j ), u ∆t (t j )} are extended from the discrete times {t j } 0≤j≤N to (0, T ) by ϕ ∆t (t) = ϕ ∆t (t j ) if t j ≤ t < t j+1 , j = 0, . . . N -1;

(42)

ϕ ∆t (t 0 ) = ϕ ∆t (0) = ϕ 0 ; (43) u ∆t (t) = u ∆t (t j ) if t j ≤ t < t j+1 , j = 0, . . . N -1. ( 44
)
Proposition 4. Let us suppose that hypothesis (25) holds true and

|O| γ/(2+γ) /|∆t| = O(|∆t|). ( 45
)
Let {ϕ ∆t , u ∆t } be a solution to (27)-(28) corresponding to the time discretization step ∆t. Then we have

||(ϕ ∆t ) + e(u ∆t )|| L ∞ (0,T ;L 2 ) + √ k||u ∆t || L ∞ (0,T ;H 1 ) ≤ C, (46)
||(ϕ ∆t ) + || L ∞ ((0,T )×(0,L) 3 ) ≤ C, (47)
||ϕ ∆t || L ∞ (0,T ;L 2 ) + ||ε∇ϕ ∆t || L ∞ (0,T ;L 2 ) + ||∂ ∆t ϕ ∆t || L 1 ((0,T )×(0,L) 3 ) ≤ C, ( 48
)
where C is a generic constant independent of ∆t and

∂ ∆t ψ(t) = ψ(t j+1 ) -ψ(t j ) ∆t
, for t j ≤ t < t j+1 , j = 0, . . . N -1.

Proof. First we recall that initial condition [START_REF] Sneddon | Crack problems in the classical theory of elasticity[END_REF] and definition of the convex sets K p imply ϕ ∆t ≤ 1 on (0, L) 3 × (0, T ) and 0 ≤ (ϕ ∆t ) + ≤ 1 on (0, L) 3 × (0, T ). ( 49)

Estimate (49) contains estimate (47).

Next we use w = u ∆t as a test function in variational equation (28). It yields

(0,L) 3 |(ϕ ∆t ) + e(u ∆t )| 2 + k(|∇u ∆t | 2 + |u ∆t | 2 ) dx ≤ C.
The above inequality implies estimate (46), with constant C independent of ∆t. At this point we recall a Meyers type result for the linear elasticity nonhomogeneous Navier equations from [START_REF] Shi | Higher integrability of the gradient in linear elasticity[END_REF]. It says that there is a γ > 0 such that,

||u ∆t (t)|| W 1,2+γ ≤ C(k, ||(ϕ ∆t ) + || ∞ , L, G) ||(ϕ ∆t ) 2 + p(t)|| L 2+γ ((0,L) 3 ) + ||τ (t) + αp(t)|| L 2(2+γ)/3 (∂ N (0,L) 3 ) . (50)
Hence under the data smoothness we have, in addition to (46), that

||u ∆t || L ∞ (0,T ;W 1,2+γ ) ≤ C, (51)
where C does not depend on the time discretization step ∆t. Consequently sup 0≤t≤T (0,L) 3

(1 -χ)Ge(u ∆t ) : e(u ∆t ) ≤ C|O| γ/(2+γ) .

Next, for j ∈ {1, . . . , N } we set u = u ∆t (t j ) and u p = u ∆t (t j-1 ) and test equation ( 28) by w = uu p . Then we subtract inequality (27), with ψ = ϕ p , from the obtained equality. It yields

(0,L) 3 (1 -k)ϕ + (ϕ -ϕ p )Ge(u p ) : e(u p ) + (1 -k)ϕ 2 + + k Ge(u) : e(u -u p ) dx +G c (0,L) 3 - 1 ε (1 -ϕ)(ϕ -ϕ p ) + ε∇ϕ • ∇(ϕ -ϕ p ) dx + (0,L) 3 (1 -α) (ϕ 2 + -(ϕ p ) 2 + )p p div u p + ϕ 2 + p p div (u -u p ) + (ϕ 2 + -(ϕ p ) 2 + )∇p p • u p + ϕ 2 + ∇p p • (u -u p ) dx - ∂ N (0,L) 3 (τ p + p p n) • (u -u p ) dS = - (0,L) 3 (1 -k)(1 -χ)ϕ + (ϕ -ϕ p )Ge(u p ) : e(u p ) dx. (53)
After inserting (38), (39), ( 40) and ( 41) into (53), we get

(0,L) 3 1 2 (1 -k)ϕ 2 + + k Ge(u) : e(u) -(1 -k)(ϕ p ) 2 + + k Ge(u p ) : e(u p ) dx -G c (0,L) 3 ϕ -ϕ p ε dx + G c (0,L) 3 1 2ε (ϕ 2 -ϕ 2 p ) + ε 2 (|∇ϕ| 2 -|∇ϕ p | 2 ) dx+ (0,L) 3 (1 -α)p div u + ∇p • u ϕ 2 + -((1 -α)p p div u p + ∇p p • u p (ϕ p ) 2 + dx+ ∂ N (0,L) 3 (τ + pn) • u -(τ p + p p n) • u p dS ≤ ∆t ∂ N (0,L) 3 |u|| τ -τ p + (p -p p )n ∆t | dS +C∆t (0,L) 3 (|ϕ + e(u)|| p -p p ∆t | + |u||∇ p -p p ∆t |) dx + C|O| γ/(2+γ) . ( 54
)
We remark that our motivation for using the discrete derivative for ϕ 2 + , instead of 2ϕ + , is to obtain (54). Namely, the Taylor remainder can not be controlled by our time estimates. Next we sum up over the time intervals (t j-1 , t j ), j = 1, . . . , M and obtain

(0,L) 3 1 2 (1 -k)(ϕ ∆t ) 2 + (t) + k Ge(u ∆t (t)) : e(u ∆t (t)) dx +G c t 0 (0,L) 3 |∂ ∆t ϕ ∆t (ξ)| ε dxdξ +G c (0,L) 3 ϕ 2 ∆t (t) 2ε + ε 2 |∇ϕ ∆t (t)| 2 dx ≤ G c (0,L) 3 ϕ 2 (0) 2ε + ε 2 |∇ϕ(0)| 2 dx +C||p|| W 1,1 (0,t;L 2 ) ||(ϕ ∆t ) + e(u ∆t )|| L ∞ (0,t;L 2 ) + C||∇p|| W 1,1 (0,t;L 2 ) ||u ∆t || L ∞ (0,t;L 2 ) + C(||τ || W 1,1 (0,t;L 2 (∂ N (0,L) 3 )) + ||p|| W 1,1 (0,t;L 2 (∂ N (0,L) 3 )) )||u ∆t || L ∞ (0,t;H 1 ) + C |∆t| |O| γ/(2+γ) + (0,L) 3 (1 -k)ϕ 2 (0) + k Ge(u(0)) : e(u(0)) dx, (55)
for t ≤ t M . Using that |O| γ/(2+γ) /|∆t| = O(|∆t|), estimate (48) follows from (55).

Our goal is now to use the estimates ( 46)-(48) and pass to the limit ∆t → 0.

Theorem 1. Let us suppose the data regularity assumptions ( 25) and (45). Let {ϕ ∆t , u ∆t } be a solution to ( 27)-( 28) corresponding to the time discretization step ∆t. Then there is a subsequence of {ϕ ∆t , u ∆t }, denoted by the same subscript, and {ϕ, u} ∈ BV ((0,

T ) × (0, L) 3 ) ∩ L 2 (0, T ; H 1 ) × L ∞ (0, T ; H 1 ), ϕ ≤ ϕ 0 a.e. on (0, T )×(0, L) 3 , ∂ t ϕ, ∂ t ϕ -and ∂ t ϕ + are nonpositive bounded measures, such that u ∆t * u weak - * in L ∞ (0, T ; H 1 ); ( 56
)
u ∆t → u strongly in L 2 (0, T ; H 1 ); ( 57
)
ϕ ∆t * ϕ weak - * in BV ((0, T ) × (0, L) 3 ); ( 58
)
ϕ ∆t → ϕ strongly in L q ((0, T ) × (0, L) 3 ), ∀q ∈ [1, 4 3
);

(59)

ϕ ∆t ϕ weakly in L 2 (0, T ; H 1 ), ( 60
)
when ∆t → 0. Furthermore, {ϕ, u} is a solution for ( 10), ( 20)- [START_REF] Tolstoy | Acoustics, elasticity, and thermodynamics of porous media. Twenty-one papers[END_REF].

Proof. We start by observing that

< ∂ t ϕ ∆t , ψ >=< ∂ ∆t ϕ ∆t , ψ >, ∀ψ ∈ C ∞ 0 ((0, T ) × (0, L) 3 ),
and by density for all ψ ∈ C 0 ((0, T ) × (0, L) 3 ). Therefore, we have

| < ∂ t ϕ ∆t , ψ > | ≤ ||∂ ∆t ϕ ∆t || L 1 ||ψ|| L ∞ , ∀ψ ∈ C 0 ((0, L) 3 × (0, T )).
Estimates ( 46)-( 48), and the above estimate imply the a priori estimates

||u ∆t || L ∞ (0,T ;H 1 ) + ||ϕ ∆t || BV ((0,T )×(0,L) 3 ) + ||ϕ ∆t || L 2 (0,T ;H 1 ) ≤ C. (61)
From ( 61) we observe that Sobolev embeddings and the weak * compactness in L ∞ (0, T ; H 1 ) and in BV ((0, T ) × (0, L) 3 ) ∩ L 2 (0, T ; H 1 ) respectively, take place. This provides existence of a subsequence of {ϕ ∆t , u ∆t }, denoted by the same subscript, with limit {ϕ, u} ∈ BV ((0, T )×(0, L) 3 )∩L 2 (0, T ; H 1 )× L ∞ (0, T ; H 1 ), satisfying convergences ( 56), ( 58)-( 60). In addition, ϕ ≤ ϕ 0 a.e. on (0, T ) × (0, L) 3 and ∂ t ϕ is a nonpositive bounded measure. Since u ∆t (t) ≤ u ∆t (t -∆t) implies the same inequality for its positive and negative parts, ∂ t ϕ -and ∂ t ϕ + are also nonpositive bounded measures. It remains to prove that {ϕ, u} is a solution for ( 10), ( 20)-( 21), satisfies inequality ( 22) and u ∆t → u strongly in L 2 (0, T ; H 1 ). Replacing in (28) u by u ∆t (t) and ϕ by ϕ ∆t (t) and passing to the limit, we obtain

T 0 (0,L) 3 (1 -k)ϕ 2 + + k Ge(u) : e(w) dxdt- T 0 (0,L) 3 (α -1)ϕ 2 + p div w dxdt + T 0 (0,L) 3 ϕ 2 + ∇pw dxdt - T 0 ∂ N (0,L) 3 (τ + pn) • w dSdt = 0, ∀w ∈ L 2 (0, T ; V U ). ( 62
)
Next we choose w = u ∆t as test function in (28) and pass to the limit ∆t → 0. It yields

lim ∆t→0 T 0 (0,L) 3 (1 -k)(ϕ ∆t ) 2 + + k Ge(u ∆t ) : e(u ∆t ) dxdt = T 0 (0,L) 3 (1 -k)(ϕ + ) 2 + k Ge(u) : e(u) dxdt. ( 63
)
Using Fatou's lemma we have

T 0 (0,L) 3 lim inf ∆t→0 (1 -k)(ϕ ∆t ) 2 + + k Ge(u ∆t ) : e(u ∆t ) dxdt ≤ lim inf ∆t→0 T 0 (0,L) 3 (1 -k)(ϕ ∆t ) 2 + + k Ge(u ∆t ) : e(u ∆t ) dxdt = T 0 (0,L) 3 (1 -k)(ϕ) 2 + + k Ge(u) : e(u) dxdt. (64)
Consequently,

u ∆t → u strongly in L 2 (0, T ; V U ), as ∆t → 0. ( 65
)
We note that u ∈ L ∞ (0, T ; W 1,2+γ ), for some γ > 0, (66) and for every nonnegative function g ∈ L 1+γ 0 ((0, T ) × (0, L) 3 ), γ 0 > 0,

T -∆t 0 (0,L) 3 |(ϕ ∆t ) + (t + ∆t) -(ϕ ∆t ) + (t)|g dxdt ≤ C|∆t| γ 0 /(1+γ 0 ) ||g|| L 1+γ 0 ((0,T )×(0,L) 3 ) ||∂ ∆t ϕ ∆t || (1+γ 0 )/γ 0 L 1 ((0,T )×(0,L) 3 ) . ( 67
)
For every ψ ∈ C ∞ 0 ((0, L) 3 × (0, T )), (65) implies

lim ∆t→0 | T -∆t 0 (0,L) 3 (ϕ ∆t ) + (t + ∆t)ψGe(u ∆t -u) : e(u ∆t -u) dxdt| → 0, as ∆t → 0, (68)
and ( 63)-( 68) yield Next we use that

T -∆t 0 (0,L) 3 (ϕ ∆t ) + (t + ∆t)ψGe(u ∆t ) : e(u ∆t ) dxdt = T -∆t 0 (0,L) 3 (ϕ ∆t ) + (t + ∆t)ψGe(u ∆t -u) : e(u ∆t -u) dxdt+
| (ϕ ∆t ) 2 + (t) -(ϕ 2 ∆t ) + (t -∆t) ϕ ∆t (t) -ϕ ∆t (t -∆t) -2(ϕ ∆t ) + (t)| ≤ C|(ϕ ∆t ) + (t) -(ϕ ∆t ) + (t -∆t)|
and convergences ( 63)-( 68) to get

lim ∆t→0 T ∆t (0,L) 3 (ϕ ∆t ) 2 + (t) -(ϕ 2 ∆t ) + (t -∆t) ϕ ∆t (t) -ϕ ∆t (t -∆t) (1 -α)p(t -∆t) div u(t -∆t)+ ∇p(t -∆t) • u(t -∆t) ψ dxdt = lim ∆t→0 T 0 (0,L) 3 2(ϕ ∆t ) + (t) (1 -α)p(t) div u(t) + ∇p(t) • u(t) ψ dxdt = T 0 (0,L) 3 2ϕ + (t) (1 -α)p(t) div u(t) + ∇p(t) • u(t) ψ dxdt. ( 70
)
Now we write equation ( 27) in the equivalent form (71)-( 72)

(0,L) 3 χ(1 -k)(ϕ ∆t ) + (t)ψGe(u ∆t )(t -∆t) : e(u ∆t (t -∆t)) dx +G c (0,L) 3 - 1 ε (1 -ϕ ∆t (t))ψ + ε∇ϕ ∆t (t) • ∇ψ dx+ (0,L) 3 (ϕ ∆t ) 2 + (t) -(ϕ ∆t ) 2 + (t -∆t) ϕ ∆t (t) -ϕ ∆t (t -∆t) (1 -α)p(t -∆t) div u ∆t (t -∆t)+ ∇p(t -∆t) • u ∆t (t -∆t) ψ dx ≤ 0, a.e. on (∆t, T ), ∀ ψ ∈ C ∞ ([0, L] 3 ), ψ ≥ 0 a.e. on (0, L) 3 ; ( 71
) (0,L) 3 χ(1 -k)(ϕ ∆t ) + (t)(ϕ ∆t (t -∆t) -ϕ ∆t (t))Ge(u ∆t (t -∆t)) : e(u ∆t (t- ∆t)) dx + G c (0,L) 3 - 1 ε (1 -ϕ ∆t )(ϕ ∆t (t -∆t) -ϕ ∆t (t))+ ε∇ϕ ∆t • ∇(ϕ ∆t (t -∆t) -ϕ ∆t (t)) dx + (0,L) 3 (ϕ ∆t ) 2 + (t) -(ϕ ∆t ) 2 + (t -∆t) ϕ ∆t (t) -ϕ ∆t (t -∆t) ( 1 -α)p(t -∆t) div u ∆t (t -∆t)+ ∇p(t -∆t) • u ∆t (t -∆t) (ϕ ∆t (t -∆t) -ϕ ∆t (t)) dx = 0, a.e. on (∆t, T ). (72)
Using convergences (56)-( 60), ( 65) and ( 69)-(70) we pass to the limit ∆t → 0 in inequality (71) and get

T 0 (0,L) 3 (1 -k)ϕ + ψGe(u) : e(u) dxdt + G c T 0 (0,L) 3 - 1 ε (1 -ϕ)ψ +ε∇ϕ • ∇ψ dxdt + 2(1 -α) T 0 (0,L) 3 ϕ + p div uψ dxdt+ T 0 (0,L) 3 2ϕ + ∇p • u ψ dxdt ≤ 0, ∀ ψ ∈ C ∞ ([0, T ] × [0, L] 3 ), ψ ≥ 0 a.e. on (0, T ) × (0, L) 3 ; (73) ϕ(x, 0) = ϕ 0 (x) on (0, L) 3 , 0 ≤ ϕ 0 (x) ≤ 1. (74)
It remains to prove Rice's equality in its weak form [START_REF] Tolstoy | Acoustics, elasticity, and thermodynamics of porous media. Twenty-one papers[END_REF]. We write estimate (54) in the form

(0,L) 3 1 2 (1 -k)ϕ 2 + + k Ge(u) : e(u) -(1 -k)(ϕ p ) 2 + + k Ge(u p ) : e(u p ) dx +G c (0,L) 3 1 2ε ((1 -ϕ) 2 -(1 -ϕ p ) 2 ) + ε 2 (|∇ϕ| 2 -|∇ϕ p | 2 ) dx+ (0,L) 3 (1 -α)p div u + ∇p • u ϕ 2 + -((1 -α)p p div u p + ∇p p • u p (ϕ p ) 2 + dx- ∂ N (0,L) 3 (τ + pn) • u -(τ p + p p n) • u p dS + ∆t ∂ N (0,L) 3 u • τ -τ p + (p -p p )n ∆t dS -∆t (0,L) 3 ϕ 2 + div u p -p p ∆t + u • ∇ p -p p ∆t ) dx = C|O| γ/(2+γ) (75)
After summing up over the time intervals (t j-1 , t j ) we obtain

(0,L) 3 1 2 (1 -k)ϕ 2 + (t) + k Ge(u(t)) : e(u(t)) dx + G c (0,L) 3 (1 -ϕ(t)) 2 2ε + ε 2 |∇ϕ(t)| 2 dx + (0,L) 3 (1 -α)p div u(t) + ∇p • u(t) ϕ 2 + (t) dx- ∂ N (0,L) 3 (τ + pn) • u(t) dS - t 0 (0,L) 3 ϕ 1+b + (1 -α)∂ t p div u + ∇∂ t p • u dxdt+ t 0 ∂ N (0,L) 3 (∂ t τ + ∂ t pn) • u dSdt - (0,L) 3 1 2 (1 -k)(ϕ 0 + ) 2 + k Ge(u 0 ) : e(u 0 )+ G c 2 ε|∇ϕ 0 | 2 + 1 ε (1 -ϕ 0 ) 2 + (ϕ 0 + ) 2 (1 -α)p 0 div u 0 + ∇p 0 • u 0 dx + ∂ N (0,L) 3 (τ 0 + p 0 n) • u 0 dS = C ∆t |O| γ/(2+γ) , (76)
for all t ∈ (0, T ). Next we multiply equality (76) by the characteristic function of any time interval and pass to the limit ∆t → 0. Using the previously established convergences and the lower semicontinuity of the L 2 -norm of ∇ϕ ∆t , we obtain inequality [START_REF] Tolstoy | Acoustics, elasticity, and thermodynamics of porous media. Twenty-one papers[END_REF].

Numerical approximation

In this section, we formulate finite element approximations for the Formulations 4 and 1 presented in Section 2. Specifically, we apply a standard Galerkin finite element method in 2D on quadrilaterals.

The displacements u are approximated by continuous bilinears and are referred to as the finite element space V h . We take ϕ to be continuous bilinear and denote this space as W h . The standard spatial approximation parameter is represented by h.

In our numerical treatment, we set b = 1. The regularized incremental problem reads:

(0,L) 3 (1 -k)ϕ h + ψGe(u p ) : e(u p ) dx + G c (0,L) 3 - 1 ε (1 -ϕ h )ψ + ε∇ϕ h • ∇ψ dx + (0,L) 3 2 (1 -α)ϕ h + p p div u p + ϕ h + ∇p p • u p ψ dx - (0,L) 3 k reg (ϕ p , ϕ h )θ δ (ϕ p -ϕ h )ψ = 0, ∀ψ ∈ W h , (77)
(0,L) 3 (1 -k)(ϕ h + ) 2 + k Ge(u h ) : e(w) dx - ∂ N (0,L) 3 (τ p + p p n) • w dS + (0,L) 3 (ϕ h ) 2 + (1 -α)p p div w + ∇p p • w dx = 0, ∀w ∈ V h . (78)
In contrast to [START_REF] Mikelić | A phase field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF] (and also [START_REF] Heister | A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF]), we use a sequential coupling algorithm in which both subproblems are solved subsequently; namely, we first solve equation (77) for ϕ h and then solve equation (78) for u h , with a given phase field ϕ h . This procedure has two advantages:

• it is computationally efficient since we use a loosely coupled scheme, which only requires a few subiterations;

• considering the fact that the elasticity problem is linear no Newton method is required.

In the numerical examples in Section 5, we deal with an isotropic poroelastic medium with

G ijkl = µ(δ il δ jk + δ ik δ jl ) + λδ ij δ kl .
The stress tensor is defined as

σ ij = k,l
G ijkl e kl (u) = µe ij (u) + λtr(e)δ ij , where µ and λ denote the Lamé parameters. Our computations are compared with the approach from [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF], where Hooke's law is modified by introducing an anisotropic energy storage functional. Here, the stress tensor is additively decomposed into a tensile part σ + and a compressive part σ -, i.e., σ := (1 -k)ϕ 2 + k σ + + σ -. We emphasize that the energy degradation only acts on the tensile part. The modified energy functional then reads:

E ε (u, ϕ) = (0,L) 3 1 2 (1 -k)ϕ 2 + + k σ + : e(u) + σ -: e(u) dx - ∂ N Ω τ • u dS- (0,L) 3 αϕ 2 + p B div u dx + G c (0,L) 3 1 2ε (1 -ϕ) 2 + ε 2 |∇ϕ| 2 dx. (79)
Here, the two stress contributions are given by: The corresponding Euler-Lagrange equations read:

Formulation 5 (Weak form of phase-field including pressure and an anisotropic energy storage functional).

(0,L) 3

(1 -k)ϕ h + ψσ + : e(u p ) dx + G c (0,L) 3 - 1 ε (1 -ϕ h )ψ + ε∇ϕ h • ∇ψ dx + (0,L) 3 2 (1 -α)ϕ h + p p div u p + ϕ h + ∇p p • u p ψ dx - (0,L) 3 k reg (ϕ p , ϕ h )θ δ (ϕ p -ϕ h )ψ = 0, ∀ψ ∈ W h , (80)
Formulation 6 (Weak form of elasticity including pressure and an anisotropic energy storage functional).

(0,L) 3

(1 -k)(ϕ h + ) 2 + k σ + : e(w) dx + (0,L) 3

σ -: e(w) dx

- ∂ N (0,L) 3 (τ p + p p n) • w dS + (0,L) 3 (ϕ h ) 2 + (1 -α)p p div w + ∇p p • w dx = 0, ∀w ∈ V h . (81)
The isotropic system (77) and ( 78) is recovered by setting σ + := Ge(u p ) and σ -≡ 0.

We formulate separately a semilinear and bilinear form for each of the two subproblems. The spatially discretized elasticity problem can be written in the following way:

Finite Element Formulation 1 (Variational FE formulation of the nonlinear phase-field function).

Given p p and u p find ϕ h ∈ W h such that:

A(ϕ h )(ψ) = (0,L) 3 (1 -k)ϕ h + ψσ + : e(u p ) dx + G c (0,L) 3 - 1 ε (1 -ϕ h )ψ + ε∇ϕ h • ∇ψ dx + (0,L) 3 2 (1 -α)ϕ h + p p div u p + ϕ h + ∇p p • u p ψ dx - (0,L) 3 k reg (ϕ p , ϕ h )θ δ (ϕ p -ϕ h )ψ = 0, ∀ψ ∈ W h .
Finite Element Formulation 2 (Variational FE formulation of elasticity). Given p p and ϕ p , find u h ∈ V h such that:

B(u h )(w) = (0,L) 3 (1 -k)(ϕ h + ) 2 + k σ + : e(w) dx + (0,L) 3 σ -: e(w) dx - ∂ N (0,L) 3 (τ p + p p n) • w dS + (0,L) 3 (ϕ + ) 2 (1 -α)p p div w + ∇p p • w dx = 0, ∀w ∈ V h .
We adopt a partitioned coupling scheme to solve the equations. First, the nonlinear problem (Finite Element Formulation 1) is solved with Newton's method. For the iteration steps m = 0, 1, 2, . . ., we have:

A (ϕ h,m )(δϕ h , ψ) = -A(ϕ h,m )(ψ), ϕ h,m+1 = ϕ h,m + ωδϕ h , (82)
with a line search parameter ω ∈ (0, 1]. Here, the Jacobian of A(ϕ h )(ψ) is denoted by A (•)(•). Then, we solve for the linear elasticity problem (Finite Element Formulation 2), which can be treated with an appropriate solver.

Numerical tests

We perform four numerical tests where we assume a dimensionless form of the equations. In the first example, we consider Sneddon's benchmark and compare it to findings in the literature. In the second and third test cases, we compare formulation (77), ( 78) with (80), (81). In the second test case, an increasing pressure is applied to a single crack, which starts propagating after some time. For the third test, we use an increasing pressure to study two joining cracks. In the final test, heterogeneous porous media are considered. The last two tests demonstrate the performance of the variational approach to fractures in which joining and branching of cracks are easily treated. The programming code is a modification of the multiphysics program template [START_REF] Wick | Solving monolithic fluid-structure interaction problems in arbitrary Lagrangian Eulerian coordinates with the deal.ii library[END_REF], based on the finite element software deal.II [START_REF] Bangerth | Differential Equations Analysis Library[END_REF].

Sneddon's benchmark

The first example is motivated by Bourdin et al. [START_REF] Bourdin | A Variational Approach to the Numerical Simulation of Hydraulic Fracturing[END_REF] and is based on Sneddon's theoretical calculations [START_REF] Sneddon | Crack problems in the classical theory of elasticity[END_REF][START_REF] Sneddon | The distribution of stress in the neighbourhood of a crack in an elastic solid[END_REF]. Specifically, we consider a 2D problem where a (constant) pressure p is used to drive the deformation and crack propagation. (required for computing the width of the crack), mesh, and zoom-in of the mesh with the crack pattern. The crack is denoted in red with ϕ = 0 and the unbroken material is blue with ϕ = 1. In the second figure, the normal displacements are positive (red) above the fracture and negative (blue) below the fracture. The displacements are of order 10 -4 (see Figure 3). The green part in the final figure shows the thickness ε of the mushy zone in the phase-field variable. The configuration is displayed in Figure 2. Here we have assumed the following geometric data: Ω = (0, 4) 2 and a (prescribed) initial crack on the line Ω C : y = 2.0 and 1.8 ≤ x ≤ 2.2. As boundary conditions we set the displacements zero on ∂Ω. In addition, we set the regularization parameter ε = h coarse = 4.4 × 10 -3 and k = 1 2 ε = 2.2 × 10 -3 . For other choices using a different solution algorithm, we refer to [START_REF] Wheeler | An augmented-Lagrangian method for the phase-field approach for pressurized fractures[END_REF][START_REF] Heister | A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF]. We set α = 0 and the fracture toughness G c = 1.0. The mechanical parameters are Young's modulus and Poisson's ration E = 1.0 and ν s = 0.2. The injected pressure is p = 10 -3 .

Our computational goal is to measure the crack opening displacement (COD). To do so, we observe u along Ω C . Specifically, the width is determined as the jump of the normal displacements:

w = COD = [u • n]= u + n + -u -n + . (83)
Here, [•] denotes the jump. Equation (83) can be written in integral form (following, e.g., the arguments presented in [START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF], p.51):

w = COD = ∞ -∞ u • ∇ϕ dy.
The crack pattern and the corresponding mesh are displayed in Figure 2. In Figure 3, the crack opening displacement is shown for different h. The solution on the finest mesh approximates well the analytical solution provided by Sneddon.

Increasing pressure

In this second example, we study a propagating pressurized-crack. In particular, we compare the standard energy with an energy split into tensile and compressive parts as derived in [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF].

We keep the geometry and all material and model parameters as in the first example. We compute 40 time steps with time step size ∆t = 1. While applying an increasing pressure, p = ∆tp with p = 0.1, we are able to study the propagation of and the crack as shown in Figure 4. Two subiterations are used to solve the coupled elasticity phase-field problem. We study time versus crack length as shown in Figure 5. As observed in the literature, we note that the crack does not grow unless G c is exceeded and brutal crack growth occurs then. On the left, we use the energy expression and in the right subfigure, the energy-split is shown. Both expressions lead to the same results from which we infer that pressurized crack propagation as modeled in this paper does not rely on compressive energy modes.

Increasing pressure leading to joining of two cracks

The third example demonstrates a major capability of the phase-field approach: joining of two cracks without any special geometry-adapted technique. A similar test has been computed in [START_REF] Wheeler | An augmented-Lagrangian method for the phase-field approach for pressurized fractures[END_REF], where it is shown that the crack propagation is insensitive to grid perturbation. As in the previous Example 5.2, we compare the energy split into tensile and compressive parts. On the left, we use the standard energy expression and in the right subfigure, the energy-split is shown. Both expressions lead to the same results from which we infer that pressurized crack propagation as modeled in this paper does not rely on compressive energy modes.

We keep all material and model parameters as in the previous example, with the exception of having two initial cracks. The first one (Crack 1) is the same as before and centered in the middle of the domain. The second crack is vertically oriented at x = 2.6 and y ∈ [1.8, 2.2]. Consequently, the shortest distance between the two cracks is 0.4. We compute 40 time steps with time step size ∆t = 1. While applying an increasing pressure, p = ∆tp with p = 0.1, we are able to study the propagation and joining of the cracks as shown in Figure 7. The unbroken material is denoted in blue with ϕ = 1. The mushy zone is yellow/green.

Our computational goal is to detect the time when the cracks meet. In Figure 6, there are two graphs corresponding to two different mesh values, showing growth of Crack 1. We note that when Crack 1 meets the vertical crack, there is merging of the two cracks followed by no growth in Crack 1 as illustrated in Figure 7.

Joining and branching in a nonhomogeneous porous medium

The final example demonstrates another important feature of the phase-field approach: joining and branching of two cracks in a nonhomogeneous porous medium. Therefore, we adapt the Lamé coefficients such where µ varies in a range of 0.42 -1.41 and λ between 0.28 -1.28. In Figure 8 the pattern of material distribution for λ and µ is provided. The initial crack is the same as in the previous example. In addition, all other material and model parameters are taken as in the previous examples. We compute 15 time steps with time step size ∆t = 1. While applying an increasing pressure, p = ∆tp with p = 0.5, we are able to study the evolution crack patterns as shown in Figure 9. where R does not depend on g, but only on ||e(u p )|| L 4 , |k reg | and p p . Let K = B H 1 (0, R). Then K is a convex, non-empty and compact set of X 1 = L 2 . Furthermore, Φ(K) ⊂ K. It remains to prove that Φ is continuous.

Let {g n } ⊂ K be such that g n → g ∈ K in X 1 . Let (k reg (ϕ p , g n ) -k reg (ϕ p , g m ))θ δ (ϕ p -ϕ m ))ψ dx, ∀ψ ∈ X, (92)

  ) = (∇u + ∇ τ u)/2 linearized strain tensor dimensionless K density kg/m 3 B f = ρ f,0 /ρ f formation volume factor dimensionless

Figure 1 :

 1 Figure 1: Configuration and notation of the reservoir domain Ω and the crack C and zoom-in to the crack boundary ∂C where the interface law (7) is prescribed. We recall that Ω = (0, L) 3 \ C and ∂Ω = ∂(0, L) 3 ∪ ∂C.

3 ( 3 (

 33 ϕ ∆t ) + (t + ∆t)ψGe(u ∆t ) : e(u) dxdtϕ ∆t ) + (t + ∆t)ψGe(u) e(u) dxdt → T 0 (0,L) 3 (ϕ) + ψGe(u) : e(u) dxdt, as ∆t → 0. (69)

σ

  + := κtr + (e)I + 2µe D , σ -:= κtr -(e)I, with κ = 2 n µ + λ and where the deviatoric part of the strain tensor e is defined as e D := e -1 n tr(e)I, n = 2, 3. Moreover, tr + (e) = max(tr(e), 0), tr -(e) = tr(e) -tr + (e).

Figure 2 :

 2 Figure 2: Example 1 (going from top left to bottom right): Crack pattern, normal displacements(required for computing the width of the crack), mesh, and zoom-in of the mesh with the crack pattern. The crack is denoted in red with ϕ = 0 and the unbroken material is blue with ϕ = 1. In the second figure, the normal displacements are positive (red) above the fracture and negative (blue) below the fracture. The displacements are of order 10 -4 (see Figure3). The green part in the final figure shows the thickness ε of the mushy zone in the phase-field variable.

Figure 3 :

 3 Figure 3: Example 1: COD for different h. Sneddon's pink line with squares corresponds to his analytical solution.

Figure 4 :Figure 5 :

 45 Figure 4: Example 2: Crack pattern (in red with ϕ = 0) for four different time steps at T = 1, 30, 35, 40.The unbroken material is denoted in blue with ϕ = 1. The mushy zone is yellow/green.

Figure 6 :

 6 Figure 6: Example 3: Length of Crack 1 versus time for three mesh levels. The propagation of this crack stops around when the centered crack reaches the vertical-oriented crack on the right.On the left, we use the standard energy expression and in the right subfigure, the energy-split is shown. Both expressions lead to the same results from which we infer that pressurized crack propagation as modeled in this paper does not rely on compressive energy modes.

Figure 7 :

 7 Figure 7: Example 3: Crack pattern (in red with ϕ = 0) for four different time steps at T = 1, 25, 30, 35.The unbroken material is denoted in blue with ϕ = 1. The mushy zone is yellow/green.

Figure 8 :

 8 Figure 8: Example 4: µ-distribution (left) and initial crack configuration (right). The Lamé parameter µ varies in a range of 0.42 -1.41 where the red parts denotes stiff material and the blue regions smoother material.

Figure 9 :

 9 Figure 9: Example 4: Crack pattern (in red with ϕ = 0): joining, branching and non-planar crack growth in a heterogeneous medium. The unbroken material is denoted in blue with ϕ = 1. The mushy zone is yellow/green.

F 3 k 3 (

 33 n = inf{1, (g n ) 1+b + -(ϕ p ) 1+b + g n -ϕ p } (1 -α)p p div u p + ∇p p • u p + (1 -k)g n Ge(u p ) : e(u p ). (89)Then Ge(u p ) : e(u p ) ∈ L r , r > 3 andg n → g in L 2 , as n → ∞, imply (1 -k)g n Ge(u p ) : e(u p ) → (1 -k)gGe(u p ) : e(u p ) in L 2 .Next, we define h 2 = (1 -α)p p div u p + ∇p p • u p ∈ L r , r > 3, and the functionz → Q(z, y) = inf{1, z 1+b + -y 1+b + z -y }, y ∈ R, is in C 0,b [0, 1] with values between 0 and 1. Consequently, we have(0,L) 3 h 2 2 |Q(g n , ϕ p ) -Q(g, ϕ p )| 2 dx ≤ C||h 2 || 2 L 3 ||g n -g|| β L 2 ,whereβ = 2/3, for b ≥ 1/3; 2b, for 0 < b < 1/3.The convergence in L 2 of the second term in F n and the above estimate, yield thatg n → g in L 2 , as n → ∞, ⇒ F n → F in L 2 . (90)Furthermore, using thatk reg (ϕ p , g) = inf{0, -G c ε∆ϕ p -G c ε (1 -ϕ p ) + (1 -k)χGe(u p ) : e(u p )g+ inf{1, g 1+b + -(ϕ 1+b p ) + g -ϕ p } (1 -α)p p div u p + ∇p p • u p } on (0, L) 3 ,we immediately observe thatg n → g in L 2 , as n → ∞, implies k reg (ϕ p , g n ) → k reg (ϕ p , g) in L 2 . (91) Next Φ(g n ) = ϕ n satisfies G c (0,L) 3 ε∇(ϕ n -ϕ m ) • ∇ψ + 1 ε (ϕ n -ϕ m )ψ dx-(0,L) reg (ϕ p , g n )(θ δ (ϕ p -ϕ n ) -θ δ (ϕ p -ϕ m ))ψ dx = (0,L) F m -F n )ψ dx -(0,L) 3

Table 1 :

 1 Unknowns and effective coefficients

2013. The research by M. F. Wheeler was partially supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences through DOE Modeling and Simulation of coupled complex multiscale subsurface phenomena under Contract No. DE-FG02-04ER25617, MOD. 010. and T. Wick was partially supported by Conoco Phillips and is currently employed by RICAM, an institute by the Austrian Academy of Sciences.

A. Proof of Proposition 1

Proof. We will use Schauder theorem to prove Proposition 1. An alternative possibility would be to use an approach based on the pseudomonotone operators .

Let X = H 1 ((0, L) 3 ), X 1 = L 2 and let g ∈ X. We introduce g = inf{1, g + } and consider the following monotone elliptic variational problem.

Problem (84) corresponds to the minimization of a strictly convex and continuous, coercive functional Γ on X, where Γ is given by

with

The basic variational calculus implies that the minimization problem has a unique solution ϕ ∈ X.

Hence, problem (84) has a unique solution as well and the map g → ϕ = Φ(g) is well defined as a map Φ :

(1 -k)gϕGe(u p ) : e(u p ) dx and we get

and we get, after choosing ψ = ϕ n -ϕ m , that {ϕ n } is a Cauchy sequence in X = H 1 ((0, L) 3 ). Therefore, there is ϕ 0 ∈ X such that Φ(g n ) = ϕ n → ϕ 0 and ϕ 0 satisfies equation (84), with g n replaced by g. Hence ϕ 0 = Φ(g) and Φ(g n ) → Φ(g), as n → ∞. Now Schauder's theorem implies existence of a fixed point in K.

Next we observe that ϕ satisfies -ε∆ϕ + 1 ε ϕ = G ∈ L r ((0, L) 3 ), r > 3 and ∂ϕ ∂n = 0 on ∂(0, L) 3 .

(93)

By a theorem of Calderon-Zygmund we get estimate (33).