Adel Blouza
email: adel.blouza@univ-rouen.fr

Linda El Alaoui
email: elalaoui@math.univ-paris13.fr

A POSTERIORI ANALYSIS OF A RICHARDS PROBLEM WITH A FREE BOUNDARY

Richards equation models the water flow in a partially saturated underground porous medium. When it rains on the surface, boundary conditions of Signorini type may be considered on this part of the boundary. We recall the well-posedness of this problem together with its discretization by an implicit Euler's scheme in time and finite elements in space. We establish quasi-optimal a posteriori estimates for the solution of the discrete problem. Some numerical experiments confirm the interest of this discretization.

Introduction

The following ∂ t Θ(ψ) -∇ • K w Θ(ψ) ∇(ψ + z) = 0, (1.1) models the flow of a wetting fluid, mainly water, in the underground medium located just under the surface, hence in an unsaturated medium, see L.A. Richards [START_REF]Richards -Capillary conduction of liquids through porous mediums[END_REF] for the introduction of this type of models. In opposite to Darcy's or Brinkman's systems (see [START_REF]On a hierarchy of approximate models for flows of incompressible fluids through porous solid[END_REF] for all these models), this equation which is derived by combining Darcy's generalized equation with the mass conservation law is highly nonlinear: This follows from the fact that, due to the presence of air above the surface, the porous medium is only partially saturated with water. The unknown ψ is the difference between the pressure of water and the atmospherical pressure.

This equation is usually provided with Dirichlet or Neuman type boundary conditions. Indeed, Neumann boundary conditions on the underground part of the boundary are linked to the draining of water outside of the domain and Dirichlet boundary conditions on the surface are introduced to take into account the rain. However, when it is raining cats and dogs, the part of the domain saturated with water reaches the boundary. Then the Dirichlet boundary conditions must be replaced by variational inequalities of the following type:

-ψ ≥ 0, v(ψ) • n ≥ v r • n, ψ v(ψ) • n -v r • n = 0, (1.2)
where v(ψ) is the flux v(ψ) = -K w Θ(ψ) ∇(ψ + z).

(1.3)

Here, n stands for the unit normal vector to the surface and v r is a given rain fall rate. We refer to the thesis of H. Berninger [START_REF] Berninger | Domain decomposition methods for elliptic problems with jumping nonlinearities and application to the Richards equation[END_REF] for the full derivation of this model from hydrology laws. It is not so easy to give a mathematical sense to the system (1.1) -(1.2). As standard, the key argument for the analysis of problem (1.1) is to use Kirchhoff's change of unknowns. Indeed, after this transformation, the new equation fits the general framework proposed in [START_REF] Alt | Luckhaus -Quasilinear elliptic-parabolic differential equations[END_REF]. On the other hand, the inequality in (1.2) can be handled by a variational formulation. In [START_REF] Bernardi | El Alaoui -The rain on underground porous media, Part I. Analysis of the model[END_REF], we prove that the corresponding variational problem is well-posed when the coefficients and the data are smooth enough, but without any restriction on their maximal value.

The discretization of problem (1.1) has been studied in many papers when provided with standard boundary conditions, see [START_REF] Abdellatif | Touihri --A priori error analysis of an Euler implicit, spectral discretization of Richards equation[END_REF], [START_REF] Gabbouhy | Analyse mathématique et simulation numérique des phénomènes d'écoulement et de transport en milieux poreux non saturés. Application à la région du Gharb[END_REF], [START_REF] Radu | Knabner -Order of convergence estimates for an Euler implicit, mixed finite element discretization of Richards' equation[END_REF], [START_REF] Schneid | Radu -A priori error estimates for a mixed finite element discretization of the Richards' equation[END_REF], [START_REF] Sochala | Piperno -Mass conservative BDF-discontinuous Galerkin/explicit finite volume schemes for coupling subsurface and overland flows[END_REF] and [START_REF] Woodward | Analysis of expanded mixed finite element methods for a nonlinear parabolic equation modeling flow into variably saturated porous media[END_REF], however it seems that it has not been treated when provided with the boundary inequality (1.2). As first proposed in [START_REF] Bernardi | El Alaoui -The rain on underground porous media, Part I. Analysis of the model[END_REF]Section 3, we write here a discretization of system (1.1) -(1.2), in two steps: (i) We first use the Euler's implicit scheme to build a time semi-discrete problem, where one of the nonlinear terms is treated in an explicit way for simplicity; (ii) We then construct a fully discrete problem that relies on the Galerkin method and finite elements in the spatial domain. We also recall the convergence properties of this discretization. We observe that the discrete problem is unconditionally stable and that its convergence does not require any limitation on the time step and mesh size. Indeed, these two properties are very important for performing adaptivity.

The aim of this paper is to prove a posteriori estimates for the discrete problem. We refer to [START_REF] Hild | Nicaise -A posteriori error estimations of residual type for Signorini's problem[END_REF] for the first a posteriori analysis of a problem with Signorini type boundary conditions but in the stationary case. Following the approach in [START_REF] Bergam | Mghazli -A posteriori analysis of the finite element discretization of some parabolic equations[END_REF], in this a posteriori analysis, we try to uncouple the errors issued from the time and space discretizations as much as possible. This leads to a simple adaptivity strategy, as first described in [START_REF] Bernardi | Süli -Time and space adaptivity for the second-order wave equation[END_REF]. Some numerical experiments confirm the interest of the discretization.

The outline of the paper is as follows.

• In Section 2, we present the variational formulation of the full system and recall its well-posedness.

• Section 3 is devoted to the description of the time semi-discrete problem and of the fully discrete problem. There also, we recall their well-posedness and their convergence.

• The a posteriori analyses of the time semi-discretization and space discretization are performed in Sections 4 and 5, respectively.

• In Section 6, we present a few numerical experiments.

The continuous problem and its well-posedness

Let Ω be a bounded connected open set in R d , d = 2 or 3, with a Lipschitz-continuous boundary ∂Ω. From now on and for simplicity, we assume that Ω is a polygon (d = 2) or a polyhedron (d = 3), and we denote by n the unit outward normal vector to Ω on ∂Ω. We assume that ∂Ω admits a partition without overlap into three parts Γ B , Γ F , and Γ G (these indices mean "bottom", "flux" and "ground", respectively) and that Γ B has a positive measure. We also assume that this partition is smooth enough, i.e. ∂Γ B , ∂Γ F and ∂Γ G are Lipschitz submanifolds of ∂Ω.

In order to perform the Kirchhoff's change of unknowns in problem (1.1), we observe that, since the conductivity coefficient K w is positive, the mapping:

x → K(x) = x 0 K w Θ(ξ) dξ,
is one-to-one from R into itself. Thus, by setting u = K(ψ), we obtain an equivalent equation to (1.1). Moreover, the Kirchhoff's change of unknowns has the further property of preserving the positivity, so that (1.2) can be written in terms of u. All this leads to the following system

α ∂ t u + ∂ t b(u) -∇ • ∇u + k • b(u)e z = 0 in Ω×]0, T [ (2.1) u = u B on Γ B ×]0, T [ (2.2) -∇u + k • b(u)e z • n = f F on Γ F ×]0, T [ (2.3) u ≤ 0, -∇u + k • b(u)e z • n ≥ q r • n, u ∇u + k • b(u)e z + q r • n = 0 on Γ G ×]0, T [, (2.4
)

u| t=0 = u 0 inΩ, (2.5)
where -e z stands for the unit vector in the direction of gravity. The unknown is now the quantity u. The data are the Dirichlet boundary condition u B on Γ B and the initial condition u 0 on Ω, together with the boundary conditions f F and q r on the normal component of the flux, with f F corresponding to the draining of water and q r corresponding to the rain. Finally, the coefficients b and k are supposed to be nonnegative, while α is a positive constant. From now on, we assume that the function b is continuously differentiable on R, with a nonnegative, bounded and Lipschitzcontinuous first derivative b , and that the function k • b is continuous, bounded, and uniformly Lipschitz-continuous on R. Remark It must be observed that the unknown u has no physical meaning. However, the pressure ψ can easily be recovered from u, for instance by a dichotomy algorithm, see the definition of the operator K.

In what follows, we use the whole scale of Sobolev spaces W m,p (Ω), with m ≥ 0 and 1 ≤ p ≤ +∞, equipped with the norm • W m,p (Ω) and seminorm | • | W m,p (Ω) , with the usual notation H m (Ω) when p = 2. For any separable Banach space E equipped with the norm • E , we denote by C 0 (0, T ; E) the space of continuous functions on [0, T ] with values in E. For each integer m ≥ 0, we also introduce the space H m (0, T ; E) as the space of measurable functions on ]0, T [ with values in E such that the mappings: v → ∂ t v E , 0 ≤ ≤ m, are square-integrable on ]0, T [.

To write a variational formulation for the problem, we introduce the time-dependent subset

V(t) = v ∈ H 1 (Ω); v| Γ B = u B (•, t) and v| Γ G ≤ 0 .
It is readily checked that each V(t) is closed and convex, see [START_REF] Berninger | Domain decomposition methods for elliptic problems with jumping nonlinearities and application to the Richards equation[END_REF]Prop. 1.5.5], when u B belongs to C 0 (0, T ; H 1 2 (Γ B )). Thus, we are led to consider the following variational problem (with obvious notation for L 2 (0, T ; V) )

Find u in L 2 (0, T ; V) such that

u| t=0 = u 0 ,
and that, for a.e.

