Thi-Thao-Phuong Hoang

Caroline Japhet
email: japhet@math.univ-paris13.fr

Michel Kern
email: michel.kern@inria.fr

AND Jean E Roberts
email: jean.roberts@inria.fr

SPACE-TIME DOMAIN DECOMPOSITION AND MIXED FORMULATION FOR REDUCED FRACTURE MODELS

Keywords: mixed formulations, domain decomposition, reduced fracture model, optimized Schwarz waveform relaxation, Ventcell transmission conditions, time-dependent SteklovPoincaré operator, convergence factor, nonconforming time grids AMS subject classifications. 65M55, 65M50, 65M60, 76S05, 35K20

In this paper we are interested in the "fast path" fracture and we aim to use globalin-time, nonoverlapping domain decomposition methods to model flow and transport problems in a porous medium containing such a fracture. We consider a reduced model in which the fracture is treated as an interface between the two subdomains. Two domain decomposition methods are considered: one uses the time-dependent SteklovPoincaré operator and the other uses optimized Schwarz waveform relaxation (OSWR) based on Ventcell transmission conditions. For each method, a mixed formulation of an interface problem on the space-time interface is derived, and different time grids are employed to adapt to different time scales in the subdomains and in the fracture. Demonstrations of the well-posedness of the Ventcell subdomain problems is given for the mixed formulation. An analysis for the convergence factor of the OSWR algorithm is given in the case with fractures to compute the optimized parameters. Numerical results for two-dimensional problems with strong heterogeneities are presented to illustrate the performance of the two methods.

1. Introduction. In many simulations of time-dependent physical phenomena, the domain of calculation is a union of domains with different physical properties and in which the lengths of the domains and the time scales may be very different. In particular, this is the case for a domain where there exist fractures and faults. In such a case, the fluid flows rapidly through these paths while it moves much more slowly through the rock matrix. As a result, the contaminants present in the porous medium that travel with the fluid are transported faster than in the case when there is no fracture. Thus the time scales in the fractures and in the surrounding medium are very different, and in the context of simulation, one might want to use much smaller time steps in the fractures than in the rock matrix. For simplicity we consider the case in which the domain is separated into two matrix subdomains by a fracture. The permeability in the fracture can be larger or smaller than that in the surrounding medium. A large permeability fracture corresponds to a fast pathway and a small permeability fracture corresponds to a geological barrier. Here we are interested in the "fast path" fracture. Modeling flow in porous media with fractures is challenging and requires a multi-scale approach: first, the fractures represent strong heterogeneities as they have much higher or much lower permeability than that in the surrounding medium; second, the fracture width is much smaller than any reasonable parameter of spatial discretization. Thus, to tackle the problem, one might need to refine the mesh locally around the fractures. However, this is well-known to be very computationally costly and is not useful at the macroscopic scale (i.e. when the fractures can be modeled individually). One possible approach is to treat the fractures as domains of co-dimension one, i.e. interfaces between subdomains (see [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF][START_REF] Angot | Asymptotic and numerical modelling of flows in fractured porous media[END_REF][START_REF] Bastian | Numerical simulation of multiphase flow in fractured porous media[END_REF][START_REF] Faille | A new fault model in geological basin modelling. Application of finite volume scheme and domain decomposition methods[END_REF][START_REF] Fumagalli | Numerical modelling of multiphase subsurface flow in the presence of fractures[END_REF][START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF][START_REF] Morales | The narrow fracture approximation by channeled flow[END_REF][START_REF] Morales | Interface approximation of Darcy flow in a narrow channel[END_REF][START_REF] Tunc | A model for conductive faults with non-matching grids[END_REF] and the references therein) so that one can avoid refining locally around the fractures. We point out that in these reduced fracture models, unlike in some discrete fracture models, interaction between the fractures and the surrounding porous medium is taken into account.

We are concerned with algorithms for modeling flow and transport in porous media containing such fractures. In particular, in this article we investigate two space-time domain decomposition methods, well-suited to nonmatching time grids. We use mixed finite elements [START_REF] Boffi | Mixed finite elements methods and applications[END_REF][START_REF] Roberts | Mixed and hybrid methods[END_REF] as they are mass conservative and they handle well heterogeneous and anisotropic diffusion tensors.

The first method is a global-in-time preconditionned Schur method (GTP-Schur) which uses a SteklovPoincaré-type operator. For stationary problems, this kind of method (see [START_REF] Mathew | Domain Decomposition Methods for the Numerical Solution of Partial Differential Equations[END_REF][START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF][START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF]) is known to be efficient for problems with strong heterogeneity. It uses the so-called balancing domain decomposition (BDD) preconditioner introduced and analyzed in [START_REF] Mandel | Balancing domain decomposition[END_REF][START_REF] Mandel | Balancing domain decomposition for problems with large jumps in coefficients[END_REF], and in [START_REF] Cowsar | Balancing domain decomposition for mixed finite elements[END_REF] for mixed finite elements. It involves at each iteration the solution of local problems with Dirichlet and Neumann data and a coarse grid problem to propagate information globally and to ensure the consistency of the Neumann subdomain problems. An extension to the case of unsteady problems with the construction of the time-dependent Steklov-Poincaré operator was introduced in [START_REF] Hoang | Space-time domain decomposition methods for mixed formulations of flow and transport problems in porous media[END_REF][START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF], where an interface problem on the space-time interfaces between subdomains is derived. However, for the time-dependent Neumann-Neumann problems there are no difficulties concerning consistency, and we are dealing with only a small number of subdomains, so we consider only a Neumann-Neumann type preconditioner, an extension to the nonsteady case of the method of [START_REF] Tallec | Domain decomposition methods for large linear elliptic three-dimensional problems[END_REF]. A Richardson iteration for the primal formulation was independently introduced in [START_REF] Gander | Dirichlet-neumann and neumann-neumann waveform relaxation for the wave equation[END_REF][START_REF] Kwok | Waveform Relaxation for the Time-Dependent Heat Equation[END_REF], and its convergence was analyzed. In the case of elliptic problems with fractures, a local preconditionner [START_REF] Amir | Décomposition de domaine pour un milieu poreux fracturé: un modèle en 3D avec fractures qui s'intersectent[END_REF] significantly improves the convergence of the method.

The second method is a global-in-time optimized Schwarz method (GTO-Schwarz) and uses the optimized Schwarz waveform relaxation (OSWR) approach. The OSWR and GTP-Schur methods are iterative methods that compute in the subdomains over the whole time interval, exchanging space-time boundary data through transmission conditions on the space-time interfaces. The OSWR algorithm uses more general (Robin or Ventcell) transmission operators in which coefficients can be optimized to improve convergence rates, see [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF][START_REF] Japhet | Optimized Krylov-Ventcell method. Application to convection-diffusion problems[END_REF][START_REF] Martin | An optimized Schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions[END_REF]. The optimization of the Robin (or Ventcell) parameters was analyzed in [START_REF] Bennequin | A homographic best approximation problem with application to optimized Schwarz waveform relaxation[END_REF] and the optimization method was extended to the case of discontinuous coefficients in [START_REF] Berthe | Méthodes de décomposition de domaine de type relaxation d'ondes pour l'équation de convection-diffusion instationnaire discrétisée par volumes finis[END_REF][START_REF] Berthe | Space-time domain decomposition with finite volumes for porous media applications[END_REF][START_REF] Blayo | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients. Part 1: the constant coefficients case[END_REF][START_REF] Blayo | Optimized Schwarz waveform relaxation algorithms with nonconforming time discretization for coupling convection-diffusion problems with discontinuous coefficients[END_REF][START_REF] Gander | A Schwarz waveform relaxation method for advectiondiffusion-reaction problems with discontinuous coefficients and non-matching grids[END_REF][START_REF] Hoang | Space-time domain decomposition methods for mixed formulations of flow and transport problems in porous media[END_REF][START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF]. Generalizations to heterogeneous problems with nonmatching time grids were introduced in [START_REF] Berthe | Méthodes de décomposition de domaine de type relaxation d'ondes pour l'équation de convection-diffusion instationnaire discrétisée par volumes finis[END_REF][START_REF] Berthe | Space-time domain decomposition with finite volumes for porous media applications[END_REF][START_REF] Blayo | Optimized Schwarz waveform relaxation algorithms with nonconforming time discretization for coupling convection-diffusion problems with discontinuous coefficients[END_REF][START_REF] Gander | A Schwarz waveform relaxation method for advectiondiffusion-reaction problems with discontinuous coefficients and non-matching grids[END_REF][START_REF] Haeberlein | Time space domain decomposition methods for reactive transport -Application to CO 2 geological storage[END_REF][START_REF] Halpern | Nonconforming in time domain decomposition method for porous media applications[END_REF][START_REF] Halpern | Discontinuous Galerkin and nonconforming in time optimized Schwarz waveform relaxation[END_REF][START_REF] Halpern | Optimized Schwarz waveform relaxation and discontinuous Galerkin time stepping for heterogeneous problems[END_REF][START_REF] Hoang | Space-time domain decomposition methods for mixed formulations of flow and transport problems in porous media[END_REF][START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF]. More precisely, in [START_REF] Blayo | Optimized Schwarz waveform relaxation algorithms with nonconforming time discretization for coupling convection-diffusion problems with discontinuous coefficients[END_REF][START_REF] Halpern | Discontinuous Galerkin and nonconforming in time optimized Schwarz waveform relaxation[END_REF][START_REF] Halpern | Optimized Schwarz waveform relaxation and discontinuous Galerkin time stepping for heterogeneous problems[END_REF], a discontinuous Galerkin (DG) method for the time discretization of the OSWR algorithm was introduced and analyzed for the case of nonconforming time grids. A suitable time projection between subdomains is defined using an optimal projection algorithm as in [START_REF] Gander | Algorithm PANG: Software for non-matching grid projections in 2d and 3d with linear complexity[END_REF][START_REF] Gander | A new cement to glue nonconforming grids with Robin interface conditions: the finite element case[END_REF] with no additional grid. The classical Schwarz algorithm for stationary problems with mixed finite elements was analyzed in [START_REF] Douglas | A parallel iterative procedure applicable to the approximate solution of second order partial differential equations by mixed finite element methods[END_REF]. An OSWR method with Robin transmission conditions for a mixed formulation was proposed and analyzed in [START_REF] Hoang | Space-time domain decomposition methods for mixed formulations of flow and transport problems in porous media[END_REF][START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF], where a mixed form of an interface problem on the space-time interfaces between subdomains was derived. In [START_REF] Hoang | Ventcell conditions with mixed formulations for flow in porous media[END_REF], an Optimized Schwarz method with Ventcell conditions in the context of mixed formulations was proposed. This method is not obtained in such a straightforward manner as in the case of primal formulations as Lagrange multipliers have to be introduced on the interfaces to handle tangential derivatives involved in the Ventcell conditions.

