Rita Riedlbeck
email: rita.riedlbeck@edf.fr

Daniele A Di

Alexandre Ern

Sylvie Granet

Kyrylo Kazymyrenko

Stress and flux reconstruction in Biot's poro-elasticity problem with application to a posteriori error analysis

Keywords: a posteriori error estimate, equilibrated stress reconstruction, Arnold-Winther finite element space, Biot's poro-elasticity problem

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Biot's poro-elasticity problem was originally proposed by von Terzaghi [START_REF] Terzaghi | Theoretical soil mechanics[END_REF] and Biot [START_REF] Biot | General theory of three-dimensional consolidation[END_REF] to describe the hydro-mechanical coupling between the displacement field u of a linearly elastic, porous material and the pressure p of an incompressible, viscous fluid saturating its pores. Let Ω ⊂ R 2 be the simply connected polygonal region occupied by the porous material, and let t F > 0 denote the simulation time. For the sake of simplicity, we assume that the material is clamped at its impermeable boundary, and fix the Biot-Willis coefficient equal to 1. We also assume that the deformation of the material is much slower than the flow rate, so that the problem can be considered in quasi-static form. Then, the displacement field u : Ω×(0, t F ) → R 2 1 and the pore pressure p : Ω × (0, t F ) → R are determined by

-∇ • σ(u) + ∇p = f
in Ω × (0, t F ), (1.1a)

∂ t (∇ • u + c 0 p) -∇ • (κ∇p) = g in Ω × (0, t F ), (1.1b)
u = 0 on ∂Ω × (0, t F ), (1.1c)

κ∇p • n Ω = 0 on ∂Ω × (0, t F ), (1.1d)
u(•, 0) = u 0 in Ω, (1.1e) where f denotes the volumetric body force acting on the material, g a volumetric fluid source (which, if c 0 = 0, is assumed to verify the compatibility condition Ω g(x, t)dx = 0 for each t ∈ (0, t F )), and the effective stress tensor σ is linked to the strain tensor through Hooke's law

σ(u) = 2µ (u) + λtr( (u))I 2 , (u) = 1 2 (∇u + ∇u T ), (1.2)
where I 2 is the two-dimensional identity matrix. The Lamé parameters λ and µ, describing the mechanical properties of the material, are assumed such that µ > 0 and λ + 2 3 µ > 0 uniformly in Ω. The scalar field κ : Ω → R describes the mobility of the fluid and we assume that there exist positive real numbers κ and κ such that κ ≤ κ ≤ κ a.e. in Ω. When the specific storage coefficient c 0 is zero, we enforce uniqueness of the pore-pressure in (1.1) by further requiring that Ω p(•, t)dx = 0 in (0, t F ).

(

On the other hand, for c 0 > 0, we complement (1.1) by the following initial condition on the pressure: p(•, 0) = p 0 in Ω. (1.4) In practice, even when c 0 = 0, the initial velocity field u 0 in (1.1e) is usually obtained by first setting p 0 equal to the solution of a hydrostatic computation and then calculating u 0 by solving (1.1a) with p = p 0 (cf. Remark 2.1 in [START_REF] Phillips | A coupling of mixed and continous Galerkin finite element methods for poroelasticity II: the discrete-in-time case[END_REF]). The well-posedness of Biot's consolidation problem has been analyzed in [START_REF] Showalter | Diffusion in poro-elastic media[END_REF][START_REF] Ženíšek | The existence and uniqueness theorem in Biot's consolidation theory[END_REF]. A suitable approximation method consists of using Taylor-Hood H1 -conforming finite elements in space (using piecewise polynomials of order k ≥ 1 for the pressure and of order (k + 1) for the displacement) and a backward Euler scheme in time. The corresponding a priori error analysis can be found in [START_REF] Murad | Improved accuracy in finite element analysis of Biot's consolidation problem[END_REF][START_REF] Murad | On stability and convergence of finite element analysis of Biot's consolidation problem[END_REF][START_REF] Murad | Asymptotic behaviour of semidiscrete finite-element approximations of Biot's consolidation problem[END_REF]. This discretization strategy is adopted in the Code Aster 1 software, which is used for the numerical examples presented in this work. Several other discretization methods have been studied in the literature, among which we cite, in particular, the fully coupled algorithm of [START_REF] Boffi | A nonconforming high-order method for the Biot problem on general meshes[END_REF], where the Hybrid High-Order method of [START_REF] Di Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF] is used for the elasticity operator, while the weighted discontinuous Galerkin method of [START_REF] Di Pietro | Discontinuous Galerkin methods for anisotropic semi-definite diffusion with advection[END_REF] is used for the Darcy operator.

The two governing equations (1.1a) and (1.1b) express, respectively, the conservation of mechanical momentum and fluid mass. In particular, the Darcy velocity φ(p) := -κ∇p and the total stress tensor θ(u, p) := σ(u) -pI 2 have continuous normal component across any interface in the domain Ω, and the divergence of these fields is locally in equilibrium with the sources (and the accumulation terms) in any control volume. It is well known that the use of H 1 -conforming finite elements does not lead to discrete fluxes φ(p n h ) and θ(u n h , p n h ) (where (u n h , p n h ) denotes the discrete solution at a given discrete time t n ) that satisfy the discrete counterpart of the above properties across mesh interfaces and in mesh cells. The first contribution of this work is to fill this gap by reconstructing equilibrated fluxes from local mixed finite element solves on cell patches around mesh vertices. The Darcy velocity reconstruction uses, as in [START_REF] Braess | Equilibrated residual error estimator for edge elements[END_REF][START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF], Raviart-Thomas mixed finite elements [START_REF] Raviart | A mixed finite element method for second order elliptic problems[END_REF] on cell patches around vertices of the original mesh. The construction we propose for the total stress tensor is, to our knowledge, novel and is based on the use of the Arnold-Winther mixed finite element [START_REF] Arnold | Mixed finite elements for elasticity[END_REF], again on the same vertex-based cell patches. This construction provides, in particular, a symmetric total stress tensor. The Darcy velocity and the total stress tensor are reconstructed at each discrete time, they have continuous normal component across any mesh interface, and their divergence is locally in equilibrium with the sources (averaged over the time interval) in any mesh cell. In steady-state linear elasticity, element-wise (as opposed to patch-wise) reconstructions of equilibrated tractions from local Neumann problems can be found in [START_REF] Ainsworth | Guaranteed computable error bounds for conforming and nonconforming finite element analysis in planar elasticity[END_REF][START_REF] Chamoin | An enhanced method with local energy minimization for the robust a posteriori construction of equilibrated stress field in finite element analysis[END_REF][START_REF] Ladevèze | Error estimation and mesh optimization for classical finite elements[END_REF][START_REF] Ohnimus | Local error estimates of FEM for displacements and stresses in linear elasticity by solving local Neumann problems[END_REF], whereas direct prescription of the degrees of freedom in the Arnold-Winther finite element space is considered in [START_REF] Nicaise | An a posteriori error estimator for the lamé equation based on H(div)-conforming stress approximations[END_REF].

The second contribution of this work is to perform an a posteriori error analysis of Biot's poro-elasticity problem using the above reconstructed fluxes to compute the error indicators. Equilibrated-flux a posteriori error estimates for poro-elasticity appear to be a novel topic (residual-based error estimates can be found, e.g., in [START_REF] Ern | A posteriori error analysis of Euler-Galerkin approximations to coupled ellipticparabolic problems[END_REF][START_REF] Meunier | Analyse d'erreur a posteriori pour les couplages hydro-Mécaniques et mise en oeuvre dans Code Aster[END_REF]). Equilibrated-flux a posteriori error estimates offer several advantages. On the one hand, error upper bounds are obtained with fully computable constants. The idea can be traced back to [START_REF] Prager | Approximations in elasticity based on the concept of function space[END_REF] and was advanced amongst others by [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF][START_REF] Braess | Equilibrated residual error estimator for edge elements[END_REF][START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF][START_REF] Ern | A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF][START_REF] Ladevèze | Comparaison de modèles de milieux continus[END_REF][START_REF] Ladevèze | Error estimate procedure in the finite element method and applications[END_REF][START_REF] Luce | A local a posteriori error estimator based on equilibrated fluxes[END_REF][START_REF] Repin | A posteriori estimates for partial differential equations[END_REF]. Another interesting property is the polynomial-degree robustness proved recently for the Poisson problem in [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]. A third attractive feature introduced in [START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF] is to distinguish among various error components, e.g., discretization, linearization, and algebraic solver error components, and to equilibrate adaptively these components in the iterative solution of nonlinear problems. This idea was applied to multi-phase, multi-components (possibly non isothermal) Darcy flows in [START_REF] Di Pietro | A posteriori error estimates, stopping criteria, and adaptivity for multiphase compositional Darcy flows in porous media[END_REF][START_REF] Di Pietro | An posteriori-based, fully adaptive algorithm for thermal multiphase compositional flows in porous media with adaptive mesh refinement[END_REF][START_REF] Di Pietro | Adaptive regularization, linearization, and discretization and a posteriori error control for the two-phase Stefan problem[END_REF]. For simplicity, we consider in the present work a global error measure which lends itself naturally to the development of equilibrated-flux error estimators, and defined as the dual energy-norm of the residual of the weak formulation. This paper is organized as follows. In Section 2, we introduce the weak and discrete formulations of Biot's poro-elasticity problem (1.1), along with some useful notation and preliminary results. In Section 3, we present the equilibrated reconstruction for the Darcy velocity and the total stress tensor. In Section 4, we derive a fully computable upper bound on the residual dual norm. We then distinguish two different error sources in the upper bound, namely the spatial and the temporal discretization, and we propose an algorithm adapting the mesh and the time step so as to equilibrate these error sources. Finally, we show numerical results in Section 5.

Setting

In this section we introduce some notation, the weak formulation, and the discrete solution of problem (1.1).

Weak formulation

We denote by L 2 (Ω), L 2 (Ω) and L 2 (Ω) the spaces composed of square-integrable functions taking values in R, R 2 and R 2×2 respectively, and by (•, •) and • the corresponding inner product and norm. We also let L 2 0 (Ω) := {q ∈ L 2 (Ω) | (q, 1) = 0}. H 1 (Ω) stands for the Sobolev space composed of L 2 (Ω) functions with weak gradients in L 2 (Ω) and H 1 0 (Ω) for its zero-trace subspace. H(div, Ω) and H(div, Ω) denote the spaces composed of L 2 (Ω) and L 2 (Ω) functions with weak divergence in L 2 (Ω) and L 2 (Ω), respectively, H s (div, Ω) the subspace of H(div, Ω) composed of symmetric-valued tensors, and