t in ]0, T [, ∀v ∈ V(t), α Ω (∂ t u)(x, t)(v -u)(x, t) dx + Ω ∂ t b(u) (x, t)(v -u)(x, t) dx + Ω ∇u + k • b(u)e z (x, t) • ∇(v -u) (x, t) dx ≥ - Γ F f F (τ, t)(v -u)(τ, t) dτ - Γ G (q r • n)(τ, t)(v -u)(τ, t) dτ.
We refer to [4, Prop. 2.2] for the proof of the following statement. Proving that problem (2.3) -(2.4) is well-posed is not at all obvious. Indeed, standard arguments yield the uniqueness of the solution, and its existence only when a rather restrictive condition on the function b is satisfied. In any case, we refer to [4, Prop. 2.3 and Thm 4.3] for the following result. For brevity, we set:

X = L 2 (0, T ; V) ∩ H 1 (0, T ; L 2 (Ω)).
Theorem 2.1. For any data u B , f F , q r and u 0 satisfying

u B ∈ H 1 (0, T ; H 1 2 00 (Γ B )), f F ∈ C 0 (0, T ; H 1 2 (Γ F )), q r ∈ C 0 (0, T ; H 1 2 (Γ G ) d ), u 0 ∈ H 1 (Ω),
together with the compatibility conditions

u 0 (x) = u B (x,

The discrete problems

We present successively the time semi-discrete problem constructed from the backward Euler's scheme, next a finite element discretization of this problem relying on standard, conforming, finite element spaces.

3.1.

The time semi-discrete problem. Since we intend to work with non uniform time steps, we introduce a partition of the interval [0, T ] into subintervals [t n-1 , t n ], 1 ≤ n ≤ N , such that 0 = t 0 < t 1 < • • • < t N = T . We denote by τ n the time step t n -t n-1 , by τ the N -tuple (τ 1 , . . . , τ N ) and by |τ | the maximum of the τ n , 1 ≤ n ≤ N .

As already hinted in Section 1, the time discretization mainly relies on a backward Euler's scheme, where the nonlinear term k • b(u) is treated in an explicit way for simplicity. Thus, the semi-discrete problem reads

Find (u n ) 0≤n≤N in N n=0 V(t n ) such that u 0 = u 0 in Ω, (3.1)
and, for 1

≤ n ≤ N , ∀v ∈ V(t n ), α Ω u n -u n-1 τ n (x)(v -u n )(x) dx + Ω b(u n ) -b(u n-1 ) τ n (x)(v -u n )(x) dx + Ω ∇u n + k • b(u n-1 )e z (x) • ∇(v -u n )(x) dx (3.2) ≥ - Γ F f F (τ, t n )(v -u n )(τ ) dτ - Γ G (q r •n)(τ, t n )(v -u n )(τ ) dτ.
It can be noted that this problem makes sense when both f F and q r are continuous in time. On the other hand, it is readily checked that the problem for each n is equivalent to the minimization of a convex functional on the convex set V(t n ), whence the following result.

Proposition 3.1. For any data u B , f F , q r and u 0 satisfying

u B ∈ H 1 (0, T ; H 1 2 00 (Γ B )), f F ∈ C 0 (0, T ; L 2 (Γ F )), q r ∈ C 0 (0, T ; L 2 (Γ G ) d ), u 0 ∈ H 1 (Ω), (3.3)
and (2.7), for any nonnegative coefficient α, problem (3.1)-(3.2) has a unique solution in N n=0 V(t n ).

3.2. The fully discrete problem. We recall that Ω is a polygon (d = 2) or a polyhedron (d = 3). For 0 ≤ n ≤ N , let (T n h ) h be a regular family of triangulations of Ω (by triangles or tetrahedra), in the sense that, for each h:

• Ω is the union of all elements of T n h ; • The intersection of two different elements of T n h , if not empty, is a vertex or a whole edge or a whole face of both of them; • The ratio of the diameter h K of any element K of T n h to the diameter of its inscribed circle or sphere is smaller than a constant σ independent of h and n. As usual, h n stands for the maximum of the diameters h K , K ∈ T n h . We make the further and non restrictive assumption that, for all n, Γ B , Γ F and Γ G are the union of whole edges (d = 2) or whole faces (d = 3) of elements of T n h . We now introduce the finite element space

V n h = v h ∈ H 1 (Ω); ∀K ∈ T n h , v h | K ∈ P 1 (K)
, where P 1 (K) is the space of restrictions to K of affine functions on R d . Let I n h denote the Lagrange interpolation operator at all the vertices of elements of T n h with values in V n h and i nB h the corresponding interpolation operator on Γ B . Assuming that u B is continuous where needed, we then define for each n, 0 ≤ n ≤ N , the subset of V n h :

V n h = v h ∈ V n h ; v h | Γ B = i nB h u B (•, t n ) and v h | Γ G ≤ 0 .
We are thus in a position to write the discrete problem, constructed from problem (3.1)-(3.2) by the Galerkin method,

Find (u n h ) 0≤n≤N in N n=0 V n h such that u 0 h = I 0 h u 0 in Ω, (3.4)
and, for 1

≤ n ≤ N , ∀v h ∈ V n h α Ω u n h -I n h u n-1 h τ n (x)(v h -u n h )(x) dx + Ω b(u n h ) -I n h b(u n-1 h ) τ n (x)(v h -u n h )(x) dx (3.5) + Ω ∇u n h + I n h k • b(u n-1 h )e z (x) • ∇(v h -u n h )x) dx ≥ - Γ F f F (τ, t n )(v h -u n h )(τ ) dτ (3.6) - Γ G (q r •n)(τ, t n )(v h -u n h )(τ ) dτ. (3.7)
We also recall the theorem below which is proved in [START_REF] Bernardi | El Alaoui -The rain on underground porous media, Part I. Analysis of the model[END_REF].

Theorem 3.1. For any data u B , f F , q r and u 0 satisfying 3.3, 2.6 and

u B ∈ C 0 (Γ B × [0, T ]), u 0 ∈ C 0 (Ω),
for any nonnegative coefficient α, problem (3.4)-(3.5) has a unique solution.

3.3.

Convergence of the discretization. We introduce a lifting u * B of an extension of u B to ∂Ω which belongs to H 1 (0, T ; V) and satisfies for a.e. x ∈ Ω, u * B (x, 0) = u 0 (x), (this requires assumption ??) together with the stability property

u * B H 1 (0,T ;H 1 (Ω)) ≤ c u B H 1 (0,T ;H 1 2 00 (Γ B ))
.

We thus write the problem satisfied by the function Lemma 3.1. For any data u B , f F , q r and u 0 satisfying 2.6 and

u * = u -u * B . Similarly, when setting u * n h = u n h -I n h u * B (t n ) (
u B ∈ H 1 (0, T ; H s (Γ B )), f F ∈ C 0 (0, T ; H 1 2 (Γ F )), q r ∈ C 0 (0, T ; H 1 2 (Γ G ) d ), u 0 ∈ H s+ 1 2 (Ω), (3.8)
for some s > d-1 2 , the sequence (u * n h ) 0≤n≤N satisfies the following bound, for

1 ≤ n ≤ N , α n m=1 τ m u * m h -u * m-1 h τ m 2 L 2 (Ω) + |u * n h | 2 H 1 (Ω) ≤ c 1 + I n h u * B 2 H 1 (0,T ;H 1 (Ω)) + f F 2 C 0 (0,T ;H 1 2 (Γ F )) + q r 2 C 0 (0,T ;H 1 2 (Γ G ) d )
.

The following result is now easily derived from this boundedness. We denote by u τ and u hτ the functions which are affine on each interval

[t n-1 , t n ], 1 ≤ n ≤ N ,

A posteriori analysis of the time semi-discretization

As already hinted, we wish to uncouple as much as possible the time and space errors. Equivalently, we use the triangle inequality

u -u hτ X ≤ u -u τ X + u τ -u hτ X , (4.1)
and we evaluate successively the two terms in the right-hand side. In this section, we are interested in the evaluation of the first term: We propose error indicators and establish their reliability (i.e., that they provide an upper bound for the error). Their efficiency (i.e. that they are bounded from above by the error) is only hinted.

To go further, we introduce a lifting u * B of an extension of u B to ∂Ω which belongs to H 1 (0, T ; V) and satisfies

u * B H 1 (0,T ;H 1 (Ω)) ≤ c u B H 1 (0,T ;H 1 2 00 (Γ B ))
.

Then it is readily checked that the function u * = u -u * B is a solution of a problem very similar to 2.1, with a modified function b and an additional term in the right-hand side, namely the quantity

-Ω F B (x, t)(v -u)(x) dx defined by Ω F B (x, t)v(x) dx = α Ω (∂ t u * B )(x, t)v(x) dx + Ω (∇u * B )(x, t) • (∇v)(x) dx.
An analogous property holds for the semi-discrete problem (3.1)-(3.2). So, from now on and for simplicity, we assume that u B is equal to zero. In a first step, in order to evaluate the norms of u -u τ , we compute the quantity, for a.e. t in [0, T ]

R(t) = α Ω ∂ t (u -u τ ) (x, t)(u -u τ )(x, t) dx + Ω ∂ t b(u) -∂ t b(u τ ) (x, t)(u -u τ )(x, t) dx + Ω ∇(u -u τ ) (x, t) • ∇(u -u τ ) (x, t) dx + Ω k • b(u) -k • b(u τ ) (x, t)e z • ∇(u -u τ ) (x, t) dx.
The interest of this quantity appears in the foliowing lemma. We refer to [4, Prop. 2.3] for the arguments for proving it. Lemma 4.1. Assume that the function b belongs to W 2,∞ (R). The quantity R satisfies for all t in [0, T ]

t 0 R(s) ds ≥ α 2 (u -u τ )(t) 2 L 2 (Ω) + 1 2 t 0 |(u -u τ )(s)| 2 H 1 (Ω) ds -c D t 0 (u -u τ )(s) 2 L 2 (Ω) ds,
where the constant c D only depends on the data.

Proof. Obviously, the first and third terms in R(t) lead to (note that u -u τ vanishes at time t = 0)

t 0 α Ω ∂ t (u -u τ ) (x, s)(u -u τ )(x, s) dx + Ω ∇(u -u τ ) (x, s) • ∇(u -u τ ) (x, s) dx ds = α 2 (u -u τ )(t) 2 L 2 (Ω) + t 0 |(u -u τ )(s)| 2 H 1 (Ω) .
So, we now investigate successively the two remaining terms. 1) We use the decomposition

Ω ∂ t b(u) -∂ t b(u τ ) (x, t)(u -u τ )(x, t) dx = Ω b (u)(x, t) ∂ t (u -u τ ) (x, t)(u -u τ )(x, t) dx + Ω b (u) -b (u τ ) (x, t)(∂ t u τ )(x, t)(u -u τ )(x, t) dx,
and integrate the first term by parts with respect to t, which gives

t 0 Ω ∂ t b(u) -∂ t b(u τ ) (x, s)(u -u τ )(x, s) dx ds = Ω b (u)(x, t) 2 (u -u τ ) 2 (x, t) dx - 1 2 t 0 Ω b (u)(x, s)(∂ t u)(x, s)(u -u τ ) 2 (x, s) dx ds + t 0 Ω b (u) -b (u τ ) (x, s)(∂ t u τ )(x, s)(u -u τ )(x, s) dx,
Next, the nonnegativity of b , the boundedness of b and the Lipschitz-continuity of b yields

t 0 Ω ∂ t b(u) -∂ t b(u τ ) (x, s)(u -u τ )(x, s) dx ds ≥ -c(u, u τ ) t 0 (u -u τ )(•, s) 4 L 4 (Ω) ds,
where c(u, u τ ) > 0 depends on the norm of u and u τ in X. However, using a further decomposition allows us to replace it by a constant only depending on u X , and it follows from Remark ?? that this norm is bounded as a function of the data. Next, we use an interpolation inequality (see [13, Chap. 1, Prop. 2.3]) and the Poincaré-Friedrichs inequality: for any v in H 1 (Ω) vanishing on Γ B ,

v L 4 (Ω) ≤ v 1-d 4 L 2 (Ω) (c|v| H 1 (Ω) ) d 4 ≤ c (1 - d 4 ) v L 2 (Ω) + c d 4 |v| H 1 (Ω) ,
and conclude by applying it with v = u -u τ and using a Young's inequality.

2) Finally, to bound the last term, we combine the Lipschitz-continuity of k • b together with a Young's inequality

t 0 Ω k • b(u) -k • b(u τ ) (x, s)e z • ∇(u -u τ ) (x, s) dx ds ≥ - 1 4 t 0 |(u -u τ )(•, s)| 2 H 1 (Ω) ds -c t 0 (u -u τ )(•, s) 2 L 2 (Ω) ds .
Combining all this yields the desired estimate. We introduce the subset V 0 :

V 0 = v ∈ H 1 (Ω); v| Γ B = 0 and v| Γ G ≤ 0 . (4.2)
,1

Thus, applying problem (2.1) with v = u τ which now belongs to V 0 , we derive

R(t) ≤ -α Ω (∂ t u τ )(x, t)(u -u τ )(x, t) dx - Ω ∂ t b(u τ ) (x, t)(u -u τ )(x, t) dx Ω (∇u τ )(x, t) • ∇(u -u τ ) (x, t) dx (4.3) - Ω k • b(u τ ) e z (x, t) • ∇(u -u τ ) (x, t) dx - Γ F f F (τ, t)(u -u τ )(τ, t) dτ (4.4) - Γ G (q r •n)(τ, t)(u -u τ )(τ, t) dτ. (4.5)
On the other hand, we define the space