In this work, we define both a GTP-Schur and a GTO-Schwarz algorithm for a problem modeling flow of a single phase, compressible fluid in a porous medium with a fracture. A straightforward application of [START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF] would be to consider the fracture as a third subdomain and to take smaller time steps there. We consider instead however a reduced model in which the fracture is treated as an interface between two subdomains.

The definition of the GTP-Schur method is a straightforward extension of that in [START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF]. However, to define the GTO-Schwarz method, something more is needed: a linear combination between the pressure continuity equation and the fracture problem is used as a transmission condition (which leads naturally to Ventcell conditions), and a free parameter is used to accelerate the convergence rate. The well-posedness of the subdomain problems involved in the first approach was addressed in [START_REF] Boffi | Analysis of finite element approximation of evolution problems in mixed form[END_REF][START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF][START_REF] Li | Mixed methods using standard conforming finite elements[END_REF], using Galerkin's method and suitable a priori estimates. In this paper, the proof of wellposedness of both the coupled model and the Ventcell subdomain problems involved in the GTO-Schwarz approach is shown to follow from a more general theorem that covers the two cases.

Note that more general reduced models that can handle both large and small permeability fractures [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF] introduce more complicated transmission conditions on the fracture-interface (in the form of Robin type conditions, where the Robin coefficient has a physical origin), and it is not yet clear how to formulate an associated domain decomposition problem with a parameter that can be optimized. This paper is organized as follows: in the remainder of the introduction (Subsection 1.1), we state an abstract existence and uniqueness theorem for evolution problems in mixed form, the proof being deferred to Appendix A. In Section 2, we consider a reduced model with a highly permeable fracture and prove its well-posedness. Then in Section 3 we consider the GTP-Schur approach, based on physical transmission conditions, for solving the resulting problem. Different preconditionners for this method are proposed. In Section 4 we consider the GTO-Schwarz method, based on more general (e.g. Ventcell) transmission conditions, for solving the resulting problem. We prove the well-posedness of the subdomain problems with Ventcell boundary conditions. In Section 5 we consider the semi-discrete problems in time using different time grids in the subdomains. Finally, in Section 6, results of two-dimensional (2D) numerical experiments comparing the different methods are discussed.

1.1. Abstract evolution problems in mixed form. The goal of this section is to give an existence and uniqueness result for evolution problems posed in mixed form, in the spirit of the well-known theorem for weak parabolic problems (see for example [14, vol. 5]).

We consider two Hilbert spaces Σ, and M (M will be identified with its dual), and assume we have continuous bilinear forms

a : Σ × Σ -→ R, b : Σ × M -→ R, c : M × M -→ R
and a continuous linear form

L : M -→ R.
We study here an abstract version of a parabolic problem in mixed form:

Find p ∈ H 1 (0, T ; M ) and u u u ∈ L 2 (0, T ; Σ) such that, a(u u u, v v v) -b(v v v, p) = 0, ∀v v v ∈ Σ, (∂ t p, µ) M + c(p, µ) + b(u u u, µ) = L(µ), ∀µ ∈ M, p(•, 0) = p 0 , (1.1)
for some p 0 ∈ M .

We make the following hypotheses on the data:

• The bilinear form a is positive definite on Σ:

a(u u u, u u u) > 0 ∀u u u ∈ Σ, u u u = 0, ( H1
)
so that a defines a norm on Σ, and we denote by Σ a the space Σ with the norm induced by the bilinear form a. Note however that this norm will not necessarily be equivalent to the initial norm on Σ.

• The bilinear form c is positive semidefinite on M c(p, p) ≥ 0 ∀p ∈ M. ( H2
)
• The bilinear forms a and b satisfy the following compatibility condition: there exists β > 0 such that

∀u u u ∈ Σ, sup µ∈M b(u u u, µ) 2 µ 2 M + u u u 2 Σa ≥ β u u u 2 Σ . (H3)
• There exists a subspace W ⊂ M (with continuous embedding) on which the bilinear form b satisfies the stronger continuity property : there exists

C b > 0 such that b(u u u, µ) ≤ C b u u u Σa µ W , ∀u u u ∈ Σ and ∀µ ∈ W. ( H4
)
In most cases, the application of hypothesis (H3) will appear in a more natural form if it is written using the operator B : W → M associated with the bilinear form b, that is such that

∀u u u ∈ W, µ ∈ M, b(u u u, µ) = (Bu u u, µ) M .
Then, hypothesis (H3) can be written in the equivalent form: there exists β > 0 such that