H 0 (div, Ω) := {ϕ ∈ H(div, Ω) | ϕ • n Ω = 0 on ∂Ω}.
We assume henceforth, for the sake of simplicity, that the volumetric body force f and the fluid source g lie in L 2 (0, t F ; L 2 (Ω)) and L 2 (0, t F ; L 2 0 (Ω)), respectively. In order to write a weak formulation of this poro-elastic problem, we define

U := H 1 0 (Ω), P := H 1 (Ω), (2.1)
where in the case c 0 = 0 we require additionally that P = H 1 (Ω) ∩ L 2 0 (Ω), and introduce the following Bochner spaces:

X := L 2 (0, t F ; U ) × L 2 (0, t F ; P ), (2.2a
)

Y := H 1 (0, t F ; U ) × H 1 (0, t F ; P ). (2.2b)
Let u, v ∈ U and p, q ∈ P . We define the bilinear forms

a(u, v) := (σ(u), (v)), (2.3a) b(v, q) := -(q, ∇ • v), (2.3b)
c(p, q) := (c 0 p, q), (2.3c) d(p, q) := (κ∇p, ∇q).

(2.3d)

Then, we consider the following weak formulation: find (u, p) ∈ Y , verifying the initial condition (1.1e) with u 0 ∈ H 1 0 (Ω) and (1.4) with p 0 ∈ H 1 (Ω) if c 0 > 0, and such that, for a.e. t ∈ (0, t F ),

a(u(t), v) + b(v, p(t)) = (f (t), v) ∀v ∈ U, (2.4a)
-b(∂ t u(t), q) + c(∂ t p(t), q) + d(p(t), q) = (g(t), q) ∀q ∈ P.

(2.4b)

The well-posedness of Biot's consolidation problem in slightly different weak formulations is shown in [START_REF] Showalter | Diffusion in poro-elastic media[END_REF][START_REF] Ženíšek | The existence and uniqueness theorem in Biot's consolidation theory[END_REF]. The uniqueness of the solution to (2.4) can be shown by energy arguments.

Assuming the existence of the solution in Y , we denote by σ(u) the resulting effective stress tensor, by θ(u, p) = σ(u) -pI 2 the total stress tensor and by φ(p) = -κ∇p the Darcy velocity. They verify the following properties:

θ(u, p) ∈ L 2 (0, t F ; H s (div, Ω)), -∇ • θ(u, p) = f, (2.5a)
φ(p) ∈ L 2 (0, t F ; H 0 (div, Ω)), ∇ • φ(p) = g -∂ t (∇ • u + c 0 p).
(2.5b)

Discrete setting

For the time discretization, we consider a sequence of discrete times (t n ) 0≤n≤N such that t i < t j whenever i < j, t 0 = 0, and t N = t F . For each 1 ≤ n ≤ N , let I n := (t n-1 , t n ) and τ n := t n -t n-1 . For a space-time function v, we denote v n := v(•, t n ) and define the backward differencing operator

∂ n t v = τ -1 n (v n -v n-1
). At each time step 1 ≤ n ≤ N , the space discretization is based on a conforming triangulation T n h of Ω, i.e. a set of closed triangles with union equal to Ω and such that, for any distinct T 1 , T 2 ∈ T n h , the set T 1 ∩ T 2 is either a common edge, a vertex or the empty set. We assume that T n h verifies the minimum angle condition, i.e., there exists α min > 0 uniform with respect to all considered meshes such that the minimum angle α T of each triangle T ∈ T n h satisfies α T ≥ α min . The set of vertices of the mesh is denoted by V n h ; it is decomposed into interior vertices V n,int h and boundary vertices V n,ext h . For any subdomain ω ⊂ Ω we denote V n ω the set of vertices in ω. For all a ∈ V n h , T n a is the patch of elements sharing the vertex a, and ω a the corresponding open subset of Ω. For all T ∈ T n h , V n T denotes the set of vertices of T , h T its diameter and n T its unit outward normal vector. For all n ∈ N and all k ∈ N, we denote by P k (T ) the space of bivariate polynomials in T ∈ T n h of total degree at most k and by P k (T n h ) = {ϕ ∈ L 2 (Ω) | ϕ| T ∈ P k (T ) ∀T ∈ T n h } the corresponding broken space over T n h . The following Poincaré's inequality holds for all T ∈ T n h :

v -Π 0,T v T ≤ C P,T h T ∇v T ∀v ∈ H 1 (T ), (2.6)
where Π 0,T : L 1 (T ) → P 0 (T ) is such that T (v -Π 0,T v)dx = 0 and C P,T = 1/π owing to the convexity of the mesh elements (see e.g. [START_REF] Bebendorf | A note on the Poincaré inequality for convex domains[END_REF]). Let

RM := {b + c(x 2 , -x 1 ) T | b ∈ R 2 , c ∈ R} (2.7)
denote the space of rigid body motions. We have the following Korn's inequality, again valid for all T ∈ T n h :

∇(v -Π RM,T v) T ≤ C K,T (v) T ∀v ∈ H 1 (T ), (2.8)
where Π RM,T :

H 1 (T ) → RM is such that T (v -Π RM,T v)dx = 0 and T rot(v -Π RM,T v)dx = 0 (with rot(v) := ∂ x 1 v 2 -∂ x 2 v 1 )
, and the constant C K,T is bounded by √ 2(sin(α T /4)) -1 (cf [START_REF] Kim | Guaranteed a posteriori error estimator for mixed finite element methods of linear elasticity with weak stress symmetry[END_REF]). Combining (2.6) and (2.8), and accounting for the bounds on the corresponding constants, we infer that, for all

T ∈ T n h , v -Π RM,T v T ≤ h T π √ 2 sin(α T /4) (v) T ∀v ∈ H 1 (T ).
(2.9)

Discrete problem

We will focus on the conforming Taylor-Hood finite element method using for each time step 1 ≤ n ≤ N the spaces

U n h := P k+1 (T n h ) ∩ U, P n h := P k (T n h ) ∩ P, (2.10)
with k ≥ 1. This method was first proposed in [START_REF] Taylor | A numerical solution of the Navier-Stokes equations using the finite element technique[END_REF] for incompressible flows and is known to provide stable pore pressure approximations (cf. [START_REF] Murad | Asymptotic behaviour of semidiscrete finite-element approximations of Biot's consolidation problem[END_REF] and [START_REF] Sandhu | Finite element analysis of seepage in elastic media[END_REF][START_REF] Yokoo | Finite element method applied to Biot's consolidation theory[END_REF]) and is a classical choice for the discretization of poro-mechanical problems by conforming finite elements.

Assumption 2.1 (Piecewise-constant-in-time source terms). For simplicity of exposition, we assume henceforth that the functions f and g are constant-in-time on each time interval I n and denote f n := f | In and g n := g| In .

Using the Taylor-Hood finite element spaces and a backward Euler scheme to march in time, the discrete problem reads: given

u 0 h and, if c 0 > 0, p 0 h , find (u n h , p n h ) ∈ U n h × P n h , for all 1 ≤ n ≤ N , such that a(u n h , v h ) + b(v h , p n h ) = (f n , v h ) ∀v h ∈ U n h , (2.11a) -b(∂ t u n hτ , q h ) + c(∂ t p n hτ , q h ) + d(p n h , q h ) = (g n , q h ) ∀q h ∈ P n h , (2.11b)
where we denote by u hτ , p hτ the discrete space-time functions which are continuous and piecewise affine in time, and such that, for each 0

≤ n ≤ N , (u hτ , p hτ )(•, t n ) = (u n h , p n h ), so that ∂ t u n hτ := ∂ t u hτ | In = τ -1 n (u n h -u n-1 h ) and ∂ t p n hτ := ∂ t p hτ | In = τ -1 n (p n h -p n-1 h ).
3 Quasi-static flux reconstructions

In contrast to (2.5), we have in general θ(u n h , p n h ) / ∈ H s (div, Ω) and φ(p n h ) / ∈ H 0 (div, Ω). In this section, we restore these properties by reconstructing H(div)-conforming discrete fluxes. These reconstructions are devised locally on patches of elements around mesh vertices. We first present the reconstructions in an abstract setting; then we apply these reconstructions to Biot's poro-elasticity problem. Since the time variable is irrelevant in devising the reconstructions, we drop the index n in the abstract presentation.

Darcy velocity

We reconstruct the Darcy velocity using mixed Raviart-Thomas finite elements of order l ≥ 0. For each element T ∈ T h , the local Raviart-Thomas polynomial spaces are defined by

W T := P l (T ) + xP l (T ), (3.1a)
Q T = P l (T ). (3.1b)
Figure 1 shows the corresponding degrees of freedom for l = 0 and l = 1. For each vertex a ∈ V h , the mixed Raviart-Thomas finite element spaces on the patch domain ω a are then defined as

W a h := {v h ∈ H(div, ω a ) | v h | T ∈ W T ∀T ∈ T a }, (3.2a) Qa h := {q h ∈ L 2 (ω a ) | q h | T ∈ Q T ∀T ∈ T a }. (3.2b)
We need to consider the following subspaces associated with the setting where a zero normal component is enforced on the velocity:

W a h := {v h ∈ W a h | v h • n ωa = 0 on ∂ω a }, (3.3a)
Q a h := {q h ∈ Qa h | (q h , 1) ωa = 0}. (3.3b)
The distribution of the degrees of freedom of functions in W a h is presented in Figure 2. Note that we enforce the zero normal condition also on patches associated with boundary vertices since a zero normal Darcy velocity is prescribed in the exact problem; see Remark 3.3 for other types of boundary conditions. For each a ∈ V h , let γ a ∈ L 2 (ω a ) be such that (γ a , 1) ωa = 0, and let Γ a ∈ L 2 (ω a ). Consider the following constrained minimization problem:

ϕ a h = argmin w h ∈W a h , ∇•w h =Π Q a h γa w h -Γ a ωa , (3.4)
where Π Q a h denotes the L 2 -orthogonal projection on Q a h . Then, extending ϕ a h by zero outside ω a , set

φ h := a∈V h ϕ a h . (3.5)
Since functions in W a h have zero normal component on ∂ω a , condition (γ a , 1) ωa = 0 is crucial for the well-posedness of the constrained minimization problem (3.4). This problem is classically solved by finding

ϕ a h ∈ W a h and s a h ∈ Q a h such that ϕ a h , w h ωa -(s a h , ∇ • w h ) ωa = (Γ a , w h ) ωa ∀w h ∈ W a h , (3.6a)
∇ • ϕ a h , q h ωa = (γ a , q h ) ωa ∀q h ∈ Q a h . (3.6b)
This problem is well-posed owing to the properties of mixed Raviart-Thomas finite elements, and we obtain the following result (cf. [START_REF] Braess | Equilibrated residual error estimator for edge elements[END_REF][START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]):