W = v ∈ H 1 (Ω); v| Γ B = 0 and v| Γ G = 0 . (4.6)
Using problem (3.2) with v = u n ± w (which also belongs to V 0 for any w in W), we obtain the equation for all w ∈ W,

α Ω u n -u n-1 τ n (x)w(x) dx + Ω b(u n ) -b(u n-1 ) τ n (x)w(x) dx (4.7)
+

Ω ∇u n + k • b(u n-1 )e z (x) • ∇w(x) dx = - Γ F f F (τ, t n )w(τ ) dτ, (4.8)
or, equivalently, by taking w successively in D(Ω) and in

D(Ω ∪ Γ F ), α u n -u n-1 τ n + b(u n ) -b(u n-1 ) τ n -∇ • ∇u n + k • b(u n-1 )e z = 0 in Ω, (4.9)
-∇u n + k • b(u n-1 )e z • n = f F (•, t n ) on Γ F . (4.10)
Inserting these equations into (3.2) also yields

∀v ∈ V 0 , Γ G ∇u n + k • b(u n-1 )e z • n(τ )(v -u n )(τ ) dτ ≤ - Γ G (q r •n)(τ, t n )(v -u n )(τ ) dτ. (4.11)
For any function f continuous on [0, T ], we agree to denote by π τ f the piecewise constant function equal to f (t n ) on each interval ]t n-1 , t n ], 1 ≤ n ≤ N . Multiplying both equations in (4.9) by u -u τ and adding them to (4.3) (note that ∂ t u τ is equal to u n -u n-1 τn ) and also using (4.11)

with v = u gives on the interval ]t n-1 , t n [ R(t) ≤ 4 j=1 R j (t) + D(t), with R 1 (t) = Ω b(u n ) -b(u n-1 ) τ n -∂ t b(u τ ) (x, t)(u -u τ )(x, t) dx, R 2 (t) = Ω ∇(u n -u τ )(x, t) • ∇(u -u τ )(x, t) dx, R 3 (t) = Ω k • b(u n-1 ) -k • b(u τ ) (x, t)e z • ∇(u -u τ )(x, t) dx, R 4 (t) = - Γ G ∇u n + k • b(u n-1 )e z + π τ q r •n(τ, t)(u n -u τ )(τ, t) dτ,
while the data depending term D is given by

D(t) = - Γ F (f F -π τ f F )(τ, t)(u -u τ )(τ, t) dτ - Γ G (q r -π τ q r )•n(τ, t)(u -u τ )(τ, t) dτ.
It also follows from (4.11) that the quantity R 4 (t) is nonpositive. Despite its complexity, the first a posteriori bound for the error can easily be derived from the last inequality. Indeed, we observe that, for all t in [t n-1 , t n ],

u τ = u n - t n -t τ n (u n -u n-1 ).
This leads to the following definition of the time error indicators:

η n τ = τ 1 2 n u n h -u n-1 h H 1 (Ω) + τ 1 2 n b(u n h ) -b(u n-1 h ) τ n -B n h u n h -u n-1 h τ n L 2 (Ω) ,
where the function B n h is defined as, for 1 ≤ n ≤ N ,

B n h = b (u n-1 h ) + b (u n h ) 2 .
Note that the first term in the indicator is fully standard for a parabolic equation, see [3, equation (3.1)] or [START_REF] Johnson | Thomée -An a posteriori error estimate and adaptive timestep control for a backward Euler discretization of a parabolic problem[END_REF], while the second one comes from the nonlinear term ∂ t b(u) and is the same as introduced in [5, equation (4.1)].

Theorem 4.1. Assume that the function b belongs to W 2,∞ (R). The following a posteriori estimate holds between the solution u of problem (2.1) and the function u τ associated with the solution

(u n ) 0≤n≤N of problem (3.1)-(3.2). sup 0≤t≤T (u -u τ )(t) L 2 (Ω) + u -u τ L 2 (0,T ;H 1 (Ω) ≤ c N n=1 (η n τ ) 2 1 2 + u τ -u hτ X + ε (D) + ε (b) ,
where the quantities ε (D) and ε (b) are given by (with obvious definition of u * Bτ )

ε (D) = f F -π τ f F L 2 (Γ F ×]0,T [) + q r -π τ q r L 2 (Γ G ×]0,T [) d + u * B -u * Bτ X ε (b) = N n=1 (ε (b)n ) 2 1 2 , with ε (b)n = b (u hτ ) -B n h u n h -u n-1 h τ n L 2 (t n-1 ,tn;L 2 (Ω)) .
Proof. Using Cauchy-Schwarz and Young inequalities on each R j (t) and D(t) together with a triangle inequality to replace u n by u n h and u n-1 by u n-1 h leads to

t 0 R(s) ds ≤ c N n=1 (η n τ ) 2 + u τ -u hτ 2 X + (ε (D) ) 2 + N n=1 (ε (b)n 2 + c t 0 (u -u τ )(s) 2 L 2 (Ω) ds.
We then apply Lemma 4.1 and use the Gronwall's lemma, see [10, chap. V, lemma 1.8] for instance, to derive the desired bound (the further term including u * B comes from standard arguments). Remark More complex arguments would also allow us to derive a bound for the term u -u τ H 1 (0,T ;L 2 (Ω)) , we prefer to skip them for brevity.

Proving the efficiency of these indicators, i.e. deriving an upper bound of them as a function of the error apparently requires much more technical arguments than previously, so that we prefer to skip them. Since this estimate holds (and is not simple to establish) for the heat equation [START_REF] Bergam | Mghazli -A posteriori analysis of the finite element discretization of some parabolic equations[END_REF]Prop. 3.6] and also for the Richards's equations when Γ G = ∅, see [START_REF] Bernardi | Mghazli -A posteriori analysis of a space and time discretization of a nonlinear model for the flow in partially saturated porous media[END_REF]Prop. 4.10], both of them being provided with the same time discretization, we only say that it is likely.

A posteriori analysis of the space discretization

We now evaluate the second term of (4.1), following the same steps as in Section 4: We exhibit error indicators and prove first their reliability, second their efficiency. 5.1. The indicators and their reliability. In order to write the residual equation applied to u τ -u hτ , we compute the quantity, for 1

≤ n ≤ N , R n h =α Ω (u n -u n h ) -(u n-1 -I n h u n-1 h ) τ n (x)(u n -u n h )(x) dx + Ω (b(u n ) -b(u n h )) -(b(u n-1 ) -I n h b(u n-1 h )) τ n (x)(u n -u n h )(x) dx + Ω ∇(u n -u n h ) (x) • ∇(u n -u n h ) (x) dx + Ω k • b(u n-1 ) -I n h k • b(u n-1 h ) (x)e z • ∇(u n -u n h ) (x) dx.
We now prove the analogue of Lemma 4.1 for this quantity.

τ m R m h ≥ α 2 u n -u n h 2 L 2 (Ω) + 1 2 n m=1 τ m |u m -u m h | 2 H 1 (Ω) -c n-1 m=1 τ m u m -u m h 2 L 2 (Ω) - n m=1 (ε (I)m ) 2 ,
where the quantity ε (I)n is defined by

ε (I)n = τ -1 2 n (Id -I n h )(u n h -u n-1 h ) L 2 (Ω) + τ -1 2 n (Id -I n h )b(u n-1 h ) L 2 (Ω) + τ n (Id -I n h )k • b(u n-1 h ) L 2 (Ω) . (5.1)
Proof. For brevity, we denote by R 1 , . . . and R 4 the four lines in R n h . As for Lemma 4.1, the first and third terms in R n h can be handled obviously

τ n R 1 ≥ α 2 u n -u n h 2 L 2 (Ω) -u n-1 -u n-1 h 2 L 2 (Ω) + (u n -u n h ) -(u n-1 -u n-1 h ) 2 L 2 (Ω) -α (Id -I n h )(u n h -u n-1 h ) L 2 (Ω) u n -u n h L 2 (Ω) , τ n R 3 = τ n |u n -u n h | 2 H 1 (Ω) .
To estimate the second term in R h (n), we use the Lipschitz-continuity of b and the fact that it is not decreasing: if κ denotes its Lipschitz constant,

τ n R 2 ≥ κ 2 u n -u n h 2 L 2 (Ω) -u n-1 -u n-1 h 2 L 2 (Ω) -c (Id -I n h )b(u n-1 h ) 2 L 2 (Ω) .
Finally evaluating the last term relies on the Lipschitz property of k • b:

τ n R 4 ≥ - τ n 2 |u n -u n h | 2 H 1 (Ω) -c τ n u n-1 -u n-1 h 2 L 2 (Ω) .
We conclude by summing these 4 inequalities and then summing up on n.

As in Section 4, in order to avoid the technicality linked to the lifting of u B , we assume that u B and thus its interpolate are zero. It follows from problem

(3.1)-(3.2) that each R n h satisfies R n h ≤ -α Ω u n h -I n h u n-1 h τ n (x)(u n -u n h )(x) dx Ω b(u n h ) -I n h b(u n-1 h ) τ n (x)(u n -u n h )(x) dx - Ω (∇u n h )(x) • ∇(u n -u n h ) (x) dx - Ω I n h k • b(u n-1 h ) (x)e z • ∇(u n -u n h ) (x) dx - Γ F f F (τ, t n )(u n -u n h )(τ ) dτ - Γ G (q r • n)(τ, t n )(u n -u n h )(τ ) dτ.
For the space W defined in (4.6), we set

W h = V n h ∩ W. Indeed, taking v h = u n h ± w h in problem (3.5) yields ∀w h ∈ W h , α Ω u n h -I n h u n-1 h τ n (x)w h (x) dx + Ω b(u n h ) -I n h b(u n-1 h ) τ n (x)w h (x) dx + Ω ∇u n h + I n h k • b(u n-1 h )e z (x) • ∇w h (x) dx = - Γ F f F (τ, t n )w h (τ ) dτ.
We use this equation by taking w h equal to v h -u n h , where v h is any element in V n h which coincides with u n h on Γ G . This yields

R n h ≤ -α Ω u n h -I n h u n-1 h τ n (x)(u n -v h )(x) dx - Ω b(u n h )) -I n h b(u n-1 h )) τ n (x)(u n -v h )(x) dx - Ω (∇u n h )(x) • ∇(u n -v h ) (x) dx - Ω I n h k • b(u n-1 h ) (x)e z • ∇(u n -v h ) (x) dx - Γ F f F (τ, t n )(u n -v h )(τ ) dτ - Γ G (q r • n)(τ, t n )(u n -v h )(τ ) dτ.
We need some further notation: For each K in T n h , let E n K , E n F K and E n GK denote the sets of edges (d = 2) or faces (d = 3) of K that are not contained in ∂Ω or lie in Γ F or lie in Γ G , respectively. For any edge or face e in any of these sets, h e stands for the diameter of e. We also need the notation, for any quantity a,

[a] + = max{a, 0}, [a] -= min{a, 0}.
As standard, we introduce an approximation

f n F h of f n F = f F (•, t n )
in the space of piecewise constant functions relying on the triangulation of Γ F by edges or faces of elements of T n h and an approximation q n rh of q r (•, t n ) in the space of piecewise constant functions on Γ G . We are thus in a position to define the error indicators. Since we prefer to treat separately the boundary term linked to the inequation, we define: (i) For each K in T n h ,

η n K = h K α u n h -I n h u n-1 h τ n + b(u n h ) -I n h b(u n-1 h ) τ n -∇ • ∇u n h + I n h k • b(u n-1 h )e z L 2 (K) (5.2) + e∈E n K h 1 2 e [ ∇u n h + I h s n k • b(u n-1 h )e z • n] e L 2 (e) + e∈E n F K h 1 2 e ∇u n h + I n h k • b(u n-1 h )e z • n + f n F h L 2 (e) , (5.3)
where [•] e stands for the jump through e (making its sign precise is not necessary).