∀u u u ∈ Σ, Bu u u 2 M + u u u 2 Σa ≥ β u u u 2 Σ , ( H3')
see the following remark. Remark 1.1. Hypothesis (H3) is not equivalent to the inf-sup condition. The inf-sup condition expresses the surjectivity of B T , whereas here we need the ellipticity of B with respect to the norm defined by a. This also implies a form of compatibility between a and B. This implies that a is elliptic on the kernel of B, i.e. if Bu u u = 0, a(u u u, u u u) ≥ β u u u 2 Σ . The basic existence and uniqueness result for problem (1.1) is the following: Theorem 1.2. Let M and Σ be Hilbert spaces, and let a, b and c be continuous, bilinear forms satisfying (H1) through (H3). Then, if L is a continuous linear form on M and p 0 ∈ W , where W ⊂ M satisfies (H4), then problem (1.1) has a unique solution, for which the following estimate holds:

u u u L 2 (0,T ;Σ) + p L ∞ (0,T ;M ) + ∂ t p L 2 (0,T ;M ) ≤ C L 2 L 2 (0,T ;M ) + p 0 2 W . (1.
2)

The proof of the theorem will be given in Appendix A.

Remark 1.3. The case c = 0 is allowed, and is actually the most common case (cf Theorem 2.1).

Remark 1.4. This result is a generalisation to the abstract setting of Lemma (3.1) in [START_REF] Li | Mixed methods using standard conforming finite elements[END_REF]. This problem has also been considered by Boffi and Gastaldi [START_REF] Boffi | Analysis of finite element approximation of evolution problems in mixed form[END_REF], but the estimates given there (without proof ) are different: they dispense with the regularity requirement p 0 ∈ W , at the expense of introducing weighted estimates in time to cope with the possibility of a singularity at the initial time. A proof in a more general setting is given in [START_REF] Arnold | Finite element exterior calculus for parabolic problems[END_REF]. However this proof uses semigroup theory (see Theorem 4.1 of [START_REF] Arnold | Finite element exterior calculus for parabolic problems[END_REF]), while the one we propose in this paper is with a priori estimates, in the same spirit as in [START_REF] Li | Mixed methods using standard conforming finite elements[END_REF].

We give a simple application of Theorem 1.2 (other applications will be given in Theorems 2.1 and 4.1 below).

We consider the heat equation with Dirichlet boundary conditions in mixed form. For a domain Ω ⊂ R d (d = 2 or 3) and T > 0, we look for p : Ω × [0, T ] -→ R, solution of:

∂p ∂t -∆p = f in Ω × [0, T ] p = 0 on ∂Ω × [0, T ] p(x, 0) = p 0 (x) in Ω.
(

To obtain the mixed form of (1.3), we define the spaces Σ = H(div , Ω) and M = L 2 (Ω), the bilinear forms a, b (here we will take c = 0) and the linear form

L a : Σ × Σ -→ R, a(u u u, v v v) = Ω u u u • v v v (1.4) b : Σ × M -→ R, b(u u u, µ) = Ω µ div u u u, (1.5)
L : M → R, L(µ) = Ω f µ. (1.6)
To apply Theorem 1.2, we check hypothesis (H1) to (H4) above. This is trivial for (H1) and (H2). To check (H3), we use the equivalent form (H3'). Operator B is simply the divergence, so that

Bu u u 2 M + u u u 2 Σa = Ω div u u u 2 + Ω u u u 2 = u u u 2 Σ ,
and (H3) is valid with β = 1. Last we check that (H4) is also valid with W = H 1 0 (Ω). Using Green's formula, we obtain

b(u u u, µ) = Ω µ div u u u = - Ω u u u • ∇µ,
from which (H4) follows.

2. A reduced fracture model. For a bounded domain Ω of R d (d = 2, 3) with Lipschitz boundary ∂Ω and some fixed time T > 0, we consider the problem of flow of a single phase, compressible fluid written in mixed form as follows:

s∂ t p + div u u u = q in Ω × (0, T ), u u u = -K K K∇p in Ω × (0, T ), p = 0 on ∂Ω × (0, T ), p(•, 0) = p 0 in Ω, (2.1)
where p is the pressure, u u u the velocity, q the source term, s the storage coefficient and K K K a symmetric, time independent, hydraulic, conductivity tensor (see e.g. [START_REF] Hoang | Space-time domain decomposition methods for mixed formulations of flow and transport problems in porous media[END_REF]). For simplicity we have imposed a homogeneous Dirichlet condition on the boundary. We suppose that the fracture Ω f is a subdomain of Ω , of thickness δ, that separates Ω into two connected subdomains (see Figure 2.1, left where for visualization purposes the size of δ is depicted as being relatively much larger than it is in reality),

Ω \ Ω f = Ω 1 ∪ Ω 2 , Ω 1 ∩ Ω 2 = ∅.
Also, for simplicity, we assume that Ω f consists of the intersection with Ω of a line or plane γ (depending on whether d = 2 or 3), together with the points x x x = x x x γ + sn n n where x x x γ ∈ γ, s ∈ -δ 2 , δ 2 and n n n is a unit vector normal to γ. We denote by γ i the part of the boundary of Ω i shared with the boundary of the fracture Ω f :

γ i = (∂Ω i ∩ ∂Ω f ) ∩ Ω, i = 1, 2,
and we denote by n n n i the unit, outward pointing, normal vector field on ∂Ω i . We use the convention that for any scalar, vector or tensor valued function φ defined on Ω, φ i denotes the restriction of φ to Ω i , i = 1, 2, f . We rewrite problem (2.1) as the following transmission problem:

Ω 1 Ω 2 n n n 1 n n n 2 γ 1 γ 2 Ω f δ Ω 1 Ω 2 n n n 1 n n n 2 γ
s i ∂ t p i + div u u u i = q i in Ω i × (0, T ), i = 1, 2, f, u u u i = -K K K i ∇p i in Ω i × (0, T ), i = 1, 2, f, p i = 0 on (∂Ω i ∩ ∂Ω) × (0, T ), i = 1, 2, f, p i = p f on γ i × (0, T ), i = 1, 2, u u u i • n n n i = u u u f • n n n i on γ i × (0, T ), i = 1, 2, p i (•, 0) = p 0,i in Ω i , i = 1, 2, f. (2.2)
In the reduced fracture model, the fracture Ω f is treated as a simple interface γ between subdomains Ω 1 and Ω 2 (see Figure 2.1, right). We use the notation ∇ τ (respectively div τ ) for the tangential gradient (respectively tangential divergence) operators along the fracture γ. We denote by s γ and K K K γ the storage coefficient and the permeability tensor in the (d -1)-dimensional fracture γ. The reduced model that we consider was derived in [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF][START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF]. It may be obtained by averaging across the transversal cross sections of the d-dimensional fracture Ω f . It consists of equations in the subdomains,

s i ∂ t p i + div u u u i = q i in Ω i × (0, T ), u u u i = -K K K i ∇p i in Ω i × (0, T ), p i = 0 on (∂Ω i ∩ ∂Ω) × (0, T ), p i = p γ on γ × (0, T ), p i (•, 0) = p 0,i in Ω i , for i = 1, 2, (2.3)
and equations in the interface fracture

s γ ∂ t p γ + div τ u u u γ = u u u 1 • n n n 1|γ + u u u 2 • n n n 2|γ in γ × (0, T ), u u u γ = -K K K γ δ∇ τ p γ in γ × (0, T ), p γ = 0 on ∂γ × (0, T ), p γ (•, 0) = p 0,γ in γ.
(2.4)

These equations are the mass conservation equation and the Darcy equation in the subdomain together with the lower dimensional mass conservation and Darcy equations in the fracture of co-dimension 1. These two systems are coupled: the fracture sees the subdomain through the source term in the conservation equation in the fracture which represents the difference between the fluid entering the fracture from one subdomain and that exiting through the other subdomain. Each subdomain sees the fracture through the Dirichlet boundary condition imposed on the part of its boundary common with the fracture. We make the hypothesis of the following compatibility conditions: p 0,i = p 0,γ on γ, for i = 1, 2. For a general mathematical treatment of this type of problem in the stationary case see [START_REF] Hung | Phénomènes de transmission à travers des couches minces de conductivité élevée[END_REF].

To prove the well-posedness of problem (2.3)-(2.4), we shall use the abstract framework of Subsection 1.1 and apply Theorem 1.2. We first write the weak formulation for problem (2.3)-(2.4), and define the appropriate function spaces, and the forms on these spaces. We use the convention that if V is a space of functions, then V V V is a space of vector functions having each component in V . For an arbitrary domain O, we denote by (•,

•) O the inner product in L 2 (O) or L 2 L 2 L 2 (O)
and by and

• O the L 2 (O)-norm or L 2 L 2 L 2 (O)-norm.
To write the weak formulation of (2.3)-(2.4), we define the following Hilbert spaces:

M = µ = (µ 1 , µ 2 , µ γ ) ∈ L 2 (Ω 1 ) × L 2 (Ω 2 ) × L 2 (γ) , Σ = v v v = (v v v 1 , v v v 2 , v v v γ ) ∈ L 2 L 2 L 2 (Ω 1 ) × L 2 L 2 L 2 (Ω 2 ) × L 2 L 2 L 2 (γ) : div v v v i ∈ L 2 (Ω i ), i = 1, 2,
and

div τ v v v γ - 2 i=1 v v v i • n n n i|γ ∈ L 2 (γ) ,
equipped with the norms

µ 2 M = 2 i=1 µ i 2 Ωi + µ γ 2 γ , v v v 2 Σ = 2 i=1 v v v i 2 Ωi + div v v v i 2 Ωi + v v v γ 2 γ + div τ v v v γ - 2 i=1 v v v i • n n n i|γ 2 γ .
We define the following bilinear forms

a : Σ × Σ -→ R (u u u, v v v) → a(u u u, v v v) = 2 i=1 K K K -1 i u u u i , v v v i Ωi + (K K K γ δ) -1 u u u γ , v v v γ γ , b : Σ × M -→ R (u u u, µ) → b(u u u, µ) = 2 i=1 (div u u u i , µ i ) Ωi + div τ u u u γ - 2 i=1 u u u i • n n n i|γ , µ γ γ , c s : M × M -→ R (η, µ) → c s (η, µ) = 2 i=1 (s i η i , µ i ) Ωi + (s γ η γ , µ γ ) γ ,
and the linear form

L q : M -→ R µ → L q (µ) = 2 i=1 (q i , µ i ) Ωi .
With these spaces and forms, the weak form of (2.3)-(2.4) can be written as follows:

Find p ∈ H 1 (0, T ; M ) and u u u ∈ L 2 (0, T ; Σ) such that, a(u u u, v v v) -b(v v v, p) = 0, ∀v v v ∈ Σ, c s (∂ t p, µ) + b(u u u, µ) = L q (µ), ∀µ ∈ M, (2.5)
together with the initial conditions

p i (•, 0) = p 0,i in Ω i , i = 1, 2, p γ (•, 0) = p 0,γ in γ, (2.6)
for p 0,i ∈ L 2 (Ω i ), i = 1, 2 and p 0,γ ∈ L 2 (γ). We also define the space

H 1 * := {µ = (µ 1 , µ 2 , µ γ ) ∈ H 1 (Ω 1 ) × H 1 (Ω 2 ) × H 1 0 (γ) : µ i = 0 on ∂Ω i ∩ ∂Ω, and µ i = µ γ on γ, i = 1, 2}, equipped with the norm µ 2 H 1 * = µ 2 M + 2 i=1 ∇µ i 2 Ωi + ∇ τ µ γ 2 γ .
The well-posedness of problem (2.5)-(2.6) is given by the following theorem: Theorem 2.1. Assume that there exist four positive constants s -and s + , K - and

K + such that • s -≤ s i (x) ≤ s + for a.e. x ∈ Ω i , i = 1, 2, • s -≤ s γ (x) ≤ s + for a.e. x ∈ γ, • ς T K K K -1 i (x)ς ≥ K -|ς| 2 , and |K K K i (x)ς| ≤ K + |ς|, for a.e. x ∈ Ω i , ∀ς ∈ R d , i = 1, 2, • η T (K K K γ (x)δ) -1 η ≥ K -|η| 2 and |(K K K γ (x)δ) -1 η| ≤ K + |η| for a.e. x ∈ γ, ∀η ∈ R d-1 . If q is in L 2 (0, T ; M ) and p 0 = (p 0,1 , p 0,2 , p 0,γ ) in H 1 * then problem (2.5)-(2.6) has a unique solution (p, u u u) ∈ H 1 (0, T ; M ) × L 2 (0, T ; Σ).
Proof. First notice that under the assumptions on s i and s γ stated in the theorem, c s defines an inner product on M × M , and that the associated norm is equivalent to the original norm on M . We will apply Theorem 1.2, in the case c = 0. The bilinear forms a and b are obviously continuous, and, with the hypotheses concerning K K K i and K K K γ , a is positive definite on Σ.

We now check hypothesis (H3). This is easiest to do using the operator form (H3'), where for this problem B is defined by

∀u u u ∈ Σ, Bu u u = div u u u 1 , div u u u 2 , div τ u u u γ - 2 i=1 u u u i • n n n i|γ .
The result follows from the hypothesis on K i and K γ and the definition of the norm on Σ.

Last, to check hypothesis (H4), we can take W = H 1 * , and use Green's formula to see that for

u u u ∈ Σ and µ ∈ W b(u u u, µ) = 2 i=1 (div u u u i , µ i ) Ωi + div τ u u u γ - 2 i=1 u u u i • n n n i|γ , µ γ γ = 2 i=1 -(u u u i , ∇µ i ) Ωi + u u u i • n n n i|γ , µ i γ -(u u u γ , ∇ τ µ γ ) γ - 2 i=1 u u u i • n n n i|γ , µ γ γ = - 2 i=1 (u u u i , ∇µ i ) Ωi -(u u u γ , ∇ τ µ γ ) γ because µ i|γ = µ γ .
To conclude, we bound the terms of the right hand side:

(u u u i , ∇µ i ) Ωi ≤ K + u u u i Ωi µ i H 1 (Ωi) , (u u u γ , ∇ τ µ γ ) γ ≤ K + u u u γ γ µ γ H 1 (γ) ,
from which hypothesis (H4) easily follows.

It is natural to use domain decomposition methods for obtaining a numerical solution of problem (2.2) or problem (2.3)-(2.4), especially as these methods make it possible to take different time steps in the subdomains and in the fracture. For problem (2.2), it would be a straightforward application of the methods introduced in [START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF] while for problem (2.3)-(2.4), we need to derive a different formulation. In the following, we present two global-in-time domain decomposition methods for solving (2.3)-(2.4) based on different transmission conditions. A space-time interface problem, which will be solved iteratively, is derived for each approach.

3. Global-in-time preconditioned Schur (GTP-Schur): using the timedependent Steklov-Poincaré operator. The Global-in-time preconditioned Schur (GTP-Schur) method is directly derived from the formulation of problem (2.3) -(2.4). To obtain the interface problem for this method, we need to introduce some notation. For a bounded domain O ∈ R d (d = 2, 3) with Lipschitz boundary ∂O containing an open subset γ ⊂ ∂O, we define the space

H 1 * ,γ (O) := µ ∈ H 1 (O) : µ = 0 on (∂O \ γ) .
Then we define the following Dirichlet to Neumann operators S DtN i , i = 1, 2 :

S DtN i : H 1 (0, T ; H 1 2 00 (γ)) × L 2 (0, T ; L 2 (Ω i )) × H 1 * ,γ (Ω i ) → L 2 0, T ; (H 1 2 00 (γ)) S DtN i (λ, q, p 0 ) → u u u i • n n n i|γ ,
where (p i , u u u i ), i = 1, 2, is the solution of the problem

s i ∂ t p i + div u u u i = q in Ω i × (0, T ), u u u i = -K K K i ∇p i in Ω i × (0, T ), p i = 0 on (∂Ω i ∩ ∂Ω) × (0, T ), p i = λ on γ × (0, T ), p i (•, 0) = p 0 in Ω i . (3.1)
Remark 3.1. A straightforward application of Theorem 1.2 shows the wellposedness of subdomain problem (3.1). See also [START_REF] Hoang | Space-time domain decomposition methods for mixed formulations of flow and transport problems in porous media[END_REF][START_REF] Li | Mixed methods using standard conforming finite elements[END_REF] for a direct proof.

Problem (2.4) is reduced to an interface problem with unknowns λ and u γ :

s γ ∂ t λ + div τ u u u γ = 2 i=1 S DtN i (λ, q i , p 0,i ) in γ × (0, T ), u u u γ = -K K K γ δ∇ τ λ in γ × (0, T ), λ = 0 on ∂γ × (0, T ), λ(•, 0) = p 0,γ in γ. (3.2)
or equivalently

s γ ∂ t λ + div τ u u u γ - 2 i=1 S DtN i (λ, 0, 0) = 2 i=1 S DtN i (0, q i , p 0,i ) in γ × (0, T ), u u u γ = -K K K γ δ∇ τ λ in γ × (0, T ), λ = 0 on ∂γ × (0, T ), λ(•, 0) = p 0,γ in γ, (3.3
) or in compact form (space-time),

S λ u u u γ = χ.
This problem is solved using an iterative solver such as GMRES since due to the time derivative the system is nonsymmetric.

To improve the convergence of the iterative algorithm, we will consider two preconditioners. The first, introduced in [START_REF] Amir | Décomposition de domaine pour un milieu poreux fracturé: un modèle en 3D avec fractures qui s'intersectent[END_REF], arises from the observation that the interface problem is dominated by the second order operator (div τ (K K K γ δ∇ τ )) since the Steklov-Poincaré operator is of lower order (first order). This is even more the case when the permeability in the fracture is much larger than that in the surrounding domain. Thus one choice for a preconditioner is P P P -1 loc defined by taking the discrete counterpart of the operator (div τ (K K K γ δ∇ τ ))

-1 . We have

P P P -1 loc : L 2 (γ) → L 2 (γ) g γ → pγ , where (p γ , ũ u u γ ) is the solution of the problem div τ ũ u u γ = g γ in γ, ũ u u γ = -K K K γ δ∇ τ pγ in γ, pγ = 0 on ∂γ.
This preconditioner was introduced for elliptic problems, and it was shown numerically [START_REF] Amir | Décomposition de domaine pour un milieu poreux fracturé: un modèle en 3D avec fractures qui s'intersectent[END_REF] that it significantly improves the convergence of the algorithm, especially, as mentionned before, for high permeability in the fracture.

A second possibility is to use the Neumann-Neumann preconditioner as was done in [START_REF] Hoang | Space-time domain decomposition methods for mixed formulations of flow and transport problems in porous media[END_REF][START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF] for ordinary domain decomposition algorithms (i.e. without fractures). The preconditioned problem is then

P P P -1
NN ϕ = χ, with

P P P -1 NN := σ 1 ( ŠDtN 1 ) -1 + σ 2 ( ŠDtN 2 ) -1 ,
where

σ i : Γ × (0, T ) → [0, 1] is such that σ 1 + σ 2 = 1. If K K K i = K i I I I and K i is constant in each subdomain then σ i = K i K 1 + K 2 .
The operator ( ŠDtN

( ŠDtN i ) -1 : L 2 0, T ; L 2 (γ) → H 1 0, T ; L 2 (γ) ŠDtN i -1 (ϕ) → p i|γ ,
where (p i , u u u i ), i = 1, 2, is the solution of the problem

s i ∂ t p i + div u u u i = 0 in Ω i × (0, T ), u u u i = -K K K i ∇p i in Ω i × (0, T ), p i = 0 on (∂Ω i ∩ ∂Ω) × (0, T ), -u u u i • n n n i = ϕ on γ × (0, T ), p i (•, 0) = 0 in Ω i . (3.4)
In Section 6, we will carry out numerical experiments and compare the performance of these two preconditioners.

4. Global-in-time optimized Schwarz (GTO-Schwarz): using optimized Schwarz waveform relaxation. While the extension of the GTP-Schur method to handle the fracture model is straightforward, the extension of the GTO-Schwarz method to the fracture problem needs something more. Indeed, instead of imposing Dirichlet boundary conditions on γ × (0, T ) when solving the fracture problem as was done for the GTP-Schur method, for the GTO-Schwarz approach one uses optimized Robin transmission conditions. Thus, we introduce new transmission conditions, that combine the equation for continuity of the pressure across the fracture with the flow equations (2.4) in the fracture. These new transmission conditions contain a free parameter, which is used to accelerate the convergence. This is an extension of the OSWR method with optimized Robin parameters studied in [START_REF] Hoang | Space-time domain decomposition methods for mixed formulations of flow and transport problems in porous media[END_REF][START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF] in which Robin-to-Robin transmission conditions are considered in mixed form. Here however, because of the fracture problem, we obtain what we will call Ventcell-to-Robin transmission conditions as described below.

4.1. Ventcell-to-Robin transmission conditions. The new transmission conditions are derived by introducing Lagrange multipliers p i,γ , i = 1, 2, with p i,γ representing the trace on the interface γ of the pressure p i in the subdomain Ω i . As the pressure is continuous across the interface, one has

p 1,γ = p 2,γ = p γ , on γ × (0, T ). (4.1)
We then rewrite the Darcy equation in the fracture associated with each p i,γ as

u u u γ,i = -K K K γ δ∇ τ p i,γ , on γ × (0, T ), i = 1, 2.
We have used the notation u u u γ,i , i = 1, 2, instead of u u u i,γ to insist on the fact that u u u γ,i is not the tangential component of a trace of u u u i on γ. In fact, u u u γ,i , i = 1, 2, represents the tangential velocity in the fracture associated with the pressure p i so that

u u u γ,1 = u u u γ,2 = u u u γ , on γ × (0, T ), i = 1, 2.
With the notation introduced above, the flow equation (2.4) in the fracture can be rewritten, for i = 1, 2, and j = (3 -i), as

-u u u i • n n n i + s γ ∂ t p i,γ + div τ u u u γ,i = -u u u j • n n n i , on γ × (0, T ), u u u γ,i = -K K K γ δ∇ τ p i,γ , on γ × (0, T ), p i,γ = 0 on ∂γ × (0, T ), p i,γ (•, 0) = p 0,γ in γ. (4.2)
In the context of domain decomposition, (4.1) and (4.2) are the coupling conditions between the subdomains. As in the case without a fracture we take a linear combination of these conditions (for a parameter α > 0), but here we obtain equivalent Ventcell-to-Robin transmission conditions (instead of Robin-to-Robin):

-u u u 1 • n n n 1 + αp 1,γ + s γ ∂ t p 1,γ + div τ u u u γ,1 = -u u u 2 • n n n 1 + αp 2,γ u u u γ,1 = -K K K γ δ∇ τ p 1,γ on γ × (0, T ), (4.3) -u u u 2 • n n n 2 + αp 2,γ + s γ ∂ t p 2,γ + div τ u u u γ,2 = -u u u 1 • n n n 2 + αp 1,γ u u u γ,2 = -K K K γ δ∇ τ p 2,γ on γ × (0, T ), (4.4
) Using these transmission conditions together with the boundary and initial conditions

p 1,γ = p 2,γ = 0 on ∂γ × (0, T ), p 1,γ (•, 0) = p 2,γ (•, 0) = p 0,γ in γ, (4.5)
the subdomain problem is obtained by imposing Ventcell boundary conditions on γ × (0, T ), i = 1, 2, j = 3 -i:

s i ∂ t p i + div u u u i = q in Ω i × (0, T ), u u u i = -K K K i ∇p i in Ω i × (0, T ), -u u u i • n n n i + αp i,γ + s γ ∂ t p i,γ + div τ u u u γ,i = -u u u j • n n n i + αp j,γ on γ × (0, T ), u u u γ,i = -K K K γ δ∇ τ p i,γ in γ × (0, T ), p i = 0 on (∂Ω i ∩ ∂Ω) × (0, T ), p i,γ = 0 on ∂γ × (0, T ), p i (•, 0) = p 0 in Ω i , p i,γ (•, 0) = p 0,γ in γ, (4.6
) where the quantity on the right hand side of the third equation will be known in the context of an iterative method for solving (2.3)-(2.4). In the next subsection we prove that problem (4.6) is well-posed.

s O ∂ t p O + div u u u O = q in O × (0, T ), u u u O = -K K K O ∇p O in O × (0, T ), -u u u O • n n n + αp γ + s γ ∂ t p γ + div τ u u u γ = θ γ on ∂γ × (0, T ), u u u γ = -K K K γ δ∇ τ p γ in γ × (0, T ), p O = 0 on (∂O \ γ) × (0, T ), p γ = 0 on ∂γ × (0, T ), p O (•, 0) = p 0,O in Ω i , p γ (•, 0) = p 0,γ in γ, (4.7)
where θ γ is a function defined on γ × (0, T ), and α ∈ R, α > 0. In order to write the weak formulation of (4.7), we need to define the following Hilbert spaces:

M O = µ = (µ O , µ γ ) ∈ L 2 (O) × L 2 (γ) , Σ O = v v v = (v v v O , v v v γ ) ∈ L 2 L 2 L 2 (O) × L 2 L 2 L 2 (γ) : div v v v O ∈ H(div , O) and div τ v v v γ -v v v O • n n n |γ ∈ L 2 (γ) ,
equipped with the norms

µ 2 M O = µ O 2 O + µ γ 2 γ , v v v 2 Σ O = v v v O O + div v v v O 2 O + v v v γ 2 γ + div τ v v v γ -v v v O • n n n |γ 2 γ .
Then define the bilinear forms

a O : Σ O × Σ O -→ R, a O (u u u, v v v) = K K K -1 O u u u O , v v v O O + (K K K γ δ) -1 u u u γ , v v v γ γ b O : Σ O × M O -→ R, b O (u u u, µ) = (div u u u O , µ O ) O + div τ u u u γ -u u u O • n n n |γ , µ γ γ , c O : M O × M O -→ R, c O (η, µ) = (αη γ , µ γ ) γ , c s,O : M O × M O -→ R, c s,O (η, µ) = (s O η O , µ O ) O + (s γ η γ , µ γ ) γ ,
and the linear form

L q,O : M O -→ R, L q,O (µ) = (q, µ O ) O + (θ γ , µ γ ) γ .
With these spaces and forms, the weak form of (4.7) can be written as follows:

For a.e. t ∈ (0, T ), find

p(t) ∈ M O and u u u(t) ∈ Σ O such that a O (u u u, v v v) -b O (v v v, p) = 0 ∀v v v ∈ Σ O , c s,O (∂ t p, µ) + c O (p, µ) + b O (u u u, µ) = L q,O (µ) ∀µ ∈ M O , (4.8)
together with the initial conditions

p O (•, 0) = p 0,O in O, p γ (•, 0) = p 0,γ in γ. (4.9)
We will also make use of

H 1,1 * (O, γ) := H 1 * ,γ (O) × H 1 0 (γ).
The well-posedness of problem (4.8)-(4.9) is given by the following theorem: Theorem 4.1. Assume that there exist positive constants s -, s + , K -, K + with

• s -≤ s O (x) ≤ s + for a.e. x ∈ O, • s -≤ s γ (x) ≤ s + for a.e. x ∈ γ, • ς T K K K -1 O (x)ς ≥ K -|ς| 2 , and |K K K i (x)ς| ≤ K + |ς|, for a.e. x ∈ O, ∀ς ∈ R d , • η T K K K -1 γ (x)δη ≥ K -|η| 2 and |(K K K γ (x)δ) -1 η| ≤ K + for a.e. x ∈ γ, ∀η ∈ R d-1 . If q is in L 2 (0, T ; M O ), p 0 in H 1,1
* (O, γ) and θ γ in L 2 (0, T ; L 2 (γ)) then problem (4.8)-(4.9) has a unique solution (p, u u u) ∈ H 1 (0, T ; M O ) × L 2 (0, T ; Σ O ).

Proof. As in the case with Dirichlet boundary conditions, we apply Theorem 1.2. Again notice that, under the hypotheses on s O and s γ , c s,O defines an inner product on M , equivalent to its usual inner product.

• It is clear that a, b and c are all continuous forms, that a is positive definite (due to the hypotheses on K O and K γ ) and that c O is positive (since α ≥ 0), so that hypotheses (H1) and (H2) in Theorem 1.2 hold.

• To verify hypothesis (H3), we define the operator B by

Bu u u = div u u u O , div τ u u u γ -u u u O • n n n |γ ,
and it follows easily that

Bu u u 2 M + a O (u u u, u u u) = div u u u O 2 O + div τ u u u γ -u u u O • n n n |γ 2 γ + K K K -1 O u u u O , u u u O O + (K K K γ δ) -1 u u u γ , u u u γ γ ≥ β u u u 2 Σ ,
for some β > 0, again because of the lower bounds on K O and K γ .

• Last, to check (H4), we proceed as in Theorem 2.1, and use Green's formula for

u u u ∈ Σ and µ ∈ H 1,1 * (O, γ). b O (u u u, µ) = (div u u u O , µ O ) O + div τ u u u γ -u u u O • n n n |γ , µ γ γ = -(u u u O , ∇µ O ) O + u u u O • n n n |γ , µ O γ -(u u u γ , ∇ τ µ γ ) γ -u u u O • n n n |γ , µ γ γ ,
from which the proof of the theorem follows.

The interface problem.

As for the GTP-Schur method, we derive an interface problem which in this case is associated with Ventcell-to-Robin transmission conditions (4.3)-(4.4). Towards this end, we define the following Ventcell-to-Robin operator S VtR i , which depends on the parameter α, for i = 1, 2; j = (3 -i):

S VtR i : L 2 (0, T ; L 2 (γ)) × L 2 (0, T ; L 2 (Ω i )) × H 1 * (Ω i ) × H 1 0 (γ) → L 2 (0, T ; L 2 (γ)) S VtR i (θ γ , q, p 0 , p 0,γ ) → -u u u i • n n n j|γ + αp i,γ ,
where (p i , u u u i , p i,γ , u u u γ,i ) is the solution of the subdomain problem with Ventcell boundary conditions

s i ∂ t p i + div u u u i = q in Ω i × (0, T ), u u u i = -K K K i ∇p i in Ω i × (0, T ), -u u u i • n n n i + αp i,γ + s γ ∂ t p i,γ + div τ u u u γ,i = θ γ on γ × (0, T ), u u u γ,i = -K K K γ δ∇ τ p i,γ in γ × (0, T ), p i = 0 on (∂Ω i ∩ ∂Ω) × (0, T ), p i,γ = 0 on ∂γ × (0, T ), p i (•, 0) = p 0 in Ω i , p i,γ (•, 0) = p 0,γ in γ. (4.10)
The interface problem with two Lagrange multipliers is then

θ γ,1 = S VtR 2 (θ γ,2 , q 2 , p 0,2 , p 0,γ ) θ γ,2 = S VtR 1 (θ γ,1
, q 1 , p 0,1 , p 0,γ ) on γ × (0, T ), (4.11) or equivalently

θ γ,1 -S VtR 2 (θ γ,2 , 0, 0, 0) = S VtR 2 (0, q 2 , p 0,2 , p 0,γ ) θ γ,2 -S VtR 1 (θ γ,1 , 0, 0, 0) = S VtR 1 (0, q 1 , p 0,1 , p 0,γ ) on γ × (0, T ). (4.12)
The discrete counterpart of this problem can be solved iteratively using Jacobi iterations or GMRES. The former choice yields an algorithm equivalent to the OSWR algorithm for the reduced fracture model (2.3) -(2.4) and is written as follows: starting with an initial guess θ 0 γ,j , j = 3 -i, on γ × (0, T ) for the first iteration,

-u u u 0 i • n n n i + α p 0 i,γ + s γ ∂ t p 0 i,γ + div τ u u u 0 γ,i = θ 0 γ,j ,
then at the k th iteration, k = 1, . . . , solve in each subdomain the time-dependent problem, for i = 1, 2; j = (3 -i),

s i ∂ t p k i + div u u u k i = q i in Ω i × (0, T ), u u u k i = -K K K i ∇p k i in Ω i × (0, T ), -u u u k i • n n n i + αp k i,γ + s γ ∂ t p k i,γ + div τ u u u k γ,i = θ k-1 γ,j on γ × (0, T ), u u u k γ,i = -K K K f,τ δ∇ τ p k i,γ on γ × (0, T ), p k i = 0 on (∂Ω i ∩ ∂Ω) × (0, T ), p k i,γ = 0 on ∂γ × (0, T ), p k i (•, 0) = p 0,i in Ω i , p k i,γ (•, 0) = p 0,γ in γ, (4.13) with θ k-1 γ,j = -u u u k-1 j
• n n n i + αp k-1 j,γ on γ × (0, T ). The convergence of algorithm (4.13) depends on the choice of the parameter α. Thus we extend the analysis for the convergence factor of the OSWR algorithm derived in the case without fractures [START_REF] Bennequin | A homographic best approximation problem with application to optimized Schwarz waveform relaxation[END_REF][START_REF] Gander | Optimized Schwarz methods[END_REF][START_REF] Japhet | Optimized Krylov-Ventcell method. Application to convection-diffusion problems[END_REF] to this algorithm and from that, one can calculate the optimal or optimized values of the parameter α.

4.4.

Convergence factor formula for computing the optimized parameter. In this section, we extend the two domain analysis [START_REF] Bennequin | A homographic best approximation problem with application to optimized Schwarz waveform relaxation[END_REF][START_REF] Halpern | Nonconforming in time domain decomposition method for porous media applications[END_REF][START_REF] Hoang | Space-time domain decomposition methods for mixed formulations of flow and transport problems in porous media[END_REF][START_REF] Japhet | Optimized Krylov-Ventcell method. Application to convection-diffusion problems[END_REF][START_REF] Martin | An optimized Schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions[END_REF] to derive the convergence factor of the OSWR algorithm introduced in Section 4.3 for a reduced fracture model for compressible flow. Towards this end, we consider the two halfspace decomposition Ω -= R -× R, Ω + = R + × R and write the OSWR algorithm, applied to the fractured model, in the primal formulation: at the k th Jacobi iteration, solve

s -∂ t p k -+ div (-K K K -∇p k -) = q in Ω -× (0, T ), K K K - ∂p k - ∂n n n - + αp k -+ s γ ∂ t p k -+ div τ (-K K K f,τ δ∇ τ p k -) = K K K + ∂p k-1 + ∂n n n - +αp k-1 + on γ × (0, T ), p k -(•, 0) = p 0 in Ω -, (4.14)
and

s + ∂ t p k + + div (-K K K -∇p k + ) = q in Ω + × (0, T ), K K K + ∂p k + ∂n n n + + αp k + + s γ ∂ t p k + + div τ (-K K K f,τ δ∇ τ p k + ) = K K K - ∂p k-1 - ∂n n n + +αp k-1 - on γ × (0, T ), p k + (•, 0) = p 0 in Ω + , (4.15
) where γ = {x = 0} is the fracture. We assume that the permeability is isotropic:

K K K pm = K ± I I I, and K K K f,τ = K f ,
where I I I is the 2D identity matrix, and that the solution of the problem decays at infinity. As the problem is linear, we only consider q = 0 and p 0 = 0, and analyse the convergence of (4.14)-(4.15) to the zero solution. We use a Fourier transform in time and in the y direction with parameters ω and η, respectively, to obtain the Fourier functions pk ± in time t and y of p k ± , as the solutions to the ordinary differential equation in x

-K ∂ 2 p ∂x 2 + siω + Kη 2 p = 0. Thus p = A(η, ω)e r + x + B(η, ω)e r -x ,
where r ± are the roots of the characteristic equation

-Kr 2 + siω + Kη 2 = 0, so r ± = ± √ ∆ 2K , ∆ = 4K siω + Kη 2 .
Here and throughout this article, we use the square root symbol √ to denote the complex square root with positive real part. In order to work with at least square integrable functions in time and space, we look for solutions which do not increase exponentially in x. Since r + > 0 and r -< 0, we obtain pk

- = A k (η, ω)e r + (s-,K-,η,ω)x , pk + = B k (η, ω)e r -(s+,K+,η,ω)x .
Substituting these formulas into the transmission conditions on the interface γ ×(0, T ) (i.e. the second equations of (4.14) and (4.15)), we find

K -r + (s -, K -, η, ω) + α + s γ iω + K f δη 2 pk -(0, η, ω) = (K + r -(s + , K + , η, ω) + α) pk-1 + (0, η, ω), -K + r -(s + , K + , η, ω) + α + s γ iω + K f δη 2 pk + (0, η, ω) = (-K -r + (s -, K -, η, ω) + α) pk-1 -(0, η, ω). (4.16)
From (4.16) using induction ans using ζ to denote s γ iω + K f δη 2 , we obtain

p2k -(0, η, ω) = K + r -(s + , K + , η, ω) + α K -r + (s -, K -, η, ω) + α + ζ p2k-1 + (0, η, ω) = K + r -(s + , K + , η, ω) + α K -r + (s -, K -, η, ω) + α + ζ -K -r + (s -, K -, η, ω) + α -K + r -(s + , K + , η, ω) + α + ζ p2k-2 - (0, η, ω) = ρ k f p0 -(0, η, ω),
where

ρ f = K + r -(s + , K + , η, ω) + α K -r + (s -, K -, η, ω) + α + ζ -K -r + (s -, K -, η, ω) + α -K + r -(s + , K + , η, ω) + α + ζ) ,
is the convergence factor of the algorithm (4.14)-(4.15). Similarly, we obtain

p2k + (0, η, ω) = ρ k f p0 + (0, η, ω),
Thus, we can calculate the parameter α in such a way as to minimize this continuous convergence factor:

min α>0 max |η|∈[ π L , π h ],|ω|∈[ π T , π ∆t ] ρ f (s + , K + , s -, K -, α, η, ω) , (4.17)
where L is the length of the fracture, h is the spatial mesh size, T is the final time and ∆t is the maximum time step of the discretization in time. Remark 4.2. One could make use of the two-sided Robin as in [START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF]. In this article, the optimized one-sided Robin parameter works well since in the test case we considered, the two subdomains Ω 1 and Ω 2 (representing the rock matrix) have similar physical properties (though a comparison of the performance of the one-sided and twosided Robin might be considered). In our applications, the fracture is assumed to have much larger permeability than the surrounding domain, which suggests that the time step inside the fracture should be small compared with that of the surrounding matrix subdomains. As both of the methods derived in Sections 3 and 4 are global in time, i.e. the subdomain problem is solved over the whole time interval before the information is exchanged on the space-time interface, we can use different time steps in the fracture and in the rock matrix. In the next section, we consider the semi-discrete problem in time with nonconforming time grids.

Nonconforming discretization in time.

Let T 1 , T 2 and T γ be three different partitions of the time interval (0, T ) into subintervals J i m = (t i m-1 , t i m ] for m = 1, . . . , M i , and i = 1, 2, γ, (see Figure 5.1). For simplicity, we consider uniform partitions only, and denote by ∆t i , i = 1, 2, γ the corresponding time steps. Assume that ∆t γ ∆t i , i = 1, 2. We use the lowest order discontinuous Galerkin method [START_REF] Blayo | Optimized Schwarz waveform relaxation algorithms with nonconforming time discretization for coupling convection-diffusion problems with discontinuous coefficients[END_REF][START_REF] Halpern | Discontinuous Galerkin and nonconforming in time optimized Schwarz waveform relaxation[END_REF][START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF], which is a modified backward Euler method. The same idea can be generalized to higher order methods.

We denote by P 0 (T i , L 2 (γ)) the space of functions piecewise constant in time on grid T i with values in L 2 (γ):

P 0 (T i , L 2 (γ)) = ψ : (0, T ) → L 2 (γ), ψ is constant on J, ∀J ∈ T i .
In order to exchange data on the space-time interface between different time grids, we use, for i, j in {1, 2, γ}, the L 2 projection Π ji from P 0 (T i , L 2 (γ)) to P 0 (T j , L 2 (γ)): for ψ ∈ P 0 (T i , L 2 (γ)), Π ji ψ | J j m is the average value of ψ on J j m , for m = 1, . . . , M j .

t x Ω 1 Ω 2 γ T = M 1 ∆t 1 = M 2 ∆t 2 = M γ ∆t γ ∆t 1 ∆t 2 Fig. 5.1.
Nonconforming time grids in the rock matrix and in the fracture.

For the GTP-Schur method. The unknown λ on the interface represents the fracture pressure, thus λ is piecewise constant in time on grid T γ . In order to obtain Dirichlet boundary data for solving subdomain problem (3.1), we project λ into P 0 (T i , L 2 (γ)), for i = 1, 2 :

p i = Π iγ (λ), on γ, i = 1, 2.
The semi-discrete counterpart of the interface problem (3.2) is obtained by weakly enforcing the fracture problem over each time sub-interval of T γ as follows

s γ λ m+1 -λ m + t m+1 γ t m γ div τ u u u m+1 γ = t m+1 γ t m γ 2 i=1 Π γi S DtN i (Π iγ (λ), q i , p 0,i ) , u u u m+1 γ = -K K K γ δ∇ τ λ m+1 in γ, λ m+1 = 0 on ∂γ, λ 0 = p 0,γ in γ, (5.1)
where λ m = λ J γ m , for m = 0, . . . , M γ -1. For a function piecewise constant in time ϕ on the fine grid T γ , the semi-discrete Neumann-Neumann preconditioner (still denoted by P P P -1 N N ) is defined by:

P P P -1 N N ϕ := 2 i=1 σ i Π γi ŠDtN i -1 (Π iγ (ϕ)) , (5.2)
where we have solved the subdomain problem with Neumann-Neumann data projected from T γ onto T i , i = 1, 2, then extracted the pressure trace on the interface and projected backward from T i onto T γ . Thus the interface problem is defined on the fracture time grid. Remark 5.1. In the nonconforming semi-discrete (in time) case, we see from (5.2) that P P P -1 N N is not strictly a preconditioner and may affect the accuracy of the scheme due the projection operators used to define it. This is indeed what we observed in the numerical experiments in Section 6.

For the GTO-Schwarz method. In the GTO-Schwarz method, there are two interface unknowns representing the linear combination of the fracture pressure and some terms from the fracture problem. Thus we let θ γ,j ∈ P 0 (T γ , L 2 (Γ)), for j = 1, 2. In order to obtain Ventcell boundary data for solving the subdomain problem (4.7), we project θ γ,j onto the T i , for i = 1, 2; j = 3 -i :

-u u u i • n n n i + αp i + s γ ∂ t p i + div τ u u u γ = Π iγ (θ γ,j ), on γ, i = 1, 2.
Remark 5.2. This setting is different from the case of usual domain decomposition (without fractures) analyzed in [START_REF] Hoang | Space-time domain decomposition methods for mixed formulations of flow and transport problems in porous media[END_REF][START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF], where the two interface unknowns represent the Robin data in each subdomain and thus are chosen to be constant on the associated subdomain's time grid, i.e. θ γ,i ∈ P 0 (T i , L 2 (Γ)), for i = 1, 2.

The semidiscrete-in-time counterpart of (4.11) is weakly enforced over each time sub-interval of the fracture time grid as follows: for all m = 0, . . . , M γ -1,

t m+1 γ t m γ θ γ,1 = t m+1 γ t m γ Π γ2 S VtR 2 (Π 2γ (θ γ,2 ), q 2 , p 0,1 , p 0,γ ) , t m+1 γ t m γ θ γ,2 = t m+1 γ t m γ Π γ1 S VtR 1 (Π 1γ (θ γ,1 ), q 1 , p 0,1 , p 0,γ ) , on γ, (5.3)
Remark 5.3. We point out that with the GTO-Schwarz method as with the GTP-Schur method preconditioned by a Neumann-Neumann preconditioner (cf. Remark 5.1), we can not hope to gain in accuracy in the fracture by using a finer grid in the fracture only since the fracture problem is actually solved on the coarser time grids of the two subdomains. We will see this in the numerical experiments.

6. Numerical results. In all of the numerical experiments, for the spatial discretization we use mixed finite elements with the lowest order Raviart-Thomas spaces on rectangles [START_REF] Boffi | Mixed finite elements methods and applications[END_REF][START_REF] Roberts | Mixed and hybrid methods[END_REF].

Remark 6.1. The subdomain problem of the GTO-Schwarz method corresponding to Ventcell boundary conditions is somewhat more complicated than that of GTP-Schur method (problem (3.1)). Consequently, for solving problem (4.6), one needs to introduce Lagrange multipliers (see e.g. [START_REF] Boffi | Mixed finite elements methods and applications[END_REF][START_REF] Roberts | Mixed and hybrid methods[END_REF]) on the interface to handle the Ventcell conditions (representing the fracture problem).

We carry out some preliminary experiments to investigate the numerical performance of the two methods proposed above. We consider the test case pictured in Figure 6.1 where the domain is a rectangle of dimension 2 × 1 and is divided into two equally sized subdomains by a fracture of width δ = 0.001 parallel to the y axis. The permeability tensors in the subdomains and in the fracture are isotropic: K K K = K i I I I, i = 1, 2, f, and K i is assumed to be constant. Here we choose K 1 = K 2 = 1 and K f = 10 3 (so that K f δ = 1). A pressure drop of 1 from the bottom to the top of the fracture is imposed. On the external boundaries of the subdomains a no flow boundary condition is imposed except on the lower fifth (length 0.2) of both lateral boundaries where a Dirichlet condition is imposed: p = 1 on the right and p = 0 on the left. See Figure 6.1.

We consider a uniform rectangular mesh with size h = 1/100. In time, we fix T = 0.5 and use uniform time partitions in the subdomains with time step ∆t i , i = 1, 2, and in the fracture with varying time step ∆t γ . We first consider the case with the same time step throughout the domain, ∆t 1 = ∆t 2 = ∆t γ = ∆t = T /300.

In Figure 6.2, snapshots at different times of the pressure field and of the flow field (on a coarse grid for visualization) are shown. The length of each arrow is proportional to the magnitude of the velocity it represents and the red arrows represent the flow in the fracture. We see that the flow field is a combination of flow in the fracture and flow going from right to left in the rest of the porous medium and there is interaction between them as some fluid flows out of the fracture (near the bottom) and some flows into it (near the top at later times). Since K f K i , i = 1, 2, the velocity is much larger in the fracture than in the surrounding medium.

Next, in order to analyze the convergence behavior of both methods, we consider the problem with homogeneous Dirichlet boundary conditions (i.e., the solution converges to zero). We start with a random initial guess on the space-time interface-fracture and use GMRES as an iterative solver and compute the error in the L 2 (0, T ; L 2 (Ω))-norm for the pressure p and for the velocity u u u. We stop the iteration when the relative error is less than 10 -6 . We consider four algorithms: the GTP-Schur method with no preconditioner, the GTP-Schur method with the local preconditioner, the GTP-Schur method with the Neumann-Neumann preconditioner and the GTO-Schwarz method with the optimized Robin parameter. We compare the convergence behavior of these four algorithms in terms of the number of iterations. Note however that for the GTP-Schur method with the Neumann-Neumann preconditioner the cost per iteration is roughly twice as large as that of the other methods.

In Figure 6.3, the error curves versus the number of iterations are shown: the error in p (on the left) and in u u u (on the right). We see that the GTP-Schur method with no preconditioner (the blue curves) converges extremely slowly (after 500 iterations, the error, both in p and in u u u, is about 10 -1 ). The performance of the GTP-Schur method with the local preconditioner (the green curves) is better but still quite slow: it requires about 350 iterations to reach an error reduction of 10 -6 . The Neumann-Neumann preconditioner (the cyan curves) further improves the convergence rate about 150 iterations are needed to obtain a similar error reduction. The GTO-Schwarz method needs only 6 iterations to reduce the error to 10 -6 and thus the convergence of the GTO-Schwarz method is much faster than that of the other algorithms (at least by a factor of 25). This is due to the use of the optimized parameter α. In Figure 6.4, we show the error in u u u (in logarithmic scale) after 10 Jacobi iterations for various values of α. We see that the optimized Robin parameter (the red star) is located close to those giving the smallest error after the same number of iterations. Also we observe that the convergence can be significantly slower if α is not chosen well.

Next, we study the behavior of three of the algorithms when nonconforming time grids are used. For this we again use the nonhomogeneous boundary conditions depicted in Figure 6.1. In all cases, we consider equal time steps for the subdomains as they have the same permeability: ∆t 1 = ∆t 2 = ∆t m . We examine three time grids as follows:

• Time grid 1 (conforming coarse): ∆t m = ∆t f = T /100.

• Time grid 2 (nonconforming): ∆t m = T /100 and ∆t f = T /500.

• Time grid 3 (conforming fine): ∆t m = ∆t f = T /500. We start with a zero initial guess on the space-time interface and stop the GMRES iterations when the relative residual is less than 10 -6 . In Figure 6.5 we show the relative residual versus the number of iterations for three schemes: the GTP-Schur method with the local preconditioner, the GTP-Schur method with the Neumann-Neumann preconditioner and the GTO-Schwarz method with an optimized Robin parameter. We see that the GTO-Schwarz method still performs better than the GTP-Schur method, and the GTP-Schur method with the Neumann-Neumann preconditioner still converges faster than with the local preconditioner. The convergence rate of both the GTP-Schur method with the Neumann-Neumann preconditioner and the GTO-Schwarz method are almost independent of the time grid (the number of iterations does not change with the time grid) while that of the local preconditioner significantly depends on the temporal grid in the fracture. We also notice that the behavior of all three methods in the cases of nonconforming and conforming fine grids are very similar.

Now we analyze the error (in time) of the three algorithms for each of the three time grids. A reference solution is obtained by solving problem (2.3) -(2.4) directly on a very fine time grid ∆t = T /2000. The L 2 -L 2 error of the difference between the multi-domain and the reference solutions at each iteration is computed. We distinguish two different errors: error in the rock matrix L 2 (0, T ; L 2 (Ω i )), i = 1, 2, and error in the fracture L 2 (0, T ; L 2 (γ)). Figures 6.6 and 6.7 show the pressure error in the subdomains and in the fracture respectively.

We first observe that the error in the subdomains after convergence (Figure 6.6) in the nonconforming case (Time grid 2) is equal to that in the conforming coarse case (Time grid 1) for all three algorithms. This is as expected as we use the same time step ∆t m = T /100 in the matrix for both of these grids. However, as already pointed out in Remark 5.3, though one might hope that the error in the fracture (Figure 6.7) in the nonconforming case is close to that in the conforming fine grid case (Time grid 3), this can only be the case for GTP-Schur method with the local preconditioner. Only for this case do we actually solve the fracture problem on the fine grid. For the other algorithms, the fracture error of the nonconforming case is equal to that of the conforming coarse grid instead (see Remark 5.1). However, none of the methods deteriorates the accuracy because of nonconforming time grids. Remark 6.2. While the GTO-Schwarz method does not make it particularly useful to use a finer time grid in the fracture, it does give a rather remarkable convergence speed. For the advection-diffusion problem with an explicit time scheme for advection, one of the main advantages of using smaller time steps in the fracture is to avoid imposing a time step in the two subdomains dictated by the CFL number of the equation in the fracture. Thus we are hopeful that this algorithm will be useful when coupled with the advection equation simply for the convergence speed that it gives. We add however that we are still pursuing some ideas for modifying this scheme to obtain an algorithm that can take advantage of smaller time steps in the fracture for the diffusion equation.

Conclusion.

We consider two domain decomposition methods for modeling the compressible flow in fractured porous media in which the fractures are assumed to be much more permeable than the surrounding medium. Two space-time interface problems are formulated using the time-dependent Dirichlet-to-Neumann and the Ventcell-to-Robin operators respectively, so that different time discretizations in the subdomains and in the fracture can be adapted. For the GTP-Schur method, two different preconditioners -the local and the Neumann-Neumann preconditioners-are considered and are first validated for a simple test case with one fracture. For the GTO-Schwarz method, the optimized parameter is used to accelerate the convergence of the associated iterative algorithm. Preliminary numerical experiments show that the GTO-Schwarz method converges much faster than the GTP-Schur method (with either the preconditioner) in terms of the number of iterations. The Neumann-Neumann preconditioner works better than the local preconditioner in the sense that its convergence is faster and is only weakly dependent on the mesh size of the discretizations. The GTO-Schwarz method also has a weak dependence on the mesh size. When nonconforming time steps are used, only the local preconditioner preserves the accuracy in time: the L 2 error in the fracture of the nonconforming time grid is close to that of the conforming fine grid. For the other algorithms, the L 2 error in the fracture of the nonconforming time grid is close to that of the conforming coarse grid instead. However, for the GTO-Schwarz method, this weak point when different time steps are used is compensated by the fast convergence of the algorithm.

Fig. 2 . 1 .

 21 Fig. 2.1. Left: The domain Ω with the fracture Ω f . Right: The domain Ω with the interfacefracture γ.

i) - 1

 1 , i = 1, 2, is the inverse of the operator ŠDtN i := S DtN i (•, 0, 0), and is defined by

4. 2 .

 2 Well-posedness of the subdomain problem with Ventcell boundary conditions. For a bounded domain O ⊂ R d (d = 2, 3) with Lipschitz boundary ∂O containing an open subset γ ⊂ ∂O, consider the following time-dependent problem written in mixed form with Dirichlet and Ventcell boundary conditions

Fig. 6 . 1 .

 61 Fig. 6.1. Geometry of the test case where the fracture is considered as an interface.

Fig. 6 . 2 .

 62 Fig. 6.2. Snapshots of the pressure field (left) and flow field (right) at t = T /300, t = T /4, t = T /2 and t = T respectively (from top to bottom).

Fig. 6 . 3 .

 63 Fig. 6.3. Convergence curves for the compressible flow: errors in p (on the left) and in u u u (on the right) -GTP-Schur method with no preconditioner (blue), with local preconditioner (green) and with Neumann-Neumann preconditioner (cyan) and GTO-Schwarz method (red).

Fig. 6 . 4 .

 64 Fig. 6.4. L 2 velocity error (in logarithmic scale) after 10 Jacobi iterations for various values of the Robin parameter. The red star shows the optimized parameters computed by numerically minimizing the continuous convergence factor.

Fig. 6 . 5 .

 65 Fig. 6.5. Relative residual with GMRES for different time grids: GTP-Schur method with the local preconditioner (green), GTP-Schur method with the Neumann-Neumann preconditioner (cyan) and GTO-Schwarz method (red).

Fig. 6 . 6 .

 66 Fig. 6.6. L 2 pressure error in the rock matrix: Time grid 1 (blue), Time grid 2 (magenta), Time grid 3 (black).

Fig. 6 . 7 .

 67 Fig. 6.7. L 2 pressure error in the fracture: Time grid 1 (blue), Time grid 2 (magenta), Time grid 3 (black).

§ Partially supported by ANDRA, the French agency for nuclear waste management * Partially supported by GNR MoMaS.

1

Appendix A. Proof of Theorem 1.2. We now give the proof of Theorem 1.2. The proof of the theorem is based on the Galerkin method, and its main steps will be given after the following lemma, which states the main energy estimates.

Remark A.1. The proof of Lemma A.2 is given in the infinite dimensional setting but some technical points (those involving u u u at time t = 0) can only be defined using their finite dimensional Galerkin approximations (as was done in detail for Dirichlet and Robin boundary conditions in [START_REF] Hoang | Space-time domain decomposition methods for mixed formulations of flow and transport problems in porous media[END_REF]). The results presented below have to be understood in that sense.

Lemma A.2. Under assumptions (H1), (H2) and (H4) above, the following a priori estimates hold.

where we recall that Σ a denotes the space Σ with the norm induced by the bilinear form a.

Proof. As usual we proceed by estimating successively p, u u u and ∂ t p.

• First, to derive an estimate for p, we take p(t) ∈ M and u u u(t) ∈ Σ as the test functions in (1.1) and add the two equations to obtain

Using the Cauchy-Schwarz inequality, we see that 1 2

Now integrating (A.2) over (0, t) for t ∈ (0, T ], we find

Then we use the non-negativity of c and apply Gronwall's lemma to obtain the first two estimates in (A.1)

• Next, to estimate ∂ t p, we differentiate the first equation of (1.1) with respect to t and take u u u as a test function. This yields

Then taking ∂ t p as a test function in the second equation of (1.1) wee see that

Now adding (A.4) and (A.5), we obtain

Integrating this inequality over (0, t) for t ∈ (0, T ], we have

where C c is the constant of continuity of the bilinear form c. There remains to bound the term u u u(0) 2 Σa . Toward this end, we use the first equation of (1.1) with v v v = u u u and for t = 0:

The (regularity) assumption that p 0 ∈ W enables us to write

and, as c(p(t), p(t)) ≥ 0, this gives the third inequality in (A.1)

• We now derive the last estimate. For this, we take µ = Bu u u as test function in the second equation of (1.1).

and the fourth inequality then follows by integrating in time and using the previous inequalities, which completes the proof of the lemma.

We now give the proof of the theorem. Proof. We first prove an estimate for u u u Σ , which follows easily from the second and fourth inequalities in Lemma A.2 and hypothesis (H3'): Note that this is the only place in the proof where (H3') was used. Lemma A.2 is independant of this hypothesis.

With the a priori estimates from Lemma A.2 and (A.8), the proof is concluded by using Galerkin's method.