Lemma 3.2 (Properties of φ h ). Let φ h be prescribed by Construction 3.1. Then, φ h ∈ H 0 (div, Ω), and letting γ h ∈ L 2 (Ω) be defined such that γ h | T = a∈V T γ a for all T ∈ T h , we have

(γ h -∇ • φ h , q) T = 0 ∀q ∈ Q T ∀T ∈ T h . (3.7)
Remark 3.3 (Other boundary conditions). Suppose that we are given a partition of the boundary as ∂Ω = ∂Ω N,P ∪ ∂Ω D,P (the subsets ∂Ω N,P and ∂Ω D,P are conventionally closed in ∂Ω, i.e., ∂Ω N,P ∩ ∂Ω D,P is the common boundary of the two subsets) and that an inhomogeneous Neumann condition is enforced on the flux φ(p) on ∂Ω N,P (and a Dirichlet condition is enforced on p in ∂Ω D,P ). Assume that the mesh is fitted to the boundary partition, so that any mesh edge on the boundary belongs to either ∂Ω N,P or ∂Ω D,P . As detailed in [START_REF] Dolejší | hp-adaption driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF], this situation can be accommodated in Construction 3.1 up to minor modifications for all a ∈ V ext h (the construction is unmodified for all a ∈ V int h ). For the flux, we consider for the trial and test spaces, respectively,

W a h,N := {w h ∈ W a h | w h • n ωa | ∂ωa\∂Ω = 0, w h • n ωa | ∂ωa∩∂Ω N,P = Φ a,N }, W a h,0 := {w h ∈ W a h | w h • n ωa | ∂ωa\∂Ω = 0, w h • n ωa | ∂ωa∩∂Ω N,P = 0},
with Φ a,N related to the Neumann condition, whereas we set Q a h := Qa h if a lies on some edge in ∂Ω D,P and Q a h as in (3.3b) otherwise.

Total stress tensor

We reconstruct the total stress tensor using mixed Arnold-Winther finite elements of order m ≥ 1. One advantage of using these elements is that the reconstructed stress tensor is symmetric.

For each element T ∈ T h , the local Arnold-Winther polynomial spaces are defined by

Σ T := P s,m+1 (T ) + {τ ∈ P s,m+2 (T ) | ∇ • τ = 0} (3.8a) = {τ ∈ P s,m+2 (T ) | ∇ • τ ∈ P m (T )}, V T := P m (T ), (3.8b)
where P s,m (T ) denotes the subspace of P m (T ) composed of symmetric-valued tensors. Figure 3 shows the corresponding degrees of freedom in the cases m = 1 and m = 2. The dimension of V T is (m + 1)(m + 2), and it is shown in [START_REF] Arnold | Mixed finite elements for elasticity[END_REF] that dim(Σ T ) = (3m 2 + 17m + 28)/2. For the lowest-order case m = 1, the 24 degrees of freedom in Σ T are

• The values of the three components of the (symmetric) stress tensor at each vertex of the triangle (9 dofs);

• The values of the moments of degree 0 and 1 of the normal components of the stress tensor on each edge (12 dofs);

• The value of the moment of degree 0 of each component of the stress tensor on the triangle (3 dofs).

For each vertex a ∈ V h , the mixed Arnold-Winther finite element spaces on the patch domain ω a are defined as

Σa h := {τ h ∈ H s (div, ω a ) | τ h | T ∈ Σ T ∀T ∈ T a }, (3.9a) Ṽ a h := {v h ∈ L 2 (ω a ) | v h | T ∈ V T ∀T ∈ T a }. (3.9b)
We need to consider subspaces where a zero normal component is enforced on the stress tensor.

Since the boundary condition in the exact problem prescribes the displacement and not the normal stress, we distinguish the case whether a is an interior vertex or a boundary vertex (see Remark 3.6 for other types of boundary conditions). For a ∈ V int h , we set

Σ a h := {τ h ∈ Σa h | τ h n ωa = 0 on ∂ω a , τ h (b) = 0 ∀b ∈ V ωa ∩ ∂ω a }, (3.10a)
V a h := {v h ∈ Ṽ a h | (v h , z) ωa = 0 ∀z ∈ RM }, (3.10b)
and for a ∈ V ext h , we set Note that, as argued in [START_REF] Arnold | Mixed finite elements for elasticity[END_REF], the nodal degrees of freedom on ∂ω a are set to zero if the vertex separates two edges where the normal stress is enforced to be zero. The distribution of the degrees of freedom in Σ a h is presented in Figure 4. Construction 3.4 (Total stress tensor θ h ). For each a ∈ V h , let λ a ∈ L 2 (ω a ) be such that (λ a , z) ωa = 0 for all z ∈ RM and all a ∈ V int h , and let Λ a ∈ L 2 (ω a ). Consider the following constrained minimization problem:

Σ a h := {τ h ∈ Σa h | τ h n ωa = 0 on ∂ω a \∂Ω, τ h (b) = 0 ∀b ∈ V ωa ∩ (∂ω a \∂Ω)}, (3.11a)
V a h := Ṽ a h . ( 3
ϑ a h = argmin τ h ∈Σ a h , ∇•τ h =Π V a h λa τ h -Λ a ωa , (3.12)
where Π V a h denotes the L 2 -orthogonal projection onto V a h . Then, extending ϑ a h by zero outside ω a , set

θ h := a∈V h ϑ a h . (3.13)
The condition on λ a for all a ∈ V int h ensures that the constrained minimization problem (3.12) is well-posed. This problem is classically solved by finding ϑ a h ∈ Σ a h and s a h ∈ V a h such that

ϑ a h , τ h ωa -s a h , ∇ • τ h ωa = Λ a , τ h ωa ∀τ h ∈ Σ a h , (3.14a)
∇ • ϑ a h , v h ωa = (λ a , v h ) ωa ∀v h ∈ V a h . (3.14b)
This problem is well-posed (see [START_REF] Arnold | Mixed finite elements for elasticity[END_REF]), and we obtain the following result:

Lemma 3.5 (Properties of θ h ). Let θ h be prescribed by Construction 3.4. Then, θ h ∈ H s (div, Ω), and letting λ h ∈ L 2 (Ω) be defined such that λ h | T = a∈V T λ a for all T ∈ T h , we have

(λ h -∇ • θ h , v) T = 0 ∀v ∈ V T ∀T ∈ T h . (3.15)
Proof. All the fields ϑ a h are in H s (div, ω a ) and satisfy appropriate zero normal conditions so that their zero-extension to Ω is in H s (div, Ω). Hence, θ h ∈ H s (div, Ω). Let us prove (3.15). Let a ∈ V int h . Since (λ a , z) ωa = 0 for all z ∈ RM , we infer that (3.14b) actually holds for all

v h ∈ Ṽ a h . The same holds true if a ∈ V ext h by definition of V a h . Hence, ∇ • ϑ a h , v h ωa = (λ a , v h ) ωa for all v h ∈ Ṽ a
h and all a ∈ V h . Since Ṽ a h is composed of piecewise polynomials that can be chosen independently in each cell T ∈ T a , we conclude that (3.15) holds. all a ∈ V ext h (the construction is unmodified for all a ∈ V int h ): For the stress tensor, we consider for the trial and test spaces, respectively,

Σ a h,N := {τ h ∈ Σa h | τ h n ωa | ∂ωa\∂Ω = 0, τ h n ωa | ∂ωa∩∂Ω N,U = Θ a,N , τ h (b) = 0 ∀b ∈ V ωa ∩ (∂ω a \∂Ω), τ h (b) = θ a,N ∀b ∈ V ωa ∩ (∂Ω \ ∂Ω D,U )}, Σ a h,0 := {τ h ∈ Σa h | τ h n ωa | ∂ωa\∂Ω = 0, τ h n ωa | ∂ωa∩∂Ω N,U = 0, τ h (b) = 0 ∀b ∈ V ωa ∩ (∂ω a \∂Ω), τ h (b) = 0 ∀b ∈ V ωa ∩ (∂Ω \ ∂Ω D,U )},
with Θ a,N and θ a,N related to the Neumann condition, whereas we set V a h as in (3.11b) if a lies on some edge in ∂Ω D,U and V a h as in (3.10b) otherwise. Remark 3.7 (Extension to 3D). The extension of Construction 3.4 to three dimensions hinges on the existence of mixed finite element spaces producing three dimensional, H(div)-conforming, symmetric tensors. These were introduced in [START_REF] Arnold | Finite elements for symmetric tensors in three dimensions[END_REF], but are complex to implement and require significant computational effort, due to the high number of degrees of freedom per element (162 for the stress tensor).

Application to Biot's poro-elasticity problem

In this section, we apply the above constructions to the discrete Biot poro-elasticity problem (2.11). The reconstructed Darcy velocity and total stress tensor are space-time functions that are piecewise constant in time, i.e. these functions are calculated at every time step. We use Constructions 3.1 and 3.4 where we now specify the data γ a , Γ a , λ a , and Λ a for each 1 ≤ n ≤ N . For this purpose, we consider for any mesh vertex a ∈ V n h , the piecewise affine "hat" function ψ a ∈ P 1 (T h ) ∩ H 1 (Ω) supported in ω a , which takes the value 1 at the vertex a and zero at the vertices lying on the boundary of ω a ; cf. Figure 5. An important property for our constructions is the following partition of unity:

a∈V n T ψ a | T = 1 ∀T ∈ T n h . ( 3
γ a = ψ a g n -ψ a ∂ t (∇ • u hτ + c 0 p hτ ) n + ∇ψ a • φ(p n h ), Γ a = ψ a φ(p n h ), (3.17a)
λ a = -ψ a f n + θ(u n h , p n h )∇ψ a , Λ a = ψ a θ(u n h , p n h ), (3.17b)
where we recall that φ( The polynomial degree k in the Taylor-Hood finite element method (2.10) corresponds to degree k for the pressure and degree k + 1 for the displacement, implying polynomial degree k -1 for φ(p h ) and k for θ(u h , p h ). Thus, it seems reasonable that the polynomial degree for the reconstruction of the velocity is one lower than for the reconstruction of the total stress tensor. It has been shown in [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF] that k -1 or k are suitable choices for the velocity reconstruction, and for k = 1 we can observe the expected convergence rates for l = 0 and m = 1 in the numerical test of Section 5.2. Lemma 3.10 (Darcy velocity and total stress reconstructions). Construction 3.8 is well defined, and the following holds:

p n h ) = -κ∇p n h and θ(u n h , p n h ) = σ(u n h )-p n h I 2 . Then define φ n h ∈ H 0 (div, Ω) and θ n h ∈ H s (div,
(-∇ • θ n h , z) T = (f n , z) T ∀T ∈ T n h , ∀z ∈ RM, (3.18a)
(∇ • φ n h , 1) T = (g n -∂ t (∇ • u hτ + c 0 p hτ ) n , 1) T ∀T ∈ T n h . (3.18b)
Proof. Construction 3.1 is well defined provided (γ a , 1) ωa = 0 holds for all a ∈ V h , and this follows by taking q h = ψ a in (2.11b) (this is possible since hat functions are contained in the discrete space for the pressure). Construction 3.4 is well defined provided (λ a , z) ωa = 0 for all z ∈ RM and all a ∈ V n,int h , and this follows by taking v h = ψ a z in (2.11a) (recall that k ≥ 1 in (2.10), so that this choice is legitimate) and using that

(θ(u n h , p n h ), (ψ a z)) = (θ(u n h , p n h ), ∇(ψ a z) T ) = (θ(u n h , p n h ), z(∇ψ a ) T ) + (θ(u n h , p n h ), ψ a (∇z) T ) = (θ(u n h , p n h ), z(∇ψ a ) T ).
Finally, the properties on the divergence follow from Lemmas 3.5 and 3.

A posteriori error analysis and space-time adaptivity

In this section, we derive an a posteriori error estimate at every time step n based on the quasiequilibrated flux reconstructions of Section 3.3 for Biot's poro-elasticity problem. Using these estimators, we devise an adaptive algorithm including the adaptive choice of the mesh size and of the time step.

A posteriori error estimate

To derive the a posteriori error estimate, we consider the residual of the weak formulation (2.4). To combine the two equations in (2.4) into a single residual, the two equations must be written using the same physical units. Therefore, we introduce a reference time scale t and a reference length scale l which we will use as scaling parameters together with the Young

modulus E = µ(3λ + 2µ)(λ + µ) -1 . Defining the bilinear map B : Y × X → L 2 (0, t F ; R) such that B((u, p), (v, q))(t) := a(u, v)(t) + b(v, p)(t) + t (-b(∂ t u, q)(t) + c(∂ t p, q)(t) + d(p, q)(t)) , (4.1)
we can restate (2.4) as follows: find (u, p) ∈ Y , verifying the initial condition (1.1e), and (1.4) if c 0 > 0, and such that for a.e. t ∈ (0, t F ),

B((u, p), (v, q))(t) = (f, v)(t) + t (g, q)(t) ∀(v, q) ∈ X. (4.2)
For any pair (u hτ , p hτ ) ∈ Y , we define the residual R(u hτ , p hτ ) ∈ X of (4.2) as R(u hτ , p hτ ), (v, q) X ,X :=

t F 0 B((u -u hτ , p -p hτ ), (v, q))(t)dt.
Its dual norm is defined as

R(u hτ , p hτ ) X := sup (v,q)∈X, (v,q) X =1 R(u hτ , p hτ ), (v, q) X ,X , with (v, q) 2 X := t F 0 (E (v) ) 2 + (l ∇q ) 2 dt.
We first derive a local-in-time a posteriori error estimate. Let the time step 1 ≤ n ≤ N be fixed. Let us set

X n := L 2 (I n ; U ) × L 2 (I n ; P ) and Y n := H 1 (I n ; U ) × H 1 (I n ; P ),
and define the norm • Xn on X n in the same way as • X on X. The local error measure for the time step n is then defined as e n := sup (v,q)∈Xn, (v,q) Xn =1 In B((u -u hτ , p -p hτ ), (v, q))dt

= sup

(v,q)∈Xn, (v,q) Xn =1 In e n U (v) + e n P (q)dt,

e n P (q) := t In ∂ t (∇ • u + c 0 p) -∂ t (∇ • u hτ + c 0 p hτ ), q n (t) -φ(p) -φ(p hτ ), ∇q (t)dt, (4.4b)
where we recall that I n = (t n-1 , t n ) and that both u hτ and p hτ are continuous, piecewise affine functions in time.

For all 1 ≤ n ≤ N , let φ n h and θ n h be the constant in time fields over I n defined by Construction 3.8. For all T ∈ T n h , we define the residual estimators η n R,T,U , η n R,T,P by

η n R,T,U := h T π √ 2 sin(α T /4) E -1 f n + ∇ • θ n h T , (4.5a) η n R,T,P := h T π t l g n -∂ n t (∇ • u hτ + c 0 p hτ ) -∇ • φ n h T , (4.5b)
where α T denotes the minimum angle of the triangle T , and the flux estimators η n F,T,U (t), η n F,T,P (t), t ∈ I n , by

η n F,T,U (t) := E -1 θ n h -θ(
e n ≤   In T ∈T n h {(η n R,T,U + η n F,T,U (t)) 2 + (η n R,T,P + η n F,T,P (t)) 2 }dt   1 /2 . (4.7)
Proof. Let (v, q) ∈ X n . Recalling (4.4a), we have

e n U (v) = In θ(u, p) -θ(u hτ , p hτ ), (v) (t)dt = In f, v (t) -θ(u hτ , p hτ ), (v) (t)dt = In f n + ∇ • θ n h , v(t) T 1 (t) + θ n h -θ(u hτ , p hτ ), (v) (t) T 2 (t) dt, (4.8)
where we have used (2.4a) to pass to the second line and we have inserted (∇ • θ n hτ , v) + (θ n hτ , (v)) = 0 inside the integral to conclude. For the first term we have, for a.e. t ∈ (0, t F ),

|T 1 (t)| = T ∈T n h (f n + ∇ • θ n h , (v -Π K,T v)(t)) T ≤ T ∈T n h η n R,T,U E (v)(t) T ,
where we have used (3.18a) to insert Π K,T v inside the integral and (2.9) to conclude. For the second term, using the Cauchy-Schwarz inequality readily yields

|T 2 (t)| ≤ T ∈T n h θ n h -θ(u hτ , p hτ )(t) T (v)(t) T = T ∈T n h η n F,T,U (t)E (v)(t) T .
Inserting these results into (4.8) and applying the Cauchy-Schwarz inequality yields

|e U (v)| ≤   In T ∈T n h (η n R,T,U + η n F,T,U (t)) 2 dt   1 /2 × In (E (v)(t) ) 2 dt 1 /2 .
Proceeding in a similar way for e n P using (3.18b) and Poincaré's inequality (2.6) in place of (3.18a) and (2.9), respectively, we obtain

|e P (q)| ≤   In T ∈T n h (η n R,T,P + η n F,T,P (t)) 2 dt   1 /2 × In (l ∇q(t) ) 2 dt 1 /2
.

Combining these results and using again the Cauchy-Schwarz inequality yields

|e U (v) + e P (q)| ≤   In T ∈T n h (η n R,T,U + η n F,T,U (t)) 2 + (η n R,T,P (t) + η n F,T,P (t)) 2 dt   1 /2 × (v, q) Xn ,
and passing to the supremum concludes the proof.

Remark 4.2 (Data oscillation). Lemmas 3.2 and 3.5, and the mixed finite element space property

∇ • Σ n h (T ) = V n h (T ) and ∇ • W n h (T ) = Q n h (T ) for any T ∈ T n h imply η n R,T,U = h T π √ 2 sin(α T /4) E -1 f n -Π V n h (T ) f T , η n R,T,P = h T π t l g n -∂ n t (∇ • u hτ + c 0 p hτ ) -Π Q n h (T ) (g n -∂ n t (∇ • u hτ + c 0 p hτ )) T .
For the sake of convenience, we assumed the source terms f and g to be piecewise constant in time. When this is not the case, an additional data time-oscillation term appears in the right-hand side of the bound (4.7). ). The details of how they are obtained for the hydraulic part are shown in [START_REF] Dolejší | hp-adaption driven by polynomial-degree-robust a posteriori error estimates for elliptic problems[END_REF] and can be directly applied to linear elasticity. For each T ∈ T n h , let E N,U T and E N,P T be the set of edges lying on ∂Ω N,U and ∂Ω N,P respectively, and let θ n h and φ n h be the flux reconstructions of Remarks 3.3 and 3.6. Then we set

η n N,T,U = e∈E N,U T h T (2h e C t ) 1 /2 E sin(α T /4)|T | 1 /2 θ n h n Ω -θ N e , η n N,T,P = e∈E N,P T t h T (h e C t ) 1 /2 l |T | 1 /2 φ n h • n Ω -φ N e ,
where C t ≈ 0.77708, and (4.7) now reads

e n ≤   In T ∈T n h {(η n R,T,U + η n F,T,U (t) + η n N,T,U ) 2 + (η n R,T,P + η n F,T,P (t) + η n N,T,P ) 2 }dt   1 /2 .
To define a global-in-time a posteriori error estimate, we additionally define the initial condition errors η IC,T,U and η IC,T,P by setting

η IC,T,U := 1 2 E -1 t σ(u 0 -u hτ (•, 0)), (u 0 -u hτ (•, 0)) T 1 /2 , (4.9a)
η IC,T,P := 1 2 E -1 (t ) 2 c 0 ((p 0 -p hτ (•, 0)), p 0 -p hτ (•, 0)) T 1 /2 , (4.9b)
and we set

e IC = η IC =   T ∈T 0 h η 2 IC,T,U + η 2 IC,T,P   1 /2 . ( 4.10)
We define the global error as e := R(u hτ , p hτ ) X + e IC .

Corollary 4.4 (Global-in-time a posteriori error estimate). The following holds:

e ≤   N n=1 In T ∈T n h {(η n R,T,U + η n F,T,U (t)) 2 + (η n R,T,P + η n F,T,P (t)) 2 }dt   1 /2 + η IC . (4.12)
Proof. For each 1 ≤ n ≤ N , let ξ n ∈ X be the Riesz-representative of J n : X n → R with J n (v, q) = In B((u -u hτ , p -p hτ ), (v, q))dt. Then the function ξ ∈ X defined by ξ |In := ξ n will be the Riesz-representative of J :

X → R with (v, q) → t F 0 B((u -u hτ , p -p hτ ), (v, q))dt, so that R(u hτ , p hτ ) 2 X = J 2 X = ξ 2 X = N n=1 ξ n 2 Xn = N n=1 J n 2 X n = N n=1 (e n ) 2 . (4.13)
Inserting this result into (4.11) and applying Theorem 4.1 concludes the proof.

Distinguishing the space and time error components

The goal of this section is to elaborate the error estimate (4.7) so as to distinguish the error components resulting from the spatial and the temporal discretization. This is essential for the development of Algorithm 4.6 below, where the space mesh and the time step are chosen adaptively. Therefore, we add and subtract the discrete fluxes in the flux estimators (4.6) and apply the triangle inequality. We obtain, for all T ∈ T n h , the following local spatial and temporal discretization error estimators:

η n sp,T,U := η n R,T,U + E -1 θ n h -θ(u n h , p n h ) T , (4.14a)
η n sp,T,P := η n R,T,P +

t l φ n h -φ(p n h ) T , (4.14b)
η n tm,T,U (t) := E -1 θ(u n h , p n h ) -θ(u hτ , p hτ )(t) T , (4.14c)
η n tm,T,P (t) :=

t l φ(p n h ) -φ(p hτ )(t) T . (4.14d)
For each of these local estimators we can define a global version by setting

η n •,{U,P } (t) :=   2 In T ∈T n h η n •,T,{U,P } (t) 2 dt   1 /2 . ( 4.15)
Inserting them into (4.7) and applying the triangle inequality yields the following result.

Adaptive algorithm

Based on the error estimate of Theorem 4.5, we propose an adaptive algorithm where the mesh size and time step are locally adapted. The idea is to compare the estimators for the two error sources with each other in order to concentrate the computational effort on reducing the dominant one. Thus, both the spatial mesh and the time step are adjusted until space and time discretization contribute nearly equally to the overall error. For this purpose, let Γ tm > 1 > γ tm > 0 be user-given weights and crit n , for all 0 ≤ n ≤ N , a chosen threshold that the error on the time interval I n should not exceed. For each of the considered error sources we define the corresponding estimator as η • := η •,U + η In order to keep computational costs in the algorithm low, the initial mesh and time step should be chosen in a way that they match the criteria crit 0 and crit 1 . This can be achieved by performing only one time step before running the whole computation, and by modifying the initial discretization if they do not.

Numerical results

In this section we illustrate numerically our theoretical results on four test cases. For all tests we use the Taylor-Hood finite elements (2.10) with k = 1 and Construction 3.8 with l = 0 and m = 1. In the first two test cases, analytical solutions are known; the first one is a purely elastic, stationary problem and the second one a Biot's poro-elasticity problem. We analyze the convergence rates of the error estimators and compare them to those of an energy-type norm of the analytical error. The third test is the quarter five-spot problem, where we compare the results of the adaptive algorithm to a "standard" solution with fixed mesh and time steps.

In the fourth test, the excavation of two parallel tunnels is simulated. It shows an industrial application of the error estimators used for remeshing and again compares the performance of Algorithm 4.6 to a standard resolution.

Purely mechanical analytical test

For this stationary, purely mechanical test we consider the mode I loading of a cracked plate, corresponding to pure tension at the top and the bottom applied at the infinity. Following [START_REF] Williams | On the stress distribution at the base of a stationary crack[END_REF], an analytical solution around the crack tip is given by

u(r, θ) = 1 + ν E √ 2π √ r -cos( θ 2 )(3 -4ν -cos(θ)) sin( θ 2 )(3 -4ν -cos(θ)) , (5.1)
leading to a singularity of the stress tensor at the crack tip. For our test, we restrain ourselves to the domain Ω = (-

1 2 , 1 2 ) × (-1 2 , 1
2 ) with a straight crack from (-1 2 , 0) (cf. Figure 6), and impose the analytical solution (5.1) as Dirichlet boundary condition on ∂Ω and the crack faces to obtain the discrete solution u h . The Young modulus and the Poisson ratio are set to E = 1 and ν = 0, leading to the Lamé parameters µ = 0.5 and λ = 0. Since in this purely mechanical 7: Error estimation (left) and analytical error (right) on an initial mesh and after three mesh refinements test case there is no need for nondimensionalization, we omit the scaling factor E -1 in the error estimators. Figure 7 compares the distribution of the error estimators and the analytical error measured in the energy norm u -u h en = a(u -u h , u -u h ) 1 /2 . Besides detecting the dominating error at the crack tip due to the singularity of σ(u), the error estimators reflect the distribution of the analytical error in the whole domain, as can be seen in the lower panel for the finer mesh.

Poro-elastic analytical test

Let Ω = (0, 1)×(0, 1). Following [START_REF] Boffi | A nonconforming high-order method for the Biot problem on general meshes[END_REF][START_REF] Ern | A posteriori error analysis of Euler-Galerkin approximations to coupled ellipticparabolic problems[END_REF], we consider the analytical solution of Biot's consolidation problem (1.1) u(t, x, y) = cos(-πt) cos(πx) sin(πy) sin(πx) cos(πy) , p(t, x, y) = sin(-πt) sin(πx) sin(πy), with κ = 1, c 0 = 0, and the Lamé coefficients µ = λ = 0.4, yielding a Young modulus E = 1 and a Poisson ratio ν = 0.25. The resulting source terms are given by f (t, x, y) = (2.4π 2 cos(-πt) + π sin(-πt)) cos(πx) sin(πy) sin(πx) cos(πy) , and g = 0.

To evaluate convergence rates under space or time uniform refinement, we measure the analytical error in the energy norm .3) where in this dimensionless test, the nondimensionalization parameters t and l are both equal to one, and we also omit the factor E.

|||(v, q)||| 2 en = t F 0 B((v, q), (t ∂ t v, q))dt = 1 2 t (a(v, v)(t F ) -a(v, v)(t 0 ))+t t F 0 d(q, q)dt, ( 5
v 2 U = t F 0 a(v, v)dt and q 2 P = t F 0 d(q, q)dt. ( 5
Tables 1 and 2 compare the convergence rates under space and time refinement of the corresponding error estimators to the analytical error in the norms defined by (5.2) and (5.3). The last column shows the effectivity index defined by I eff := η sp,U + η sp,P + η tm,U + η tm,P |||(u -u hτ , p -p hτ )||| en .

(5.4)

For both the spatial and the temporal refinement, we obtain the expected convergence rates of the Taylor-Hood finite element method (2.10) with k = 1, and a backward Euler scheme in time. The last value of the analytical mechanical error under space refinement is due to the error in time discretization which starts playing a role for the finest mesh. We observe that the orders of magnitude of the mechanical and hydraulic part are comparable in this test, and that the effectivity index is dominated by the hydraulic part under space refinement, and by the mechanical part under time refinement.

Quarter five-spot problem

In this standard configuration considered in petroleum engineering, the injection of water at the center of a square domain and the production at the four corners is simulated on a quarter of the domain. In our test, this quarter is a square of 100m side length, divided into two parts with different mobilities; a circle around the injection point of radius 50m with κ = 8•10 -9 m 2 Pa -1 s -1 , and κ = 10 As shown in the first test, the stress reconstruction and a posteriori estimate presented in this work are directly applicable to pure linear elasticity problems. The second test shows that the presented error estimate also delivers sharp bounds (as reflected by moderate effectivity indices) of more accessible error measures computed using energy-type norms. In the third and fourth tests, comparing the proportions of the estimators for the hydraulic and the mechanical parts reflects the physical properties of the problem: in the quarter five-spot test, the dominating estimators are those for the hydraulic part; for the excavation damage test, they are approximately of the same order of magnitude, with the mechanical estimator dominating in regions of stress concentration.

Figure 1 :

 1 Figure 1: The mixed Raviart-Thomas finite element for l = 0 (left) and l = 1 (right)

Figure 2 :

 2 Figure 2: The degrees of freedom of the space W a h in the case l = 0 on a patch for a ∈ V int h

Figure 3 :

 3 Figure 3: Element diagrams for the pair (Σ T , V T ) in the cases m = 1 (left) and m = 2 (right)

Remark 3 . 6 (∂ΩFigure 4 :

 364 Figure 4: The degrees of freedom of the space Σ a h in the case m = 1 on a patch for a ∈ V int h

Figure 5 :

 5 Figure 5: Hat function

  2 and the partition of unity (3.16) which implies that a∈V n T ∇ψ a | T = 0. Remark 3.11 (Other boundary conditions). If inhomogeneous Neumann boundary conditions φ(p) • n Ω = Φ N or θ(u, p)n Ω = Θ N are imposed on ∂Ω N,P and ∂Ω N,U , respectively, we modify the constructions using Remarks 3.3 and 3.6. In particular, in Remark 3.3, Φ a,N is the L 2projection of ψ a Φ N onto W a h • n Ω , whereas in Remark 3.6, Θ a,N is the L 2 -projection of ψ a Θ N onto Σa h n Ω .

  u, p) -θ(u hτ , p hτ ), (v))(t)dt,(4.4a)

Remark 4 . 3 (

 43 Other types of boundary conditions). If we consider inhomogeneous Neumann boundary conditions θ(u, p)n Ω = θ N on ∂Ω N,U ⊆ ∂Ω and φ(p) • n Ω = φ N on ∂Ω N,P ⊆ ∂Ω, two more error estimators appear in (4.7

Theorem 4 . 5 (

 45 A posteriori error estimate distinguishing the error components). Let 1 ≤ n ≤ N . Let (u, p) be the weak solution of (2.4) and let (u hτ , p hτ ) ∈ Y n be the discrete solution of(2.11). Let θ n h , φ n h be the equilibrated fluxes of Construction 3.8. Then the following holds for the error e n defined by (4.3) with estimators defined by (4.14) and (4.15): e n ≤ η n sp,U + η n sp,P + η n tm,U + η n tm,P . (4.16)

Figure 6 :

 6 Figure 6: Loading of a cracked plate

Figure

  Figure 7: Error estimation (left) and analytical error (right) on an initial mesh and after three mesh refinements

-9 m 2

 2 Pa -1 s -1 in the rest of the domain. The Young modulus and the Poisson ratio are given by E = 10 9 Pa, ν = 0.3, and we set c 0 = 0. The initial state is given by θ 0 = 0, φ 0 = 0 and p 0 = 10 5 Pa. During the computation time of 30 days, we set p = p 0 in the top right corner and p = 4 • 10 5 Pa in the bottom left corner, simulating the production and the injection respectively. The nondimensionalization parameters are l = 140m and t = 1h. The problem is dominated by hydraulic processes.

Figure 9 :

 9 Figure 9: Initial mesh (left) and meshes at the end of the first (center) and the second (right) excavation in the adaptive algorithm

Figure 10 :

 10 Figure 10: Spatial discretization error estimators at t F on a fixed mesh (left, 29,275 dofs) and on the last mesh of an adaptive algorithm (right, 15,064 dofs)

  u hτ , p hτ )(t) T , Theorem 4.1 (Local-in-time a posteriori error estimate). Let (u, p) ∈ Y be the weak solution of (2.4) and let (u hτ , p hτ ) ∈ Y be the discrete solution of(2.11). Let 1 ≤ n ≤ N . Let e n be defined by (4.3) with estimators defined by (4.5) and (4.6). Then the following holds:

				(4.6a)
	η n F,T,P (t) :=	t l	φ n h -φ(p hτ )(t) T .	(4.6b)

  Choose an initial triangulation T 0 h , an initial time step τ 0 , and set t 0 := 0 (b) Initial mesh adaptation loop • Calculate η IC,U and η IC,P • Refine or coarsen the mesh T 0 h such that the local initial condition error estimators η IC,T,• are distributed equally End of the loop if η IC ≤ crit 0

	Owing to (4.11), (4.13), (4.17) and (4.18) the obtained discrete solution satisfies
		N	1 /2
	e ≤	(crit n ) 2	+ crit 0 .	(4.19)
		n=1
			•,P , so that (4.16) becomes
		e n ≤ η n sp + η n tm		(4.17)
	Algorithm 4.6 (Adaptive algorithm).
	1. Initialisation
	(a) 2. Time loop
	(a) Set n := n + 1, T n h := T n-1 h	, and τ n := τ n-1
	(b) Calculate (u n h , p n h ) and the estimators η n sp and η n tm
	(c) Space refinement loop
	i. Space and time error balancing loop
	A. if γ tm η n sp > η n tm : Set τ n := 2τ n
	B. if Γ tm η n sp < η n tm : Set τ n := 1 2 τ n
	End of the space-time error balancing loop if
	γ tm η n sp ≤ η n tm ≤ Γ tm η n sp	or	τ n ≤ τ min
	ii. Refine or coarsen the mesh T n h such that the local spatial error estimators η n sp,T
	are distributed equally
	End of the space refinement loop if
		η n sp + η n tm ≤ crit n	(4.18)
	End of the time loop if t n ≥ t F

Table 1 :

 1 Error estimators and analytical errors under space refinement with t F = 0.5, τ = 5e-5

	.2)

Table 2 :

 2 Error estimators and analytical errors under time refinement with t F = 0.5, h = 1/128 and the mechanical and hydraulic parts separately in the following norms:

http://web-code-aster.org

Acknowledgements

The work of D. A. Di Pietro was supported by ANR grant HHOMM (ANR-15-CE40-0005) References

 3. Right: An adapted mesh at time t = 10 days

We use Algorithm 4.6 to perform space-time adaptivity. We start with an initial mesh of 10,638 vertices and with an initial time step of τ 0 = 12h. For the space-time error balancing, we set γ tm = 0.8 and Γ tm = 1.3 and fix the error limit for each time step to crit n = 0.005τ n . We compare the performance of the adaptive algorithm to two static computations (i.e. with fixed meshes and time steps), one, called equivalent, where the discretization is chosen in a way to have approximately the same number of space-time unknowns as in the adaptive algorithm, and one where the discretization is very fine, so its solution can be taken as a reference solution. Table 3 compares the number of space-time unknowns and performed iterations (i.e. the number of time steps, counting repetitions in the adaptive algorithm), and the values of the error estimators of the three computations.

The left graphic in Figure 8 shows the discrete pressure along the diagonal going from the bottom left to the top right of the domain (as indicated in the right graphic) at three different times obtained by the static computations (solid and dotted lines) and the adaptive algorithm (dashed lines). The loosely dotted vertical line marks the edge between the two parts of Ω with different permeabilities. At each of these times, the discrete solution of the adaptive algorithm is closer to the reference solution than the equivalent computation using a fixed mesh and time step. At the last time step, all the results get closer as the solution converges in time to a constant state. The performance of Algorithm 4.6 is tested on four different initial meshes with crit n = 7•10 -3 τ n for the coarsest one and, with the mesh getting finer, crit n = 4 • 10 -3 τ n , crit n = 2 • 10 -3 τ n and crit n = 1 • 10 -3 τ n . In all the calculations we fix γ tm = 0.8, Γ tm = 1.5 and τ 0 = 3.9d. Figure 9 illustrates the evolution of the second coarsest mesh with crit n = 4 • 10 -3 τ n . During the first excavation, the refinement takes only place around the left tunnel, whereas the area around the right tunnel is only refined after the beginning of the second excavation. The calculations resulting from the adaptive algorithm are compared to calculations with fixed meshes and time steps. Each of these meshes is slightly finer around the tunnels than in the rest of Ω, and the time steps are chosen in a way that η sp ≈ η tm . Figure 10 compares the spatial discretization error estimators at the final time t F of the static algorithm to those of the adaptive algorithm, which are much more evenly distributed over the domain. Furthermore, the left graphic in Figure 11 shows that in our test, the use of the adaptive algorithm reduces the number of space-time-unknowns for a similar value of the error estimator.

In the right graphic of Figure 11, we plot the evolution of the error estimators in the two computations circled in the left graphic. Each mark stands for an iteration and shows the error estimate e n of the current time interval divided by τ n . For the plain algorithm, an iteration is equal to a time step. The adaptive algorithm recalculates the solution at a time step whenever the error estimate lies over crit n (illustrated by the dashed line) by refining τ n or the mesh (or both). Thus, only the square shaped points in the graphic contribute to the overall error estimate. In the consolidation phase between the two excavations (from t = 17.4 days to t = 29 days), the mesh is slightly coarsened and the time step considerably increased, since the dominating error source in this phase is the spatial discretization.

Conclusion

The analytical test cases show that the distribution and convergence rates of our error estimators reflect those of the analytical error. The efficiency of Algorithm 4.6 has been illustrated in industrial tests, where the number of space-time unknowns is considerably decreased for a comparable overall error estimate. We also observe that the price for computing the flux reconstructions can be substantially reduced by pre-processing, a task that is fully parallelizable.


C Legrand

J Beugin

E.-M El-Koursi

J Marais

M Berbineau

Approach for evaluating the safety of a satellite-based train localisation system through the extended integrity concept

The integrity concept, safety quality criterion for satellite-based localisation systems used in aeronautics, is described in terms of levels (protection and alert levels), time (Time To Alarm) and probability (integrity risk). In land transport applications, the requirements in terms of integrity differ from aeronautics in their definition and values. Global Navigation Satellite Systems (GNSS) in railways suffer from additional weaknesses i.e. multipath and masking phenomena, which can degrade the localisation integrity. This situation cannot be tolerated in safety-related applications like train control and signalling. To mitigate these weaknesses, GNSS is usually combined with other localisation systems like inertial sensors. However, existing integrity monitoring processes are designed for GNSS integrity evaluation, i.e. to estimate the risk allocated to the position given by the GNSS receivers only. Our research work aims, first, to extend the integrity concept to such systems, and, secondly, to demonstrate how to evaluate, with this concept, the safety of a localisation system as expected in railways. The safety of GNSS-Based Localisation System is formalised and quantitatively evaluated.

INTRODUCTION

Global Navigation Satellite Systems (GNSS) are now well-spread in air and road transports for non-safety-relevant applications as, as well-known, fleet management, driver assistance or passenger information. In these services, the quality of the localisation information is not a vital parameter regarding the safety [START_REF] Beugin | Simulation-based evaluation of dependability and safety properties of satellite technologies for railway localization[END_REF] and the performances provided by mass market GNSSbased systems are sufficient, as the one embedded in smartphones.

To guarantee their high level of safety, railways are ruled by standards (especially EN50126 (CENELEC 2007)) that require the degree of confidence, which the user can place in the delivered service through a dependability evaluation. Especially, a safety demonstration is required to put the satellite-based train localisation system into service and has to be delivered in safety documentation required by certification bodies and national safety authorities. GNSS-based positioning systems will have to go through these steps before being generalised. As classical done in robotics or other vehicular applications, some of these railway consortia as 3inSat [START_REF] Salvatori | Augmentation and integrity monitoring network and EGNOS performance comparison for train positioning[END_REF]) developed hybridised GNSS-based multi sensor systems based on different choice of sensors and different strategies to ensure safety but exploiting the advantages and drawbacks of each sensor used. This paper will focus on the safety evaluation of these kind of systems.

GNSS performance evaluation does not only rely on accuracy when speaking about railway safety. Integrity, linked to a confidence granted to the system, is of prime importance (Le [START_REF] Marchand | Vehicle localization integrity based on trajectory monitoring[END_REF][START_REF] Liu | Research on integrity monitoring for integrated gnss/sins system[END_REF][START_REF] Ochieng | User level integrity monitoring and quality control for high accuracy positioning using gps/ins measurements[END_REF]. This concept permits us to estimate the risk linked to the position error exceeding an alert limit without being detected called also the integrity risk. However, initially expressed only for GNSS, the integrity must be extended to the positioning systems except GNSS. Thanks to the extended integrity definition and the use of a risk probability related to the railway safety requirements, it is possible to formalise the link between the probability of this risk and the integrity risk.

The paper is organised as follow. In section 2, each part of the considered satellite localisation system hybridised is described in order to make explicit the potential sources of hazard and, finally, the GNSS integrity concept is introduced. In the section 3, the extended integrity of a satellite localisation system hybridised is defined. The section 4 to propose an approach for evaluating the safety using the extended integrity concept thanks to a mathematical formalism of the link between integrity and safety. Paper conclusions are given in section 5.

RISK IN SATELLITE-BASED TRAIN LOCALISATION SYSTEM

In its ERTMS (European Rail Traffic Management System) Memorandum of Understanding, the European Railway Agency indicates that GNSS particularly the future Galileo (associated with the satellite augmentation system) can play a major role in the rail sector, both for fleet management and rail safety (signalling and train control) (European Railway Agency 2012). ERTMS, particularly, the ETCS (for European Train Control System) is a new train control system and signalling designed to replace the 27 existing systems in Europe, which generate interoperability problems. In an ERTMS context, satellite technologies can contribute to reducing the costs of the infrastructure (for example, the balise) and enhancing the performances of the ETCS odometer in new trains without impact on equipped lines.

Global Navigation Satellite System risk and risk mitigation

GNSS technology is based on the trilateration concept, which defines the process of determining absolute position by measurement of distances, called pseudoranges using spheres. Typically, these pseudoranges are the distances between a user receiver and satellites. For a single satellite signal, a user can be located on the surface of a sphere of a radius r centred on that satellite. With two signals, the user position can be located on the intersection of two spheres of radii r 1 and r 2 . It is necessary to use an additional signal (a third satellite is requested) in order to limit the user position to two points. For terrestrial applications, one of these two points is possible (the other point is in space). However, the user receiver and the satellite clocks are not exactly synchronised due to relativist effects and lower quality of the receiver clock compared to the satellites' clock.

In the three-equations system, the synchronisation issue adds a fourth unknown variable, the clock offset.

GNSS standalone receiver is particularly efficient in terms of availability and safety in open sky areas and precise in the long term [START_REF] Lu | Performance evaluation of GNSS for train localization[END_REF]. However, it suffers from multiple weaknesses, such as multipath effect, signals blocked by environmental items (trees, buildings, terrain relief). These signals are also disrupted by RF interference. Several mechanisms on board (use of very precise receivers, map-matching method, etc...), on the ground (relay masts, ground-based augmentation system, etc...) and at the level of satellites (use of multi-constellation of satellites (GPS (American) + Glonass (Russian) + Galileo (European)), satellite-based augmentation system, etc...) exist to mitigate these weaknesses. Besides, techniques for detecting and correcting faults can be implemented in a localisation system (integrity monitoring -integrity concept is described in subsection 2.2). This last solution is the only one on which railway stakeholders can act on, so it is chosen in this paper. Subsection 2.1.2 deals with the different kind of integration. Before explaining this, the inertial systems that can be hybridised to GNSS are introduced, especially their errors.

Dead reckoning and risk

Dead reckoning is a computation process of the current position by using a previously determined position. The Inertial Navigation Systems (INS) and odometer are examples of navigation systems using dead reckoning. First, a INS [START_REF] Woodman | An introduction to inertial navigation[END_REF]) are commonly used in air, road transport and railways. These systems are composed of two hardware and software parts and provide a navigation solution i.e. position, speed and acceleration. The acceleration is determined by the measurement of the specific force and angular rate in body frame (roll, pitch and yaw). The hardware part corresponds to the sensor part called Inertial Measurement Unit (IMU), which is composed of three accelerometers providing an acceleration measurement and three gyroscopes providing a rotation (angular position) measurement (one for each axis in a inertial frame). The software part corresponds to a computer unit, which resolves inertial equations in order to give a navigation solution. At this level, the measurements coming from accelerometers and gyroscopes are merged. Second, a odometer is a device fitted on the bogie axles, which provides a distance travelled by a vehicle. Classically, it is composed of an incremental encoder, which measures elementary motions of the vehicle. These systems suffer from slipping phenomena (odometer) and gyroscope bias and initial pitch error (INS), which lead to cumulative errors at each time. Further, these errors are called Slowly Growing Errors showing the slow and insidious nature of his errors. At their occurrence, Slowly Growing Errors are not enough significant to be detected. Without calibration means (note that, in railway context, this calibration is performed by balises, device regularly placed on a railway), these Slowly Growing Errors become too high to be tolerated for all the duration of a mission.

GNSS/INS integration and risk

A standalone GNSS receiver associated with an augmentation system does not meet railway safety requirements [START_REF] Rispoli | Recent progress in application of GNSS and advanced communications for railway signaling[END_REF]. The strategy of hybridisation appears as a realistic solution in difficult environments [START_REF] Gioia | Stand-alone and hybrid positioning using asynchronous pseudolites[END_REF]. The hybridisation is a combination of different technologies of localisation (odometry, inertial system, etc.), which provides additional navigation signals. Different kinds of hybridisation exist (loosely, tightly, deeply coupled architecture) according to the input data chosen [START_REF] Groves | Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems[END_REF]). These integration architectures will not be fully developed in this paper. It is important to know that the more the input data will be taken upstream of the GNSS subsystem, the more the architecture will be tightly coupled. In this case, these input data are the pseudoranges and not a navigation solution computed by a software part. The main advantage of the INS/GNSS architecture is that each part of the system mitigates the weaknesses others'. Indeed, a proprioceptive localisation technique associated with an exteroceptive one is particularly efficient in several aspects (accuracy, availability, safety, etc.). The use of two identical localisation techniques is less judicious: they suffer from the same phenomena and they can mitigate their effects. In the case of inertial/GNSS integration, the GNSS accuracy in the long term mitigates the Slowly Growing Errors affecting the inertial system. So, in order to keep the capacity for excluding a localisation source at the occurrence of failure(s), a loose architecture is not an optimal solution. In that case, if a failure is detected on the GNSS side, the navigation solution, which it provides, is entirely excluded and the train continues to operate with only INS, a situation to avoid in the long term (cf subsection 2.1.1). In consequence, a tightly coupled GNSS/INS is more resilient to failures. Instead of excluding all the GNSS navigation solution in the failure case, it is possible to exclude only one pseudorange on the condition that the number of satellite is higher than 4, the minimum to determine a navigation solution (cf subsection 2.1). Knowing that more than four satellites are in view at any point of earth, GNSS receiver is capable of computation a navigation solution in spite of failures.