(ii) For each e in E n GK ,

η n Ge = η n Ge1 + η n Ge2 , with η n Ge1 = q n rh + ∇u n h + I n h k • b(u n-1 h )e z ) • n + L 2 (e) , (5.4
)

η n Ge2 = e q n rh + ∇u n h + I n h k • b(u n-1 h )e z (τ ) • n - u n h (τ ) dτ.
(5.5)

Despite their apparent complexity, these indicators only involve polynomials of low degree, so they are easy to compute. Moreover, the indicators η n Ge are not zero only for the e contained in Γ G , hence are not so many, and are the analogues of those proposed in [START_REF] Hild | Nicaise -A posteriori error estimations of residual type for Signorini's problem[END_REF] in a simpler case. With each element K in T n h , we associate the real number s K , which is positive in dimension d = 2 if K intersects Γ G and equal to zero otherwise. We are in a position to prove the next result.

max 0≤n≤N u n -u n h 2 L 2 (Ω) + 1 2 N m=1 τ m u m -u m h 2 H 1 (Ω) ≤ c N m=1 τ m K∈T m h h -2s K K ( η m K ) 2 + e∈E m GK η m Ge (5.6) + N m=1 (ε (I)m ) 2 + ε (D) h , (5.7)
where the quantity ε

(D) h is given by ε (D) h = n m=1 e∈E m KF τ m h 1-2s K e f m F -f m F h 2 L 2 (e) 1 2 + n m=1 e∈E m GK τ m h -2 e q m r -q m rh 2 L 2 (e) 1 2
, and the ε (I)m are defined in (5.1).

Proof. By integrating by parts on each K and using Cauchy-Schwarz inequalities, we derive

R n h ≤ K∈T n h α u n h -I n h u n-1 h τ n + b(u n h ) -I n h b(u n-1 h ) τ n -∇ • ∇u n h + I n h k • b(u n-1 h ) L 2 (K) u n -v h L 2 (K) + e∈E n K [ ∇u n h + I n h k • b(u n-1 h ) • n] e L 2 (e) u n -v h L 2 (e) + e∈E n F K ∇u n h + I n h k • b(u n-1 h )e z • n + f n L 2 (e) u n -v h L 2 (e) - e∈E n GK e q n r + ∇u n h + I n h k • b(u n-1 h e z ) (τ ) • n(τ )(u n -u n h )(τ ) dτ.
We treat separately the first three lines of R n h that we denote by R n 1 and the fourth line that we denote by R n 2 . 1) We insert the approximation of the datum f n , use the trace inequality (obtained by schwitching to the reference element and using [8, Thm 1.6.6])

u n -v h L 2 (e) ≤ c h 1 2 K |u n -v h | H 1 (K) + h -1 2 K u n -v h L 2 (K) ,
and also triangle and Young's inequalities. We also choose v h equal to u n h + w h , where w h is the function defined in Lemma A.1 of the Appendix with u raplaced by u n -u n h . All this gives

u n -v h L 2 (K) ≤ c h 1-s K K u n -u n h H 1 (K) , u n -v h H 1 (K) ≤ c h -s K K u n -u n h H 1 (K) .
Thus, we derive, with obvious definition of ε

(D)n h , R n 1 ≤ c K∈T n h (h -2s K K (η n K ) 2 + (ε (D)n h ) 2 ) 1 2 u n -u n h H 1 (Ω) .
2) To evaluate R n 2 , we observe that, since both u n and u n h are nonpositive on Γ G ,

R n 2 ≤ - e∈E n GK e q n r + ∇u n h + I n h k • b(u n-1 h e z ) (τ ) • n(τ )] + u n (τ ) dτ + e∈E n GK e q n r + ∇u n h + I n h k • b(u n-1 h ) (τ ) • n(τ )] -u n h (τ ) dτ,
whence, for a constant c that there also only depends on the data,

R n 2 ≤ c η n Ge1 + η n Ge2 .
In a final step, we replace n by m, multiply this equation by τ m and sum up on m. When combining the result with Lemma 5.1 and taking h 0 smaller than 1, we derive

α 2 u n -u n h 2 L 2 (Ω) + 1 2 n m=1 τ m u m -u m h 2 H 1 (Ω) ≤ c n m=1 τ m u m -u m h 2 L 2 (Ω) + R n ,
where the quantity R n is given by

R n = n m=1 (ε (I)m ) 2 + τ m 2 K∈T m h h -2s K K (η m K ) 2 + (ε (D)m h ) 2 ) + e∈E m GK η m Ge .
Thus, applying the discrete Gronwall's lemma (see [10, Chap. V, Lemma 2.4]) and noting that n m=1 τ m is smaller than T , we obtain the desired result.

Efficiency of the indicators.

We successively prove an upper bound first for the indicators η n K , next for the indicators η n Ge . For the set V 0 defined in (4.2), we compute the residual term R n h (w) for any

w in V 0 R n h (w) =α Ω (u n -u n h ) -(u n-1 -I n h u n-1 h ) τ n (x)w(x) dx + Ω (b(u n ) -b(u n h )) -(b(u n-1 ) -I n h b(u n-1 h )) τ n (x)w(x) dx + Ω ∇(u n -u n h ) + (k • b(u n-1 ) -I n h k • b(u n-1 h )e z (x) • ∇w(x) dx.
Indeed, the following result follows from Cauchy-Schwarz inequalities and the Lipschitz-continuity of b and k • b.

For brevity, we use the notation δu n for the quantity u n -u n-1 τn and its analogues δu n , δ(u -u h ) n .

Lemma 5.2. For all n, 1 ≤ n ≤ N , and for all w in W, the quantity R n h (w)) satisfies

|R n h (w)| ≤ |( δ(u -u h ) n L 2 (ωw) + ε (I)n 1 w L 2 (ωw) + |u n -u n h | H 1 (ωw) + c u n-1 -u n-1 h L 2 (ωw) + ε (I)n 2 |w| H 1 (ωw) ,
where ω w stands for the support of w and the quantities ε

(I)n 1
and ε

(I)n 2
are defined by

ε (I)n 1 = τ -1 n (Id -I n h )(u n h -u n-1 h ) L 2 (ωw) + τ -1 n (Id -I n h )b(u n-1 h ) L 2 (ωw) , ε (I)n 2 = (Id -I n h )k • b(u n-1 h ) L 2 (ωw) . (5.8)
To bound the indicators η n K , we take w in the space W defined in (4.6). The residual equation can easily be derived from (4.7):

R n h (w) = -α Ω u n h -I n h u n-1 h τ n (x)w(x) dx - Ω b(u n h ) -I n h b(u n-1 h ) τ n (x)w(x) dx + Ω ∇u n h + I n h k • b(u n-1 h )e z (x) • ∇w(x) dx + Γ F f F (τ, t n )w(τ ) dτ, ,1
or equivalently, by integrating by parts on each K,

R n h (w) = K∈T n h - K α u n h -I n h u n-1 h τ n + b(u n h ) -I n h b(u n-1 h ) τ n -∇ • ∇u n h + I n h k • b(u n-1 h )e z (x)w(x) dx - 1 2 e∈E n K e [∇u n h + I n h k • b(u n-1 h ) e z • n] e (τ )w(τ ) dτ - e∈E n F K e ∇u n h + I n h k • b(u n-1 h )e z (τ ) • n(τ ) + f n F (τ ) w(τ ) dτ .
The bound for η n K can easily be derived from this equation thank to standard choices of the function w.

Proposition 5.1. For 1 ≤ n ≤ N and for all K in T n h , the following bound holds for the indicators η n K defined in (5.2)

(??)η n K ≤ c h K δ(u -u h ) n L 2 (K) + ε (I)n 1 + |u n -u n h | H 1 (K) + u n-1 -u n-1 h L 2 (K) + ε (I)n 2 + ε (D)n , (5.9)
where ω K is the union of elements of T n h that share at least an edge (d = 2) or a face (d = 3) with K, the quantities ε (I)n 1 and ε (I)n 2 are defined in (5.8), and the quantity ε (D)n is given by

ε (D)n = e∈E KF h 1 2 e f n F -f n F h L 2 (e) .
Proof. We bound successively the three terms in η n K . 1) The idea consists in taking in (??) w equal to the function w K defined by

w K = -α u n h -I n h u n-1 h τ n + b(u n h ) -I n h b(u n-1 h ) τ n -∇ • ∇u n h + I n h k • b(u n-1 h )e z ψ K on K,
where ψ K is the bubble function on K (equal to the product of the barycentric coordinates associated with the vertices of K). Indeed, this yields

R n h (w K ) = α u n h -I n h u n-1 h τ n + b(u n h ) -I n h b(u n-1 h ) τ n -∇ • ∇u n h + I n h k • b(u n-1 h )e z ψ 1 2 K 2 L 2 (K) .
By combining this with the inverse inequalities, valid for any polynomial ϕ of fixed degree [START_REF]Verfürth -A Posteriori Error Estimation Techniques for Finite Element Methods[END_REF]Prop. 3.37]

ϕ L 2 (K) ≤ c ϕ ψ 1 2 K L 2 (K) , ϕ ψ K L 2 (K) ≤ ϕ L 2 (K) |ϕ ψ K | H 1 (K) ≤ c h -1 K ϕ L 2 (K)
, and using Lemma 5.3, we derive

h K α u n h -I n h u n-1 h τ n + b(u n h ) -I n h b(u n-1 h ) τ n -∇ • ∇u n h + I n h k • b(u n-1 h )e z L 2 (K) (5.10) ≤ c h K δ(u -u h ) n L 2 (K) + ε (I)n 1 + |u n -u n h | H 1 (K) + u n-1 -u n-1 h L 2 (K) + ε (I)n 2
.

(5.11)

2) For each edge or face e of an element K, we introduce a lifting operator L e,K which • maps polynomials on e vanishing on ∂e into polynomials on K vanishing on ∂K \e, • and is constructed from a fixed operator on the reference element. Next, for each edge or face e shared by two elements K and K , we set

w e = L e,κ -[ ∇u n h + I n h k • b(u n-1 h )e z • n] e ψ e onκ ∈ {K, K }, = 0 on Ω \ (K ∪ K ),
where ψ e is now the buble function on e. By using the property of the operator L e,κ which is easily derived by switching to the reference element (recall that ϕ is any polynomial of fixed degree)

L e,κ ϕ L 2 (κ) + h κ |L e,κ ϕ| H 1 (κ) ≤ h 1 2
e ϕ L 2 (e) , the new inverse inequalities [START_REF]Verfürth -A Posteriori Error Estimation Techniques for Finite Element Methods[END_REF]Prop. 3.37]

ϕ L 2 (e) ≤ c ϕ ψ 1 2 e L 2 (e) ,
ϕ ψ e L 2 (e) ≤ ϕ L 2 (K) , and (5.10), we derive

h 1 2 e [ ∇u n h + I n h k • b(u n-1 h ) • n] e L 2 (e) ≤c κ∈{K,K } h κ δ(u -u h ) n L 2 (κ) + ε (I)n 1 + c (|u n -u n h | H 1 (K∪K ) + u n-1 -u n-1 h L 2 (K∪K ) + ε (I)n 2
.

3) To bound the terms h

1 2 e ∇u n h + I n h k • b(u n-1 h )e z • n + f n F h L 2 (e)
, if e is contained in K, we take w equal to we , with