To put into service this kind of system, it is necessary to perform a dependability evaluation through reliability, availability, maintainability and safety (RAMS). However, GNSS performances are evaluated in terms of another class of attributes: accuracy, availability, continuity and integrity (Beugin and Marais 2012). In the next part, we propose to focus on the integrity attribute because, in the aeronautic domain, the integrity is well defined and several data processing systems to evaluate it, exist (For example, the RAIM (Receiver Automonious Integrity Monitoring) algorithms [START_REF] Liu | Research on integrity monitoring for integrated gnss/sins system[END_REF]). In railways, dependability methods are few [START_REF] Lu | Performance evaluation of GNSS for train localization[END_REF] for RAMS assessment of GNSS and we propose to counterbalance this observation by revision of integrity (definition and mechanisms) in order to evaluate the safety of these systems.

GNSS integrity concept

The integrity concept for GNSS is derived from aeronautical recommendations. It defines "the measurement of the degree of confidence that can be placed in the correctness of the navigation information". It includes "the ability of a localisation system to provide timely warnings to the user when the system should not be used for navigation" [START_REF] Peters | Federal radionavigation plan[END_REF]. Aeronautics rules GNSS localisation systems by aeronautical standards (ICAO 2006) especially in terms of integrity according to the following requirements:

• Alert Limit (AL): the maximum allowable position error beyond which the localisation system should be declared unavailable for the intended application

• Time To Alert (TTA): the maximum allowable time elapsed from the onset of the navigation system being out of tolerance until the equipment enunciates the alert

• Integrity Risk: risk related to the position error exceeding the alert limit without the user being informed under TTA period of time Integrity requirements are given for each application/operation/mission. However, to verify their achievement, the position error has to be known and is function of the true position, an unobservable quantity for the user. In consequence, it is necessary to use • Protection Level (PL): information guaranteed to the user. It is a variable quantity and computed online from measurement estimations (cf subsection 3.2)

• Uncertainty Level (UL): interval around true position. The uncertainty level is a variable dependon current measurements.

To apply the concept to railways, the effort shall be put on adaptation to the terrestrial transport specificities and to railway rules. In this objective, some automotive (Le Marchand 2010) and railway [START_REF] Jonas | Integrity enhancement of the GNSS position solution for the railway applications[END_REF]) publications attempted to enhance integrity of position solution but it is always centred on GNSS. In a hybridised GNSS, the other positioning systems other than satellite-based ones are considered failurefree.

Looking at the integrity definition, this concept is compatible with other transport modes other than aeronautics except for the specific parameters UL, PL and the triplets of requirements (AL, TTA and integrity risk). In one hand, vertical error components of AL/UL/PL are less significant and can be neglected for rail and road applications. On the other hand, the integrity requirements depend on the application and the phase of a mission obviously very different according to the transport domain. In addition, these parameters do not consider the Slowly Growing Errors and the errors related to the masking or multipath effects on satellite signals due to varied environments encountered by trains. These considerations permit us to propose a transposition of integrity criteria inspired by the Along the Track Protection Interval method [START_REF] Nikiforov | Integrity equations for safe train positioning using GNSS[END_REF] that are more adapted to the railway domain.