we = L e,κ -∇u n h + I n h k • b(u n-1 h )e z ) • n + f n F h ψ e on K, 0 on Ω \ K.
(5.12)

We use exactly the same arguments as in the previous step and obtain the same estimate only with the further term h e) in the right-hand side. All these bounds give estimate (??).

1 2 e f n F -f n F h L 2 (

Finally, to bound the indicators η n

Ge , we use the quantity R n h (w) for w in V 0 . For simplicity, we define the set Ω n Gh as the domain made of Γ G and the interior of all triangles K in T n h having at least an edge (d = 2) or a face (d = 3) lying on Γ G . We denote T n Gh the set of triangles contained in Ω n Gh and we assume (this is not restrictive) that two elements in T n Gh share at most a vertex (d = 2) or an edge (d = 3). Thus, using problem (3.2) and integrating by parts on each K yields for any w having its support in Ω

n Gh R n h (w) ≥ K∈T n Gh - K α u n h -I n h u n-1 h τ n + b(u n h ) -I n h b(u n-1 h ) τ n -∇ • ∇u n h + I n h k • b(u n-1 h e z ) (x)w(x) dx - e∈E n GK e ∇u n h + I n h k • b(u n-1 h e z + q n r ) (τ ) • n(τ )w(τ ) dτ .
Proposition 5.2. For 1 ≤ n ≤ N and for all K in T n h and each e in E n GK , the following bound holds for the indicators η n Ge1 defined in (5.4)

η n Ge1 ≤c h K δ(u -u h ) n L 2 (K) + ε (I)n 1 + |u n -u n h | H 1 (K) + u n-1 -u n-1 h L 2 (K) + ε (I)n 2 + ε (D)n , (5.13)
where the quantities ε

(I)n 1 and ε (I)n 2
are defined in (5.1), and the quantity ε (D)n is given by

ε (D)n = e∈E GK h 1 2
e (q n rq n rh ) • n L 2 (e) .

(5.14)

Proof. With the same notation as in the previous proof, we take w equal to we , with we = L e,K -

q n rh + ∇u n h + I n h k • b(u n-1 h )e z • n + ψ e on K, 0 on Ω \ K. (5.15)
It can be noted that w is nonpositive on e, so that it belongs to V 0 . By combining the same inverse inequalities as previously together with a Cauchy-Schwarz inequality and Proposition 5.4, we obtain the estimate for η n Ge1 .

A slightly different argument is needed to bound η n Ge2 . Proposition 5.3. For 1 ≤ n ≤ N and for all K in T n h and each e in E n GK , the following bound holds for the indicators η n Ge2 defined in (5.4) is well-defined).

η n Ge2 ≤c h K δ(u -u h ) n L 2 (K) + ε (I)n 1 + c u n -u n h H 1 (K) + u n-1 -u n-1 h H 1 (K) + ∂ n (u n -u n h ) H -1 2 (e)
To conclude, we deduce from the formula a = [a] + + [a] -that η n Ge2 ≤ E + c η n Ge1 , and we use estimate (5.13).

Due to the presence of s K , estimate (5.6) is not fully optimal, but the lack of optimality seems negligible (s K is equal to zero in dimension d = 3 and also in nearly all the triangles in dimension d = 2 and is any positive number elsewhere). Estimate (5.16) is also not optimal because of the term involving the normal derivative. However, estimates (??) (5.13) and (5.16) are local (i.e. only involve the error in a neighborhood of K), so that they are an efficient tool for mesh adaptivity.

Numerical results

In this section we consider the well-known test case of Abdull and Gilham ([]). The experiment consists in a Plexiglas trapezoidal sandbox filled with medium-fine sand. The domain is represented at figure 6 where the dimension are indicated. Note that the surface slope is uniform equal to 14 o . The initial water table is assigned at the height of 0.78 m.

The volumetric water content given by the Van Genuchten formula The conductivity is given by

θ(ψ) = θ s (1 + (ε|h|) n ) -m if h < 0, θ s if h ≥ 0,
K(h) = Ks(1-(ε|h|) n-1 )(1+(ε|h|) n ) -m ) 2 (1+(ε|h|) n ) m/2 if h < 0, K s if h ≥ 0,
with K s = 5.10 -4 cm s -1 . A constant rainfall rate of 0.1 K s is applied on the whole surface domain, and no flow boundary is enforced on the bottom and on the side walls. We denote by p the percentage of the underground porous media which is unsaturated and we set P = {25%, 15%, 10%, 5%, 1%}.

The time-space adaptive algorithm we propose is the following.

Set T 0 h , τ 0 , u 0 h , n = 1. While t ≤ T : i) Compute u n h , η n τ , η n , η n K , η n Ge , η n Ge2 on T n h .

ii) if η n τ ≥ K∈T n We present below the refined meshes obtained thanks to our time-mesh adaptive algorithm.

It also follows from a Hölder's inequality that, with 1 q + 1 q = 1,

|α c | ≤ c h -d Kc h d q
Kc u L q (Kc) = c h -d q Kc u L q (Kc) . Using the imbedding of H 1 (K c ) into L q (K c ) for all q < ∞ in dimension d = 2 and for q = 6 in dimension 3 gives the desired result.

Lemma 5 . 1 .

 51 Assume that the function b is of class C 2 with bounded derivatives. The quantities R n h satisfy for all n, 1 ≤ n ≤ N , n m=1

Theorem 5 . 1 .

 51 The following a posteriori estimate holds between the solution (u n ) 0≤n≤N of problem (3.1)-(3.2) and the solution (u n h ) 0≤n≤N of problem (3.4)-(3.5)

1 hH 1 ( 2 +

 112 (5.1), and the quantity ε (D)n is defined in(4.11).Proof. We first bound the quantityE = e q n rh + ∇u n h + I n h k • b(u n-1 h )e z (τ ) • n(τ ) u n h (τ ) dτ.We deduce from (4.11) (with v = 2u n )) that E ≤ e q n rh + ∇u n h + I n h k • b(u n-1 h )e z (τ ) • n(τ )u n h (τ ) dτe q n r + ∇u n + k • b(u n-1 )e z (τ ) • n(τ )u n (τ ) dτ ≤ e q n rh + ∇u n h + I n h k • b(u n-1 h )e z (τ ) • n(τ )(u n h -u n )(τ ) dτ + e q n rhq rn + ∇(u n h -u n ) -I n h k • b(u n-1 h + k • b(u n-1 )e z (τ ) • n(τ )u n (τ ) dτCauchy-Schwarz inequalities combined with the Lipschitz-continuity of k • b yieldE ≤ c u n -u n h H 1 (K) + u n-1 -u n-K) + ∂ n (u n-1 -u n-1 h ) ε (D)n ,(it can be noted from problem (3.2) that ∆u n belongs to L 2 (Ω), so that the quantity ∂ n (u n-1u

Figure 1 .

 1 Figure 1. The computational domainθ s = 0.55, n = 1.9, m = 1 -1/n, ε = 0.036cm -1 .The conductivity is given by

  + η n Ge2 , then τ n+1 = 0.5τ n else τ n+1 = τ n . iii) if p ∈ P, T n+1 h is derived from T nh by 2 iterations of the standard space adaptive algorithm based on an equirepartion of the error indicators. Else T n+1 h = T n h . iv) t := t + τ n , n := n + 1.

  and equal to the solutions u n and u n h of problems (3.1)-(3.2) and (3.4)-(3.5), respectively, in t n , 0 ≤ n ≤ N . Theorem 3.2. For any data u B , f F , q r and u 0 satisfying (2.6) and (3.8), there exists a subsequence of the (u hτ ) τ,h which converges to the solution u of problem (2.1) weakly in X.

,1 Figure 2. Mesh at time t = 0, t = 200 and t = 500 ,1

Appendix

We prove here a technical lemma. Recall that the quantity s K is positive in dimension d = 2 if K intersects Γ G and equal to zero otherwise. Lemma A.1. For any function u in H 1 (Ω) vanishing on Γ B , there exists a function w h in W h such that, for any

where ∆ K is the union of all elements of T n h that intersect K. Proof. We choose w n equal to w 1 h + w 2 h , where w 1 h is the image of u by any type of Clément regularization operator preserving the nullity on Γ B and w 2 h has its support contained in the set Ω n Gh (made of Γ G and all triangles K in T n h having at least a corner in Γ G ). Thus estimate (A.1) holds (with s K = 0) on all triangles K that does not intersect Γ G and also on the other triangles with w h replaced by w 1 h . On the other hand, let C n be the set of the corners of elements of T n h which belong to Γ G . For each c in C n , we denote by K c an element of T n h containing c and by ϕ c the Lagrange function associated with c (with values in V n h ). We assume without restriction that

and we take

Thus, w h belongs to W h . To evaluate the norms of w 2 h we first observe by switching to the reference element that ϕ c L 2 (Ω) ≤ c h


Laura Ciobanu

Alexander Kolpakov

Geodesic growth of right-angled Coxeter groups based on trees

Keywords: 2000 Mathematics Subject Classification: 20E08, 20F65 geodesic growth, Coxeter group, RACG, regular languages. "And

Introduction

The geodesic growth function of a group G with respect to a finite generating set S counts, for each positive integer n, the number of geodesics of length n starting at the identity 1 G in the Cayley graph of G with respect to S. The geodesic growth series is the formal power series that takes the values of the geodesic growth function as its coefficients (see Definition 2.4).

The groups that we consider in this paper, right-angled Coxeter groups, or RACGs, are known to have a regular language of geodesics with respect to the standard generating sets, and therefore rational geodesic growth series (see [START_REF] Loeffler | Graph products and Cannon pairs[END_REF] or [START_REF] Bjorner | Combinatorics of Coxeter groups[END_REF]Theorem 4.8.3] for proofs of these facts). Our goal here is to obtain more specific data concerning the geodesic growth of RACGs. Namely, we are interested in extracting information about the geodesic growth series from the defining graph of the group. It is known that non-isomorphic graphs define non-isomorphic RACGs [START_REF] Radcliffe | Rigidity of graph products of groups[END_REF], and that non-isomorphic RACGs can have equal geodesic growth [START_REF] Antolín | Geodesic growth in right-angled and even Coxeter groups[END_REF]. However, we are interested in knowing how much the similarities or differences between two defining graphs influence the geodesic growth of the corresponding RACGs.

Two RACGs G i = G(Γ i ), i = 1, 2, with non-isomorphic defining graphs Γ i , may have equal standard growth (see Definition 2.4 (2)); this can be determined by computing the fpolynomials of the graphs Γ i [START_REF] Davis | The geometry and topology of Coxeter groups[END_REF]Proposition 17.4.2]. However, their geodesic growth exhibits more subtle properties [START_REF] Antolín | Geodesic growth in right-angled and even Coxeter groups[END_REF], and in general it is not known which graph theoretic conditions completely determine the geodesic growth of a RACG. The examples of non-isomorphic graphs Γ i defining RACGs with equal geodesic growth in [START_REF] Antolín | Geodesic growth in right-angled and even Coxeter groups[END_REF] are degree-regular and have cycles. In this paper we consider the case when the Γ i 's are trees. Although these are some of the simplest classes of graphs, we already encounter a phenomenon that shows a great difference between the behaviour of the standard growth and that of the geodesic growth. The standard growth of each RACG G i is determined solely by the number of vertices (or edges) in the respective tree Γ i , while the geodesic growth can be distinct even for two co-spectral trees (which might have co-spectral complements as well). Recall that two graphs are co-spectral if the characteristic polynomials of their adjacency matrices are the same.