Reminders about Safety, attribute of RAMS

By definition, integrity is an attribute close to safety for a particular event occurring on a navigation solution. At this point, it is necessary to recall the safety definition as described in the standards (especially in EN 50126). Safety is defined by the "absence of unacceptable risk" linked to a probability of safety related failures, f S (t) (cf Appendice C in (CENELEC 2007)).

We considered as a "safety related failures" in a GNSS-based localisation system the following event: "the position error P E exceeds an alert limit AL, without being detected at the t i , time of the position i". This paper does not present the failure modes leading to this catastrophic event. However, a complete list of them for a GNSS/INS system can be found in [START_REF] Bhatti | Failure modes and models for integrated gps/ins systems[END_REF]. The probability of this event is linked to integrity risk mentioned in subsection 2.2. Only the considered interval of time is not the same. With the integrity definition, a temporal criterion, T T A has been introduced as a particular condition of the failure detection (T T A period of time).

In the IEC 61508-4 standard, Safety Integrity Level (SIL) is also defined as a quantified objective to reach related to the safety. Be careful, "Integrity" in SIL acronym does not have the same meaning of the integrity seen in this paper. The safety integrity concept is associated to the ability of a safety-related system satisfactorily performing the required safety functions in all specified conditions within a stated period of time. It is a general definition and the integrity point of view is considered for a specific safety related function: the localisation. A Tolerable Hazard Rate (THR) is also indicated in safety standards. In a safety-related systems evaluation context, these systems are often represented by configurations composed of channels. The GNSS/INS system performing the localisation function considered in the paper is a "1-out-of-1 with diagnosis" (1oo1D) voting where GNSS and INS represent a channel and the diagnostic is typically the monitoring integrity (cf figure 3).

EXTENDED INTEGRITY CONCEPT FOR A GNSS/INS LOCALISATION TRAIN SYSTEM

In order to give a justified confidence in the localisation system, the integrity concept must be defined, not only for the GNSS part, but also for the other parts, which constitute the localisation system (Here, the INS part). The extended integrity on a navigation solution provided by the GNSS/INS system depends on the integrity of the navigation solutions provided by the GNSS and the INS system parts. The integrity monitoring already exists in the hybridised systems (GNSS receiver and inertial system) but the main objective is only to evaluate the GNSS integrity. In this case, the other systems provide additional information to GNSS. In the paper, GNSS is coupled with an Inertial Navigation System. Inertial measurements permit us to ameliorate GNSS signal acquisition when the satellite signals are not available. In the integrity monitoring of GNSS/INS, the INS measurements' validity are not discussed and are supposed to internally guaranteed. However, the assumption about guaranteed validity is not conceivable in the case of a lack of compensation of the INS drifts. In consequence, this validity can be determined through the integrity for other systems other than GNSS. This leads to an extended integrity concept.

State representation of GNSS/INS localisation system

Before introducing the modification especially on P L computation, it is necessary to lay down some equations especially the considered state representation given by the Equation 2.

ẋ(t) = F (t)x(t) + G(t)w s (t) (1) z(t) = H(t)x(t) + w m (t)
with, F (t) is the (n × n) state matrix of system where the eigenvalues correspond to the poles of the system, x(t) is the (n × 1) state vector representing the system state at time t and its derivative ẋ(t), G(t) is the (n × 1) continuous system noise distribution matrix, w s (t) and w m (t) are, respectively, the system and measurement noise considered as white Gaussians with zero mean, z(t) is the (m × 1) measurement vector corresponding to measurements received by the sensors at time t, H(t) is the (m × n) measurement matrix determined by known properties of the system such as kinematics user and GNSS satellite geometry.

For the sake of brevity, the different matrix and vectors are not developed in this paper (In simulation, n = 17 and m ≥ 17 -it depends on the number of pseudoranges i.e. viewed satellites -). The evolution of the m -n difference is important to supervise for integrity monitoring especially in computation of the failure detecting threshold.

Extended integrity criteria determination for GNSS-based localisation systems

The important point here is that the error covariance matrix P of the system and the innovation vector δz are computed in a Kalman filter. An error covariance matrix is a measurement of the estimated accuracy of the estimated state. An innovation vector shows the difference between the observation of the system state and the estimated state (cf Equation 2). These two quantities are useful respectively for P L and minimum detectable bias, M DB determination.

δz(t) = z(t) -h(x(t)) (2)
The computation of P L related to GNSS measurements does not change. However, the contribution of other positioning systems must be taken into account in P L. In the hybridized systems, a positioning system like the inertial one makes corrections of GNSS measurements possible. In consequence, we do not see, in current P L, the error contribution of inertial measurements as they are considered fault-free. By extension, we do not fully evaluate the global integrity. Equation 3 represents the P L guaranteed by GNSS measurements (Le Marchand 2010).

P L GN SS = M AX(H(1, i) 2 + H(2, i) 2 ) × M DB GN SS (3)
with, M DB GN SS is the minimum detectable bias and H(1, i) and H(2, i) are linked to horizontal coordinates (x and y) of the i-th measurement given by the GNSS part only. For an inertial system, the P L computation is not immediate (cf Equation 4). Indeed, the information of the inertial system do not appear in H matrix because these data are used to correct GNSS measurements and then they are merged. To find inertial measurements, pseudoranges ρ and pseudorange rates δρ must be subtracted from the residuals. M DB IN S is a minimum detectable bias for INS part. These M DBs depend on the performances of bias detection methods i.e. the probability that a bias is detected on GNSS or INS part. In this paper, M DB GN SS = 5.49 (meters) based on M DB calculation in (Le Marchand 2010) and M DB IN S = 0.05 (meters) based on the performances of a Slowly Growing Errors detection [START_REF] Ochieng | User level integrity monitoring and quality control for high accuracy positioning using gps/ins measurements[END_REF].

P L IN S = M AX((δz(1, i) -ρ(1, i) -δρ(1, i)) 2 +(δz(2, i) -ρ(2, i) -δρ(2, i)) 2 ×M DB IN S (4)
Finally, a single measurement of P L must be considered. For safety reasons, the P L chosen is simply the maximum between P L GN SS and P L IN S (cf Equation 5).

P L System = M AX(P L GN SS , P L IN S ) (5)
It remains the time to alert (T T A) and integrity risk criteria for us to define. They are constant quantities and represent integrity requirements (recommended values given in the subsection 4.2.1).

INTEGRITY APPROACH FOR EVALUATING SAFETY

In subsection 2.3, some similarities exist between integrity and safety attributes according to the safety related failure. This failure is linked to the integrity risk criterion, defined in subsection 2.2. The T T A period of time is also defined in this part and can be seen as an occurrence condition of the event specified as catastrophic. That is why the integrity risk is considered as a conditional probability and the Bayes theorem can be used.

Demonstration and formalism of the link between integrity and safety

To determine the link between integrity and safety for a localisation system, let consider the following events:

• A t i , the event that U L(t i ) > AL with U L(t i ),
observable quantity reflecting the current position error P E at t i (cf subsection 2.2).

• B t i is the undetection at t i , an unsafe situation linked to the integrity risk.

• B [t i ;t i +T T A] is the undetection on [t i ; t i + T T A] interval, it is the B t i event associated with the T T A period of time

As we can see in subsection 2.3, the safety attribute for a localisation system can be characterised by the probability of "P E exceeds AL, without being detected at t i " event. In consequence, this probability can be expressed by Equation 6and, for a given mission of a T m duration in seconds and T e, a sample time, by Equation 7.

f S (t i ) = p(A t i |B t i ) (6)
P F H = 3600 T m × int( T m T e ) i=0 f S (t i ) (7

Te

Figure 2: T T A period of time seen as a floating windows T T A period of time must be taken into account (cf Equation 8).

p(IR(t i )) = p(A [t i ;t i +T T A] |B [t i ;t i +T T A] ) (8)
Thanks to Equation 6, f S (t i ) can be linked with p(IR(t i )) (cf Equation 9). T T A, t i and T e are expressed in seconds.

p(IR(t

i )) = i+int( T T A T e ) j=i f S (t j ) (9)
All along a mission, T T A period of time is seen as a floating windows (cf figure 2). It permits us to make explicit a limit case occurring at t n , previous time before end of mission T m. In consequence,

t n ≤ T m < t n+1 . f S (t i ) = 1 int( T T A T e ) × p(IR(t i )) (10)
Equation 10 permits to express the safety attribute (characterised by f S (t i )) throughout the integrity (characterised by p(IR(t i )) assuming that f S (t i ), ..., f S (t i + int( T T A T e )) are equiprobables. Finally, thanks to Equation 7, a P F H can be computed.

Safety quantitative analysis with integrity results

Thanks to the extended integrity concept and the link between safety and integrity, it is possible to make a quantitative analysis of safety. Before introducing these results, the experimental protocol must be posed.

Experimental protocol

In order to generate integrity results, a simulated GNSS/INS localisation system is necessary. Its architecture is given in the figure 3 taken from [START_REF] Groves | Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems[END_REF]. To evaluate its integrity performances, some requirements must be followed: alert limit, time to alarm and risk integrity. These requirements come from GNSS rail user forum [START_REF] Wiss | GNSS Rail User Forum: Requirements of Rail Applications[END_REF] and Center). In consequence, four seconds is chosen for T T A. Concerning AL, the safety margin (< 20 m) constitutes a justified quantification. This margin refers to the reliability of braking system in ETCS system (Note that it should not be mixed up with the reliability of ETCS system) [START_REF] Hougardy | Introduction to ETCS braking curves[END_REF].

At the level of simulation, the probabilities are difficult to obtain. It is preferable to choose the term of occurrence rate, but that is easily linked to probability. Consider an event, E, that occurs in some trials but not others. The probability of its occurrence in a trial is given by Equation 11.

P b(E) = lim n→+∞ n E n (11)
With, any E event, n E , the number of trials where event E occurs and n, all the possible trials. About the simulation details, it is Matlab R routine simulating a mission with a duration of 8059s with a sampling of 0.20s.

Safety quantitative results

In failure per hour. The following definitions seen in subsection 2.3, the SIL associated with this P F H is SIL1. However, SIL1 is not a sufficient level. The SIL required must be SIL4 to ensure the system is certified for safety-related railway applications. To achieve a SIL4 without compromising the availability, a 2oo3D (triple redundancy) or 2oo4D (quadruple redundancy) voting with diagnostic and map matching seems to be the best solutions for fault tolerance against dangerous detected or undetected failures [START_REF] Ding | A novel method for SIL verification based on system degradation using reliability block diagram[END_REF].

CONCLUSION

In this paper, a quantitative safety evaluation of a localisation system based on GNSS has been proposed through the determination of its integrity risk. In a context of railway infrastructure optimisation (reduction of balises), the question of inertial systems performances must be asked as well as GNSS ones. In consequence, the integrity, firstly defined as a GNSS performance attribute, must be extended to inertial systems. After formalising this extended integrity, protection levels and integrity risk are computed about a tightly coupled GNSS/INS architecture. The architecture envisaged today in a lot of railway projects. Thanks to a probabilistic demonstration, a link between safety and integrity is proposed permitting us to eliminate the constraint of safety evaluation for the GNSS based localisation system through the incurred integrity risk.

  Figure 1: Relations between xUL, xPL and xAL in different cases in railway domain

Figure 4 :

 4 Figure 4: Integrity requirements proposed in ERTMS point of view

)

  P F H represents the probability of dangerous failure per hour. This notation, coming from IEC 61508, permits us to allocate a Safety Integrity Level (cf 4.2.2) to the localisation system. The int(x) function for any number x is the integer part of x. Indeed, T m can not be necessary a multiple of T e (same observation with T T A in relation to T e in Equation6). Now, to express the probability linked to the integrity risk (further, this probability is written by p(IR)),

			TTA
	0	t 1	t 2	t 3	t 4		t n -int( Te TTA )	t n	Tm	t
	0	t 1	t 2	t 3	t 4	t 5	t n -int( Te TTA )	t n	Tm	t
	0						t n -int( Te TTA )	t n	Tm	t

Table 1 :

 1 Integrity recommendations for safety related railways applications[START_REF] Wiss | GNSS Rail User Forum: Requirements of Rail Applications[END_REF]

	Applications				Integrity requirements AL (m) TTA (s) Integrity risk (10 -x /h)
	ATC on high density lines / Station / Parallel track 2.5	<1	1 × 10 -7
	Train Control on medium density lines	20	<1	1 × 10 -7
	Train Control on low density lines		50	<1	4.8 × 10 -6
	GNSS		INS
	GNSS receiver	IMU
	Pseudoranges computation	Resolution of inertial navigation equations
		GNSS/INS integration
		by Extended
		Kalman Filter	Software
			Integrated navigation
	solution Integrity monitoring	Hardware
	IR required AL			data flux
	TTA	PL/UL	IR achieved	Open-loop
				correction
	They are not accepted requirements in railways just
	recommended). However, these values of require-
	ments seem to be difficult to reach, particularly, time
	to alert, which appear very restrictive. In this paper,
	other values are proposed to be compatible with
	the European Railway Traffic Management System
	(ERTMS). According to the ERTMS performance
	requirements, T T A criteria such as those considered
	in Table 1 should be reviewed and sized in realistic
	way. Figures 4a and 4b propose to quantify the AL
	and T T A depending on emission ERTMS/ETCS
	messages delays from train to RBC (Radio Block

Requirements

Figure 3: Architecture of GNSS/INS navigation system and integrity monitoring summarise in table 1 (Note that, for the moment.

  Table 2, the results useful for p(IR) and, thank to the link described before, for safety (throughout f S (t i ) are introduced. The different quantities in the table are means values during the simulation. According to this Table, the GNSS/INS localisation system reaches the integrity risk of 9.7214 × 10 -5 . Thanks to Equations 7 and 10, the f

S (t) and P F H deduced from p(IR) are estimated respectively at 4.8607 × 10 -6 (dimensionless quantity) and 2.1713 × 10 -6

Table 2 :

 2 Integrity and safety quantitative evaluation

	System part	Integrity results PL(part) PL(system) p(IR)	f S (t)	P F H
	GNSS INS	15.60 6.61	15.60	9.7214 × 10 -5 4.8607 × 10 -6 2.1713 × 10 -6

ACKNOWLEDGEMENTS

The authors would like to thank the Technological Research Institute Railenium and the IFSTTAR for their support in this work and the International Campus on Safety and Intermodality in Transportation (CISIT) and the European Commission via the European Regional Development Fund (ERDF) for making possible this work. The acknowledgements go also to the GEOLOC team from IFSTTAR-Nantes for GNSS and INS data and their help in using them in the present work.