In this paper we count and compute with the help of an automaton generating the geodesic language in a RACG based on a tree T . This automaton reflects some of the path information from T . In general, a good deal of combinatorial information about T can be extracted from its spectrum. However, we encounter two rather different behaviours: on one hand we produce (infinitely many) pairs of trees T 1 and T 2 which are non-isomorphic, co-spectral, with co-spectral complements, and the RACGs G(T i ), i = 1, 2, based on them have the same geodesic growth; on the other hand we obtain (infinitely many) pairs of trees S 1 and S 2 which are non-isomorphic and co-spectral, with co-spectral complements, whose respective RACGs G(S i ) have distinct geodesic growth. This shows that the spectrum of a tree alone does not determine the geodesic growth of the RACG based on that particular tree.

The following two theorems, proved in Sections 3 and 5, respectively, are the main results of the paper: Theorem 1.1. There exist two families of trees We note that McKay showed in [START_REF] Mckay | On the spectral characterisation of trees[END_REF] that the trees T 1 i and T 2 i in Theorem 1.1 are simultaneously co-spectral, have co-spectral complements and co-spectral line graphs. On the contrary, the trees S 1 i and S 2 i from Theorem 1.2 might have line graphs with distinct spectra. Moreover, by a result of McKay [START_REF] Mckay | On the spectral characterisation of trees[END_REF], in both Theorem 1.1 and Theorem 1.2 the cardinality of each family T i n and S i n tends asymptotically to the cardinality of the set Υ n of all trees on n vertices:

T 1 n = {T 1 1 , T 1 2 , . . . } and T 2 n = {T 2 1 , T
(1) card T i n card Υ n → 1 and card S i n card Υ n → 1, as n → ∞, i = 1, 2.
The two theorems above, together with McKay's results, lead us to the following conjecture:

Conjecture 1.3. If two trees are simultaneously co-spectral, have co-spectral complements and co-spectral line graphs, then they have the same geodesic growth.

Section 4 contains details about the computation of several kinds of geodesics, and their numbers, for the trees introduced in Section 3.

The authors have created a Python code Monty 1 [START_REF] Kolpakov | SAGE" worksheet, available on-line at[END_REF], which performs the computations needed for the results in this paper by using either SAGE standard routines for symbolic computation of determinants and rational expressions of growth series, or the Berlekamp-Massey algorithm for restoring the rational expression for a growth series from a sufficient number of its coefficients, as a more efficient approach. Our Python code constructs finitestate automata that accept the geodesic languages in RACGs based on triangle-free graphs (in this case, trees), and then proceeds to determine the respective growth series.

Definitions and notation

Let S be a finite set and S * the free monoid on S. We identify S * with the set of words over S, that is, finite sequences of elements of S. We use |.| to denote word length.

Let G = S be a group generated by S. For an element g of G, denote by |g| S the word length of g with respect to S. Given w ∈ S * , we denote by w the image of w in G under the natural projection S * → G. Definition 2.1. A word w over an alphabet S is geodesic in G = S if |w| = |w| S . The set of geodesics in G with respect to S will be denoted by Geo(G) or Geo(S).

Let Γ = Γ(G, S) be a simple (no loops, no multiple edges) graph with vertex set V (Γ) = S (or simply V ) and edge set E(Γ) (or simply E), where E ⊆ V × V . The RACG based on Γ is given by the presentation s ∈ S | s 2 = 1 ∀s ∈ S, and (ss

) 2 = 1, ∀{s, s } ∈ E .
It is easy to see that for any two involutions s and s the relation (ss ) 2 = 1 implies ss = s s. This leads to another possible presentation for RACGs: s ∈ S | s 2 = 1, ss = s s, ∀{s, s } ∈ E .

In the present paper we use the same letters for the vertices of Γ and the corresponding generators of the group G(Γ).

The star of a vertex v ∈ V in Γ, denoted by St Γ (v), or St(v) if the ambient graph is clear in the given context, is the set of vertices in Γ that are adjacent to v. That is,

St Γ (v) = {w ∈ V | {v, w} ∈ E}.
We now describe a finite deterministic automaton that recognises geodesics in RACGs (see [START_REF] Hopcroft | Introduction to automata theory, languages, and computation[END_REF] for definitions of languages and automata). We define such automata in the standard way, as quintuples (Q, Σ, δ, q 0 , F ), where Q is the finite set of states, Σ the input alphabet, δ the transition function, q 0 the initial or start state, and F the set of final or accepting states. The following definition is a simplified version of Proposition 4.1 in [START_REF] Antolín | Geodesic growth in right-angled and even Coxeter groups[END_REF]. Definition 2.2. Let Γ = (V, E) be a tree. The deterministic finite state automaton recognising the geodesics in

G(Γ) is A = (Q, S, δ, {∅}, F ), where the set of states is Q = {∅} ∪ V ∪ E ∪ {ρ},
ρ is the unique "fail" state, and {∅} is the start state. The input alphabet is S = V , the set of accept states F is ∅ ∪ V ∪ E, and the transition function

δ : Q × S → Q is given by (1) δ(σ, s) = (St(s) ∩ σ) ∪ {s}, for s / ∈ σ; (2) δ(σ, s) = ρ, otherwise.
Definitions 2.3. Any set L of words over an alphabet Σ gives rise to a strict growth function

f L : N → N, defined by f L (n) := |{W ∈ L | |W | = n}|.
1 this is not the given name of the code, which would be obviously too posh for such a petty thing, but a reference name, which is seemingly good for any Python code.

Definitions 2.4. Let G be a group generated by S. Then we define the following:

(1) The geodesic growth function f Geo(G) : N → N is given by

f Geo(G) (r) := f Geo(G,S) (r) = |{w ∈ S * | |w| = |w| S = r}|,
and the geodesic growth series of G equals

G (G,S) (t) = ∞ r=0 f Geo(G) (r) t r .
(2) The spherical (standard) growth function σ (G,S) : N → N is given by

σ (G,S) (r) = |{g ∈ G | |g| S = r}|,
and the spherical (standard) growth series of G equals

Σ G (t) := Σ (G,S) (t) = ∞ r=0 σ (G,S) (r) t r .
We recall that the f -polynomial of a graph Γ is the generating function for the number of cliques (that is, complete subgraphs) of size i in Γ:

f (t) = f 0 + f 1 t + f 2 t 2 + . . .
, where f i is the number of i-cliques in Γ. We consider the empty set to be a clique on zero vertices, and therefore f 0 = 1. We remark that the spherical (standard) growth function of a RACG is determined by the f -polynomial of its defining graph [START_REF] Davis | The geometry and topology of Coxeter groups[END_REF]Proposition 17.4.2]. We want to contrast this to the fact that the geodesic growth function is not uniquely determined by the f -polynomial, or even by the spectrum of the defining graph, as the following sections show.

RACGs with equal geodesic growth

In this section we prove Theorem 1.1 by giving an explicit construction of the families of trees T i n , i = 1, 2. Definition 3.1. We define the coalescence τ • σ of two rooted trees τ and σ to be the tree which results from merging τ and σ at their roots. The tree τ • σ has as vertex set the union of the vertex sets of τ and σ, and as root the identification of the roots of τ and σ, as shown in Fig. 1.

Consider the trees T 1 and T 2 , both rooted at 0, first described by McKay in [START_REF] Mckay | On the spectral characterisation of trees[END_REF], as shown in Fig. 2. In [START_REF] Mckay | On the spectral characterisation of trees[END_REF] the following fact is proved: Theorem 3.2. Let T be a rooted tree with at least two vertices and with root labelled 0. Then the trees Γ i = T • T i , i = 1, 2, are not isomorphic, but are co-spectral. Also, their complements Γ i and line graphs

L(Γ i ), L(Γ i ), L(Γ i ), i = 1, 2, are respectively co-spectral. Now let G 1 := G(T 1
) and G 2 := G(T 2 ) be the RACGs associated to the trees T 1 and T 2 defined above. First note that since the trees T 1 and T 2 are isomorphic as graphs, the groups G 1 and G 2 are isomorphic. However, T 1 and T 2 are not isomorphic as rooted trees.

Let τ be a tree with n vertices. Fix a labelling {0, . . . , n-1} of the vertices of τ , and suppose that 0 represents τ 's root. Define Γ 1 = τ • T 1 and Γ 2 = τ • T 2 to be the trees obtained as the coalescence of τ with T i at vertex 0, and let G(Γ 1 ) and G(Γ 2 ) be the RACGs based on Γ 1 and Γ 2 , respectively. Since Γ 1 and Γ 2 are non-isomorphic, G(Γ 1 ) and G(Γ 2 ) are non-isomorphic, as well.

The following lemma is the key ingredient of Theorem 1.1. In order to simplify the exposition in the proof of Lemma 3.3 we use the notation below.

Notation 3.4. Let G be a group with generating set T containing the letter 0. Denote the set of words in Geo(T ) starting with 0 by Geo 0 (T ), the set of words in Geo(T ) ending in 0 by Geo 0 (T ), and the set of words in Geo(T ) ending and starting with 0 by Geo 0 0 (T ).

Proof. (of Lemma 3.3) Notice that f Geo(Γ 1 ) (r) (respectively f Geo(Γ 2 ) (r)) is equal to the number of all words of length r in Γ * 1 (respectively, Γ * 2 ) minus the number of those words of length r in Γ * 1 that are not geodesics. We denote the number of non-geodesics by f Geo(Γ 1 ) (r) (and

f Geo(Γ 2 ) (r), respectively). Since |Γ 1 | = |Γ 2 |, clearly f Geo(Γ 1 ) (r) = f Geo(Γ 2 ) (r) if and only if f Geo(Γ 1 ) (r) = f Geo(Γ 2 ) (r).
We now show that f Geo(Γ 1 ) (r) = f Geo(Γ 2 ) (r) for all r ≥ 1. Any word in Γ * 1 can be written as w 1 u 1 w 2 . . . u n , where u i ∈ (τ \ {0}) * and w i ∈ T * 1 . A non-geodesic word w in Γ * 1 belongs to one of the following sets (or is of the type), depending on its form:

A : w contains non-geodesics u i ∈ (τ \ {0}) * or w j ∈ T * 1 , or both, for some 1 ≤ i, j ≤ n, or B : all w i and u i are geodesic on their respective alphabets, and there exists 1 ≤ j < n such that 0u j 0 is a subword of w with u j ∈ (St Γ 1 (0) ∩ τ ) * , or C : all w i and u i are geodesic on their respective alphabets, and w contains a subword of the form s0s, where s

∈ (St Γ 1 (0) ∩ τ ) * .
Notice that the set of non-geodesics is then

A ∪ B ∪ C, where A ∩ (B ∪ C) = ∅ and B ∩ C = ∅.
The number of words in A depends only on the geodesic growth series of τ and T 1 , and thus it will be equal to the number of words of type A in Γ * 2 . The computations in Section 4 show that the following identities hold:

f Geo(T 1 ) (r) = f Geo(T 2 ) (r), f Geo 0 (T 1 ) (r) = f Geo 0 (T 1 ) (r), f Geo 0 (T 1 ) (r) = f Geo 0 (T 2 ) (r) and f Geo 0 0 (T 1 ) (r) = f Geo 0 0 (T 2 ) (r
) for all r ≥ 1, as a result of (2). This means that there is a length-presenting bijection φ 0 between Geo 0 (T 1 ) and Geo 0 (T 2 ), i.e. for each geodesic w = 0v ∈ Geo 0 (T 1 ) there is a geodesic w = 0v = φ 0 (w) ∈ Geo 0 (T 2 ), and |w| = |w |. Analogously, there is a length-preserving bijection φ 0 : Geo 0 (T 1 ) → Geo 0 (T 2 ), and a bijection φ 0 0 : Geo 0 0 (T 1 ) → Geo 0 0 (T 2 ). This means there is a length-preserving bijection φ : Geo 0 (T 1 ) ∪ Geo 0 (T 1 ) → Geo 0 (T 2 ) ∪ Geo 0 (T 2 ) between those geodesics starting or ending with 0, for which we have that φ| Geo 0 = φ 0 , φ| Geo 0 = φ 0 and φ| Geo 0 0 = φ 0 0 . By the computations in (2) we have that f (Geo 0 ∪Geo 0 )(T 1 ) (r) = f (Geo 0 ∪Geo 0 )(T 2 ) (r) by the standard formula for the cardinality of the union of two sets, and since f Geo(T 1 ) (r) = f Geo(T 2 ) (r), we also have that f Geo\(Geo 0 ∪Geo 0 )(T 1 ) (r) = f Geo\(Geo 0 ∪Geo 0 )(T 2 ) (r). Thus, there is a length-preserving bijection ψ between Geo \ (Geo 0 ∪ Geo 0 )(T 1 ) and Geo \ (Geo 0 ∪ Geo 0 )(T 2 ).

The bijection φ can be extended to the set Geo(T 1 ) of all geodesics on T 1 by letting φ(w) = ψ(w) for all w ∈ Geo \ (Geo 0 ∪ Geo 0 )(T 1 ), and then furthermore extended to Geo(T 1 ) ∪ Geo(τ \ {0}) by letting φ(w) = w for all w ∈ Geo(τ \ {0}).

Then φ provides a bijection between the words of type B in Γ * 1 and the words of type B in Γ * 2 . To see this, associate to each w 1 u 1 w 2 . . . u n the word φ(w 1 )φ(u 1 ) . . . φ(w n )φ(u n ) = φ(w 1 )u 1 . . . φ(w n )u n . Then w i ends in 0, w i+1 starts with 0, and u i ∈ St τ (0) * if and only if φ(w i ) ends in 0, φ(w i+1 ) starts with 0, and φ(u i ) = u i ∈ St τ (0) * , by definition.

It is immediate to see that φ also provides a bijection between the words of type C, and between the words of type B ∩ C in Γ 1 and Γ 2 , respectively. Thus, there is a length-preserving bijection between the words of type A ∪ B ∪ C (i.e. all non-geodesic words) in Γ 1 and Γ 2 . This concludes the proof of the lemma.

Proof. (of Theorem 1.1) Let Υ k = {τ 1 , τ 2 , . . . } be the set of non-isomorphic trees on k ≥ 2 vertices. For i = 1, 2 the two families of trees

T i n = {τ • T i |τ ∈ Υ k }, k = n -15, satisfy Theorem 1.1.
From Theorem 3.2 we already know that Γ 1 = τ • T 1 and Γ 2 = τ • T 2 are co-spectral, for any τ ∈ Υ k . By Lemma 3.3 the groups G(Γ 1 ) and G(Γ 2 ) have the same geodesic growth series. Theorem 3.5 (Lemma 4.3 in [START_REF] Mckay | On the spectral characterisation of trees[END_REF]). Let p i (n), i = 1, 2, be the proportion of trees on n vertices that have T i as a limb. Then p 1 (n) = p 2 (n) for all n and lim n→∞ p i (n) = 1. From Theorem 3.5, we obtain that card T i n card Υn → 1, as n → ∞.

Computing the numbers of special geodesics

In this section we prove, by concrete computations, the equalities between the numbers of special geodesics required in the proof of Lemma 3.3.

Let T 1 and T 2 be as in Figure 2, and recall Notation 3.4. For each group G i = G(T i ), i = 1, 2, we construct a finite automaton A i accepting the geodesic language of G i , as described in Definition 2.2, and with the help of A i we compute the growth series

γ G 1 (t) = r≥0 f Geo(T 1 ) (r) t r , 0 γ G 1 (t) = r≥0 f Geo 0 (T 1 ) (r) t r , 0 γ G 1 (t) = r≥0 f Geo 0 (T 1 ) (r) t r , and 0 0 γ G 1 (t) = r≥0 f Geo 0 0 (T 1 ) (r) t r .
We also compute the analogous growth series γ G 2 (t), 0 γ G 2 (t), 0 γ G 2 (t) and 0 0 γ G 2 (t) for G 2 , where the series coefficients are given by the sequences Geo(T 2 )(r), Geo(T 2 ) 0 (r), Geo(T 2 ) 0 (r) and Geo(T 2 ) 0 0 (r), respectively. Our computations, which we elaborate upon below, show that

(2) γ G 1 (t) = γ G 2 (t), 0 γ G 1 (t) = 0 γ G 2 (t), 0 γ G 1 (t) = 0 γ G 2 (t) and 0 0 γ G 1 (t) = 0 0 γ G 2 (t).
These identities prove the equality of the corresponding numbers of geodesics.

Let A be a deterministic finite-state automaton with accepting states q i , i = 0, 1, . . . , N , where q 0 is the start state and the "fail" state is denoted by q. Let M = M (A) be the transition matrix of the automaton A. Computing the generating function γ A (t) of A is a standard technique, and the formula for γ A (t) is

(3) γ A (t) = e T M w det(I -tM ) ,
where I is the N × N identity matrix, M is the adjoint matrix of I -tM , e = (1, 0, . . . , 0) T and w = (1, 1, . . . , 1) T are two vectors in Z N . Now let A i be the deterministic finite-state automata accepting the language of geodesics in G i , i = 1, 2, and let M i = M (A i ) be the transition matrix of A i . Then the geodesic growth series γ G i (t) is a rational function (determined by the equality (3)), and the coefficients of its numerator and denominator can be easily computed, see [START_REF] Epstein | Growth functions and automatic groups[END_REF]. The Python code Monty [START_REF] Kolpakov | SAGE" worksheet, available on-line at[END_REF] may perform the above computation either by finding the symbolic determinant det(I -tM i ) (usually slow), or by applying the Berlekamp-Massey algorithm (a faster one). This Python code starts by creating the finite-state automaton A i , given a triangle-free graph T i (in this case, a tree), and then proceeds to determining γ G i (t).

Since G 1 and G 2 are isomorphic, the equality

γ G 1 (t) = γ G 2 (t) is immediate.
Here we provide an explicit formula for these identical growth series. The output of Monty for both trees T 1 and T 2 consists of two finite-state automata A i which are isomorphic, as one only renumbers the vertices of T 2 in order to obtain T 1 . Each A i has 32 states and 466 transition arrows, and the corresponding growth series are

γ G 1 (t) = γ G 2 (t) = (1 + t)(1 + 2t -2t 3 -4t 4 -t 5 )(1 + 5t + 10t 2 + 9t 3 -5t 4 -26t 5 -34t 6 -22t 7 -t 8 + 7t 9 + 4t 10 )(1 -8t -85t 2 -243t 3 -222t 4 + 332t 5 + 1194t 6 + 1349t 7 + 132t 8 -1510t 9 -2008t 10 -1088t 11 + 28t 12 + 359t 13 + 170t 14 + 15t 15 ) -1 .
Now we compute the functions 0 γ G i (t), i = 1, 2, which will turn out to be equal, as well. However, in this case the equality is not known to hold beforehand: even though the groups G 1 and G 2 are isomorphic, the image of vertex 0 in T 1 under the canonical isomorphism does not correspond to vertex 0 in T 2 . We shall explicitly compute the generating function 0 γ G i (t) for the number of geodesic words starting with 0 accepted by each A i , i = 1, 2, and then compare these functions. Suppose that the start state of each A i is q 0 = ∅. Let the corresponding transition function δ i be so that δ i (q 0 , 0) = q k i = {0}. If a word w labels a path from q k i to an accept state, then the word 0w is a geodesic word in G i starting with 0. Thus, we have to compute the generating function 0 α G i (t) for the number of words starting at q k i and ending at an accept state:

0 α G i (t) = e T M i w det(I -tM i )
,

where e = (0, . . . , 0, 1

k i , 0, . . . , 0) T and w = (1, 1, . . . , 1) T . Then we use the fact that 0 γ G i (t) = t • 0 α G i (t). By symmetry, we get 0 γ G i (t) = 0 γ G i (t), i = 1, 2.
By using Monty we obtain

0 γ G 1 (t) = 0 γ G 2 (t) = t(1 + t)(1 + 2t -2t 3 -4t 4 -t 5 )(1 + 4t + 4t 2 -3t 3 -9t 4 -5t 5 + 3t 6 + t 7 -3t 8 -3t 9 )(1 -8t -85t 2 -243t 3 -222t 4 + 332t 5 + 1194t 6 + 1349t 7 + 132t 8 -1510t 9 -2008t 10 -1088t 11 + 28t 12 + 359t 13 + 170t 14 + 15t 15 ) -1 .
Finally, it remains to compute the growth series 0 0 γ G i (t). There are two kinds of words forming disjoint subsets of the geodesic language of G i that we are interested in:

(I) the words w starting at the accept state q k i = δ(q 0 , 0) = {0} of A i and coming back to it: then the geodesic word 0w starts and ends with a "0" (since there are only 0-transitions leading to q k i = {0} and no 0-transition coming out of q k i ), by Definition 2.2 (i). (II) the words w starting at the state q k i = {0} and ending at a state q l i such that δ(q l i , 0) = q m i = q k i and q m i is an accept state: then the word 0w0 will be a geodesic word starting and ending with a "0". Let 0 0 α G i (t) be the generating function for the words of type (I), and 0 0 β G i (t) be that for the words of type (II). Then,

0 0 γ G i (t) = t • 0 0 α G i (t) + t 2 • 0 0 β G i (t). We have that 0 0 α G i (t) = e T M i w det(I -tM i )
,

with e = w = (0, . . . , 0, 1

k i , 0, . . . , 0) T .
Analogously,

0 0 β G i (t) = e T M i w det(I -tM i )
,

with e = (0, . . . , 0, 1 k i , 0, . . . , 0) T and w having a "1" at position l i for all states q l i as described above, and zeroes at all other places.

By using Monty we obtain

0 0 γ G 1 (t) = 0 0 γ G 2 (t) = t • (1 -5t -94t 2 -374t 3 -456t 4 + 955t 5 + 4275t 6 + 5652t 7 -1617t 8 -16773t 9 -24255t 10 -7337t 11 + 26583t 12 + 45100t 13 + 26181t 14 -12789t 15 -34553t 16 -24957t 17 -3147t 18 + 8130t 19 + 6288t 20 + 1398t 21 -458t 22 -284t 23 -24t 24 ) • (1 + 3t + 2t 2 -3t 3 -9t 4 -8t 5 + 4t 7 + 3t 8 -t 9 ) -1 (1 -8t -85t 2 -243t 3 -222t 4 + 332t 5 + 1194t 6 + 1349t 7 + 132t 8 -1510t 9 -2008t 10 -1088t 11 + 28t 12 + 359t 13 + 170t 14 + 15t 15 ) -1 .
Thus, the numbers of special geodesics in G 1 and G 2 coincide.

RACGs with different geodesic growth

In this section we give a proof of Theorem 1.2. Our construction will be analogous to that in Section 3, although we shall use different trees, S 1 and S 2 , in order to construct the families

S i n = {τ • S i |τ ∈ Υ k }, k = n -9.
Namely, we will use the rooted trees in Fig. 3. We set Γ 1 = τ • S 1 and Γ 2 = τ • S 2 , where τ is an arbitrary tree with n vertices labelled {0, . . . , n -1}. Suppose that 0 represents τ 's root. Then the following holds: Theorem 5.1 (page 27 of [START_REF] Godsil | Are Almost All Graphs Cospectral?[END_REF]). The trees Γ 1 and Γ 2 are not isomorphic, but are co-spectral and have co-spectral complements. Their line graphs L(Γ i ), L(Γ i ) and L(Γ i ) are not necessarily co-spectral.

Proof. The trees Γ 1 and Γ 2 are not isomorphic, since Γ 1 has the rooted tree σ depicted in Fig. 4 as a limb more times than Γ 2 does.

The characteristic polynomials of τ , S i and Γ i (i.e. the characteristic polynomials of the adjacency matrices of these graphs) are related by

φ Γ i (t) = φ τ (t) • φ S i \{0} (t) + φ τ \{0} (t) • φ S i (t) -t • φ τ \{0} (t) • φ S i \{0} (t),
according to [13, Lemma 2.2 (i)]. Given that the trees S 1 and S 2 are isomorphic as graphs (though not as rooted trees), we obtain that Γ 1 and Γ 2 share the same characteristic polynomial (i.e. are co-spectral).

The fact that Γ 1 and Γ 2 are co-spectral follows from the above and [13, Theorem 3.1 (i)]. Now, by letting τ be the line graph on two vertices we obtain that L(Γ 1 ), L(Γ 1 ) and L(Γ 1 ) have different characteristic polynomials from L(Γ 2 ), L(Γ 2 ) and L(Γ 2 ), respectively. Lemma 5.2. Let τ be a rooted tree with at least two vertices, and let

Γ 1 = τ •S 1 and Γ 2 = τ •S 2 .
The RACGs G(Γ 1 ) and G(Γ 2 ) have distinct geodesic growth series.

Proof. We will prove that f Geo(Γ 1 ) (10) = f Geo(Γ 2 ) [START_REF] Hopcroft | Introduction to automata theory, languages, and computation[END_REF]. We use the same argument as in the proof of Lemma

3.3, that is, f Geo(Γ 1 ) (r) = f Geo(Γ 2 ) (r) if and only if f Geo(Γ 1 ) (r) = f Geo(Γ 2 ) (r)
, where f Geo(Γ i ) (r) denotes the numbers of non-geodesics of length r, r ≥ 1. In this proof we use the result of our explicit computations with Monty, which shows that f Geo

0 (Γ 1 ) (r) = f Geo 0 (Γ 2 ) (r) for r < 8, but f Geo 0 (Γ 1 ) (8) = f Geo 0 (Γ 1 ) (8) = 8919523 for G 1 and f Geo 0 (Γ 2 ) (8) = f Geo 0 (Γ 2 ) (8) = 8919522 for G 2 .
Any word in Γ * 1 has the form w 1 u 1 w 2 . . . u n , where u i ∈ (τ \ {0}) * , w i ∈ S * 1 , and w i , u i non-empty except for perhaps w 1 and u n . A non-geodesic word w in Γ * 1 belongs to one of the following sets (or is of the type), depending on its form: (A) it either contains non-geodesics w i ∈ (τ \{0}) * or u j ∈ S * 1 , or both, for some 1 ≤ i, j ≤ n, or (B) all w i and u i are geodesic on their respective alphabets, and there exists 1 ≤ j < n such that 0u j 0 is a subword of w and u j ∈ (St Γ 1 (0) ∩ τ ) * , or (C) all w i and u i are geodesic on their respective alphabets, and w contains a subword of the form s0s, where s ∈ (St Γ 1 (0) ∩ τ ) * .

For the remaining discussion we only consider words w of length 10, and call the number of non-empty subwords u i or w i the syllable length of w. The number of words of type (A) depends only on the geodesic growth series of τ and S 1 , and thus it will be equal to the number of words of type (A) in Γ 2 . Notice that if the syllable length of w is ≤ 2 there are no words of type (B) or (C), so we obtain the same numbers of words in Γ 1 and Γ 2 . If the syllable length of w is > 4 then all w i , u i are shorter than or equal to 6, and our computations show that the numbers of special geodesics of length less than or equal to 6 in Γ 1 and Γ 2 coincide. Thus, a discrepancy may appear only when the syllable length is 3 or 4. The words of type (B) or (C) with syllable length 4 have the form w = w 1 u 1 w 2 u 2 (or w = u 1 w 2 u 2 w 3 ), where w 1 ends in 0 and w 2 starts with 0, or w 2 = 0. But in this case |w 1 |, |w 2 | ≤ 7, and our computations show that f Geo 0 (Γ 1 ) (r) = f Geo 0 (Γ 2 ) (r), up to r = 7.

Next, we consider the non-geodesics w of syllable length 3. If they have the form w = u 1 w 1 u 2 , then w 1 = 0, and the numbers of geodesics u i is the same since they are all written over the same alphabet determined by the tree τ , so no discrepancy in the numbers of non-geodesics occurs. Thus, it remains to count the words of the form w = w 1 u 1 w 2 , where |w 1 | ≤ 8, w 1 ends in 0, |u 1 | ≥ 1, u 1 ∈ (St Γ 1 (0) ∩ τ ) * , and w 2 starts with 0, |w 2 | ≤ 8. Again, if |w i | ≤ 7 we obtain the same numbers of non-geodesics. However, a discrepancy occurs when |w 1 | = 8 or |w 2 | = 8.

Indeed, the number of such words is 2(f Geo 0 (Γ 1 ) (8) deg τ (0)) in Γ * 1 and 2(f Geo 0 (Γ 2 ) (8) deg τ (0)) in Γ * 2 . By using Monty, we obtained f Geo 0 (Γ 1 ) (8) = f Geo 0 (Γ 1 ) (8) = 8919523 for G 1 , but f Geo 0 (Γ 2 ) (8) = f Geo 0 (Γ 2 ) (8) = 8919522 for G 2 , so the numbers of non-geodesics in these two groups are distinct. Now we can finish the proof of our second main result.

Proof. (of Theorem 1.2) This is a straightforward consequence of Theorem 5.1 and Lemma 5.2.

Figure 1 .Figure 2 .

 12 Figure 1. Trees τ and σ with marked roots (on the left) and their coalescence τ • σ (on the right)

Figure 3 .

 3 Figure 3. Godsil's rooted trees: S 1 on the left, and S 2 on the right

Figure 4 .

 4 Figure 4. The rooted tree σ with root 0

Acknowledgments

The authors gratefully acknowledge the support received from the Swiss National Science Foundation: L.C. was supported by PP00P2-144681 (SNSF Professorship); A.K. was supported by P300P2-151316 (Advanced Post-Doc Mobility) and P300P2-151316/2 (CH-Link). The authors are also thankful to Prof. Ruth Kellerhals (University of Fribourg, Switzerland), Prof. Michelle Bucher (University of Geneva, Switzerland), Prof. Robert Young (Courant Institute of Mathematical Sciences, New York, USA) and Dr. Alexey Talambutsa (Steklov Mathematical Institute, Moscow, Russia) for fruitful discussions.

Appendix

In this section we give the Python code "Monty" that we used in our computations, with comments and remarks. Its on-line copy [START_REF] Kolpakov | SAGE" worksheet, available on-line at[END_REF] can be downloaded as a SAGE worksheet from the second author's web-page.

We begin by defining the automaton A that recognizes the language of geodesics of a given RACG G whose defining graph Γ is given as input.

def Automaton(t): # takes a triangle-free graph, returns the automaton recognising # the resp. RACG as a digraph CliqueComplex = t.clique_complex(); C = list(); for s in CliqueComplex.faces().values(): for f in s: C.append(set(f)); n = len(C); a = DiGraph(); a.add_vertices(range(n)); for i in range(n):

for v in t.vertices(): Given the automaton A, we then compute the growth function of its accepted language (in this case, the geodesic language for G with respect to S).

denominator().factor();

A clique c in the defining graph Γ of the RACG G = G(Γ) corresponds to a state q c in the automaton A. We need the following auxiliary function in order to determine the index of q c represented as a vertex of the digraph A (the automaton) created by the procedure Automaton. Now we define a function that takes as input the geodesic automaton A for a RACG G = G(Γ), a list of cliques l = [c 0 , c 1 , . . . , c k ] in the respective defining graph Γ and returns the growth function for geodesic words starting at the state q = δ(q 0 , 0) that bring A to any of the states described by the cliques in l.

def GrowthFuncStart0End(a, l): # takes an automaton and a list of cliques 'l' in the resp. defining graph # as input, returns the growth function for geodesic words starting at # state q, where $\delta(Start, 0) = q$, and ending at any By using a suitable list of cliques l we can compute the functions 0 0 α(t) G and 0 0 β(t) G . Namely, in the proof of Theorem 1.1, we find # the function ${_0^0}\alpha(t)_{G_1}$ a001 = GrowthFuncStart0End(a1, [[0]]); # the function ${_0^0}\alpha(t)_{G_2}$ a002 = GrowthFuncStart0End(a2, [[0]]); # the function ${_0^0}\beta(t)_{G_1}$ b001 = GrowthFuncStart0End(a1, [ [START_REF] Antolín | Geodesic growth in right-angled and even Coxeter groups[END_REF], [START_REF] Antolín | Geodesic growth in right-angled and even Coxeter groups[END_REF][START_REF] Bjorner | Combinatorics of Coxeter groups[END_REF], [START_REF] Antolín | Geodesic growth in right-angled and even Coxeter groups[END_REF][START_REF] Charney | Geodesic automaton and growth functions for Artin groups of finite type[END_REF], [START_REF] Bahls | The isomorphism problem in Coxeter groups[END_REF], [START_REF] Bahls | The isomorphism problem in Coxeter groups[END_REF][START_REF] Davis | The geometry and topology of Coxeter groups[END_REF]

The list l = [[1], [START_REF] Antolín | Geodesic growth in right-angled and even Coxeter groups[END_REF][START_REF] Bjorner | Combinatorics of Coxeter groups[END_REF], [START_REF] Antolín | Geodesic growth in right-angled and even Coxeter groups[END_REF][START_REF] Charney | Geodesic automaton and growth functions for Artin groups of finite type[END_REF], [START_REF] Bahls | The isomorphism problem in Coxeter groups[END_REF], [START_REF] Bahls | The isomorphism problem in Coxeter groups[END_REF][START_REF] Davis | The geometry and topology of Coxeter groups[END_REF]] above contains the cliques of Γ 1 corresponding to the accept states q of the geodesic automaton A 1 for G 1 = G(Γ 1 ) such that δ(p, 0) = q, for a state p. The list l = [[1], [START_REF] Antolín | Geodesic growth in right-angled and even Coxeter groups[END_REF][START_REF] Bjorner | Combinatorics of Coxeter groups[END_REF], [START_REF] Bahls | The isomorphism problem in Coxeter groups[END_REF], [START_REF] Bahls | The isomorphism problem in Coxeter groups[END_REF][START_REF] Charney | Geodesic automaton and growth functions for Artin groups of finite type[END_REF], [START_REF] Bahls | The isomorphism problem in Coxeter groups[END_REF][START_REF] Davis | The geometry and topology of Coxeter groups[END_REF]] contains the cliques of Γ 2 with analogous properties, corresponding to the states of the geodesic automaton A 2 for the RACG G 2 = G(Γ 2 ).

The above described Python procedures are also used to perform the computations in the proof of Theorem 1.2.

The on-line version of Monty [START_REF] Kolpakov | SAGE" worksheet, available on-line at[END_REF] contains a variation of the GrowthFunc procedure, called GrowthFuncBM, that uses the Berlekamp-Massey algorithm for faster computing.