Daniele Antonio

Di Pietro

Martin Vohralík
email: martin.vohralik@inria.fr

Soleiman Yousef
email: yousef@ann.jussieu.fr

Daniele A Di Pietro

AND Martin Vohral Ík

AND H R A L Ík

Adaptive regularization, linearization, and discretization and a posteriori error

Keywords: Mathematics Subject Classification. Primary 65N08, 65N15, 65N50, 80A22

Introduction

The two-phase Stefan problem models a phase change process which is governed by the Fourier law; cf. Friedman [START_REF] Friedman | The Stefan problem in several space variables,T r a n s[END_REF]. The two phases, typically solid and liquid, are separated by a moving interface, whose motion is governed by the so-called Stefan condition. Let Ω ⊂ R d , d ∈{ 2, 3},b ea no p e nb o u n d e dp o l y g o n a lo rp o l y h e d r a l domain, not necessarily convex, and let T>0. The mathematical statement of the problem is as follows: given an initial enthalpy u 0 and a source function f , find the enthalpy u such that

∂ t u -∇•(∇β(u)) = f in Ω × (0,T), (1.1a) u(•, 0) = u 0 in Ω, (1.1b) β(u)=0
on∂Ω × (0,T). (1.1c) An inevitable tool in practical simulations seems to be an a posteriori error estimate-driven adaptive mesh refinement. One of the first works on a posteriori error estimates for the steady Stefan problem is that of Picasso [START_REF] Picasso | An adaptive finite element algorithm for a two-dimensional stationary Stefanlike problem,C o m p u t .M e t h o d sA p p l[END_REF]. Therein, the author derives estimates based on the residual equation for a two-dimensional regularized Stefan problem and proposes a space adaptive finite element algorithm. A posteriori indicators for unsteady phase change problems were derived by Nochetto et al. in [START_REF] Nochetto | An adaptive finite element method for two-phase Stefan problems in two space dimensions. I. Stability and error estimates[END_REF][START_REF] Nochetto | An adaptive finite element method for twophase Stefan problems in two space dimensions. II. Implementation and numerical experiments[END_REF], together with an adaptive algorithm which equilibrates space and time discretization errors. Many other adaptive refinement algorithms such as that of Beckett et al. [START_REF] Beckett | Am o v i n gm e s hfi n i t ee l e m e n tm e t h o d for the solution of two-dimensional Stefan problems[END_REF] have also been proposed. Rigorous a posteriori error estimates for nonlinear parabolic problems seem much less developed. In nondegenerate cases, Verfürth [START_REF] Verfürth | Ap o s t e r i o r ie r r o re s t i m a t e sf o rn o n l i n e a rp r o b l e m s .L r (0,T; L ρ (Ω))-error estimates for finite element discretizations of parabolic equations[END_REF][START_REF] Verfürth | r i o r ie r r o re s t i m a t e sf o rn o n l i n e a rp r o b l e m s :L r (0,T; W 1,ρ (Ω))error estimates for finite element discretizations of parabolic equations,N u m[END_REF] was able to obtain an estimator which is both reliable and efficient. A pioneering contribution for degenerate parabolic problems has been obtained by Nochetto et al. in [START_REF] Nochetto | i o r ie r r o re s t i m a t i o na n da d a p t i v i t yf o r degenerate parabolic problems[END_REF]. Therein, L ∞ (0,T; H -1 (Ω)) estimates for the error in the enthalpy and L 2 (0,T; L 2 (Ω)) estimates for the error in the temperature are obtained. The approach is based on the relation of these errors to the residual of (1.1a) obtained through the corresponding dual partial differential equation and subsequent use of the Galerkin orthogonality of the finite element discretization. Recently, rigorous a posteriori error analysis in a space-time dual norm, including some degenerate cases, was given in [START_REF] Dolejš´ı | Af r a m e w o r kf o rr o b u s tap o s t e r i o r ie r r o rc o n t r o li n unsteady nonlinear advection-diffusion problems[END_REF].

The aim of this paper is to derive fully computable a posteriori error estimates and adaptive strategies for the two-phase Stefan problem (1.1) for conforming spatial discretization schemes such as finite element, co-volume, or vertexcentered finite volume schemes with backward Euler time stepping. As in Nochetto et al. [START_REF] Nochetto | i o r ie r r o re s t i m a t i o na n da d a p t i v i t yf o r degenerate parabolic problems[END_REF], our approach is based on the dual norm of the residual. However, we proceed differently in order to have fully and easily computable estimates not featuring any undetermined constants. This is achieved by introducing H(div; Ω)-conforming and locally conservative flux reconstructions following Prager and Synge [START_REF] Prager | Approximations in elasticity based on the concept of function space[END_REF], Ladevèze [START_REF] Ladevèze | Comparaison de modèles de milieux continus.P h . D .t h e s i s[END_REF], Destuynder and Métivet [START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF], Luce and Wohlmuth [START_REF] Luce | a lap o s t e r i o r ie r r o re s t i m a t o rb a s e do ne q u i l ibrated fluxes[END_REF], Braess and Schöberl [START_REF] Braess | Equilibrated residual error estimator for edge elements[END_REF], Repin [START_REF] Repin | AP o s t e r i o r iE s t i m a t e sf o rP a r t i a lD i ff e r e n t i a lE q u a t i o n s ,R a d o nS e r i e so nC o mputational and Applied Mathematics[END_REF], and [START_REF] Dolejš´ı | Af r a m e w o r kf o rr o b u s tap o s t e r i o r ie r r o rc o n t r o li n unsteady nonlinear advection-diffusion problems[END_REF][START_REF] Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF][START_REF] Ern | Ap o s t e r i o r ie r r o re s t i m a t i o nb a s e do np o t e n t i a la n dfl u xr e c o nstruction for the heat equation[END_REF][START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF][START_REF] Vohralík | Ap o s t e r i o r ie r r o re s t i m a t i o ni nt h ec o n f o r m i n gfi n i t ee l e m e n tm e t h o d based on its local conservativity and using local minimization (English, with English and French summaries)[END_REF]; see also the references therein.

In Section 2 we give a weak formulation, introduce a regularized problem with ar e g u l a r i z a t i o np a r a m e t e rϵ>0, and fix the notation for temporal and spatial meshes. In Section 3, we identify the residual and its dual norm and we derive an a posteriori error estimate on this problem-dependent error measure. We next split this estimate into estimators characterizing the space, time, regularization, linearization, and quadrature errors.

Section 4 subsequently presents a criterion for the choice of the regularization parameter ϵ and a stopping criterion for an iterative linearization such as the Newton method. The former is designed to facilitate the treatment of the degeneracy while not spoiling the accuracy, whereas the latter is designed to avoid performing an excessive number of nonlinear solver iterations. These criteria are inspired mainly from [START_REF] Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF][START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF][START_REF] Jiránek | Ap o s t er i o r ier r o re s t i m a t esi n c l ud i n ga l g e b r a i cer r o r and stopping criteria for iterative solvers[END_REF]. We then propose an adaptive algorithm which uses these criteria while simultaneously performing the usual local mesh refinement and equilibration of the spatial and temporal errors. This algorithm is inspired from [START_REF] Beckett | Am o v i n gm e s hfi n i t ee l e m e n tm e t h o d for the solution of two-dimensional Stefan problems[END_REF][START_REF] Nochetto | An adaptive finite element method for two-phase Stefan problems in two space dimensions. I. Stability and error estimates[END_REF][START_REF] Nochetto | An adaptive finite element method for twophase Stefan problems in two space dimensions. II. Implementation and numerical experiments[END_REF][START_REF] Nochetto | i o r ie r r o re s t i m a t i o na n da d a p t i v i t yf o r degenerate parabolic problems[END_REF][START_REF] Picasso | An adaptive finite element algorithm for a two-dimensional stationary Stefanlike problem,C o m p u t .M e t h o d sA p p l[END_REF] and from the work [START_REF] Dolejš´ı | Af r a m e w o r kf o rr o b u s tap o s t e r i o r ie r r o rc o n t r o li n unsteady nonlinear advection-diffusion problems[END_REF][START_REF] Ern | Ap o s t e r i o r ie r r o re s t i m a t i o nb a s e do np o t e n t i a la n dfl u xr e c o nstruction for the heat equation[END_REF][START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF][START_REF] Hilhorst | Ap o s t e r i o r ie r r o re s t i m a t e sf o rc o m b i n e dfi n i t ev o l u m efinite element discretizations of reactive transport equations on nonmatching grids[END_REF]. We conclude Section 4 by proving that, under these criteria, our estimators are also efficient while representing a lower bound for the dual norm of the residual.

In Section 5, we show how to bound the L 2 (0,T; H -1 (Ω))-type error in the enthalpy and L 2 (0,T; L 2 (Ω))-type error in the temperature by the above dual norm of the residual. In particular, we focus on the use of the Gronwall lemma with as small an overestimation as possible and no appearance of the exponential term e T other than in the approximation of the initial condition. Guaranteed and fully computable a posteriori error estimates on these natural norms immediately follow.

Section 6 presents the application of all these developments to the vertex-centered finite volume (or, equivalently, finite element with mass lumping and numerical quadrature) discretization in space, backward Euler discretization in time, and Newton linearization. Illustrative numerical results fill up Section 7 and, finally, Appendix A collects the more involved proofs of the various theorems of the paper.

Continuous and discrete settings

This section fixes the basic continuous and discrete settings. More precisely, Section 2.1 presents the continuous problem and the regularization, whereas the basic assumptions on the discretization are introduced in Section 2.2.

2.1. Continuous setting.

2.1.1.

The continuous problem. The starting point for our a posteriori analysis is the weak form of problem (1.1). To give it, we need to introduce the assumptions on the data and set up some notation. We suppose that: (i) the enthalpy-temperature function β : R → R is a Lipschitz continuous function such that β(s) = 0 in (0, 1), β is strictly increasing in R -and R + \ (0, 1), and there exist c, C > 0s u c ht h a t , for all s ∈ R \ (0, 1), sign(s)β(s) ≥ c|s|-C (see Figure 1); the Lipschitz constant of β is denoted by L β ; (ii) the initial enthalpy u 0 is such that u 0 ∈ L 2 (Ω); (iii) the source term is such that f ∈ L 2 (0,T; L 2 (Ω)). We will rep eatedly use throughout the pap er the two following spaces:

(2.1) X := L 2 (0,T; H 1 0 (Ω)),Z := H 1 (0,T; H -1 (Ω)). We will also need the dual space X ′ of X,

X ′ = L 2 (0,T; H -1 (Ω)),
and equip the space X with the norm

∥ϕ∥ X := T 0 ∥∇ϕ(•,t)∥ 2 L 2 (Ω) dt 1 2
.

We denote by ⟨•, •⟩ the duality pairing between H -1 (Ω) and

H 1 0 (Ω), while (•, •) S is the usual scalar product in L 2 (S)o r[ L 2 (S)] d ,

w i t ht h es u b s c r i p to m i t t e dw h e n S =Ω.

The weak formulation of problem (1.1) can now be stated. It reads: find

(2.2a) u ∈ Z with β(u) ∈ X such that (2.2b) u(•, 0) = u 0 in Ω
and, for a.e. s ∈ (0,T),

(2.2c) ⟨∂ t u(•,s),ϕ⟩ +(∇β(u(•,s)), ∇ϕ)=(f (•,s),ϕ) ∀ϕ ∈ H 1 0 (Ω).
I(t) := {x ∈ Ω:β(u)(x,t)=0} .
This fact may hinder both the design and the convergence analysis of a discretization method. Additionally, the lack of smoothness in the dependency of the solution on the problem data can severely affect the convergence of nonlinear iterations. A possible and often employed approach [START_REF] Beckett | Am o v i n gm e s hfi n i t ee l e m e n tm e t h o d for the solution of two-dimensional Stefan problems[END_REF][START_REF] Nochetto | Error estimates for multidimensional singular parabolic problems,J a p a nJ[END_REF][START_REF] Nochetto | The combined use of a nonlinear Chernoff formula with a regularization procedure for two-phase Stefan problems[END_REF][START_REF] Picasso | An adaptive finite element algorithm for a two-dimensional stationary Stefanlike problem,C o m p u t .M e t h o d sA p p l[END_REF] to overcome these difficulties consists in regularizing the problem (2.2) by replacing the function β by a smooth, strictly increasing regularized function

β ϵ ∈ C 1 (R), β ′ ϵ ≥ ϵ,f o
rap a r a m e t e rϵ>0; see Figure 1 for an example.

The regularized problem reads as follows: find

(2.3a) u ϵ ∈ Z with β ϵ (u ϵ ) ∈ X such that (2.3b) u ϵ (•, 0) = β -1 ϵ (β(u 0 )) in Ω,
and, for a.e. s ∈ (0,T),

(2.3c) ⟨∂ t u ϵ (•,s),ϕ⟩ +(∇β ϵ (u ϵ (•,s)), ∇ϕ)=(f (•,s),ϕ) ∀ϕ ∈ H 1 0 (Ω).

Discrete setting.

We describ e here the basic discrete setting that will b e sufficient for the developments of Sections 3-5. Further details are given in Section 6.

Time mesh.

Our focus is on first-order time discretizations based on the backward Euler scheme. Let {τ n } 1≤n≤N denote a sequence of positive real numbers corresponding to the discrete time steps such that T = N n=1 τ n . We let t 0 :=0 and, for 1 ≤ n ≤ N , we introduce the discrete times t n := n i=1 τ i and the time intervals License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

I n :=(t n-1 ,t n ).
6D A N I E L E A . D I P I E T R O , M A R T I N V O H R A L ÍK, AND SOLEIMAN YOUSEF

Space meshes.

Let {K n } 0≤n≤N denote a family of matching simplicial meshes of the space domain Ω. The initial mesh K 0 is used to approximate the initial condition u 0 ,w h i l eK n is used to march in time from t n-1 to t n . The meshes can be refined or coarsened as time evolves. For the developments of Section 4.3 below, we are led to suppose that K n ,1≤ n ≤ N , is obtained from K n-1 by refining some elements and coarsening (a limited number of times) some other ones. We denote by K n-1,n the coarsest common submesh (overlay) of both K n and K n-1 and, once again for the developments of Section 4.3, suppose that the meshes {K n-1,n } 1≤n≤N are shape-regular in the sense that there exists a constant κ K > 0s u c ht h a t min

K∈K n-1,n ρ K h K ≥ κ K
for all 1 ≤ n ≤ N ,w h e r eρ K denotes the diameter of the largest ball inscribed in the element K and h K the diameter of K.F o r0≤ n ≤ N ,w ed e n o t eb yΠ n 0 the L 2 -orthogonal projection onto the space of piecewise constant functions on K n .

For 0 ≤ n ≤ N , let F n denote the set of mesh faces. Boundary faces are collected in the set F b,n := {F ∈F n ; F ⊂ ∂Ω} and we let F i,n := F n \F b,n . For a given face F ∈F i,n we fix an arbitrary orientation and denote the corresponding unit normal vector by n F ;f o rF ∈F b,n , n F coincides with the exterior unit normal n Ω of Ω. A similar notation for the faces F n-1,n of the meshes K n-1,n will also be used.

An a posteriori error estimate for the dual norm of the residual

In this section we derive an a posteriori estimate for the error measured by the dual norm of the residual that we first identify. We then give a basic estimate that we subsequently refine to distinguish the space, time, linearization, regularization, and quadrature errors.

3.1.

Dual norm of the residual. As in Picasso [START_REF] Picasso | An adaptive finite element algorithm for a two-dimensional stationary Stefanlike problem,C o m p u t .M e t h o d sA p p l[END_REF] or Nochetto et al. [START_REF] Nochetto | i o r ie r r o re s t i m a t i o na n da d a p t i v i t yf o r degenerate parabolic problems[END_REF], our key for deriving a posteriori error estimates for the Stefan problem (1.1) will be the residual and its dual norm. Recall that u denotes the weak solution of the Stefan problem given by (2.2) and the definition of the space X (2.1). Let u hτ ∈ Z such that β(u hτ ) ∈ X is arbitrary. In practice, u hτ will be the result of the numerical simulation. We define the residual R(u hτ ) ∈ X ′ such that

⟨R(u hτ ),ϕ⟩ X ′ ,X (3.1) := T 0 {⟨∂ t (u -u hτ ),ϕ⟩ +(∇β(u) -∇β(u hτ ), ∇ϕ)} (s)ds, ϕ ∈ X.
Using (2.2c), we can infer the following alternative expression for (3.1):

⟨R(u hτ ),ϕ⟩ X ′ ,X = T 0 {(f, ϕ) -⟨∂ t u hτ ,ϕ⟩-(∇β(u hτ ), ∇ϕ)} (s)ds, ϕ ∈ X.
The norm of the residual in the dual space X ′ is then given by

(3.2) ∥R(u hτ )∥ X ′ :=s u p ϕ∈X, ∥ϕ∥ X =1 ⟨R(u hτ ),ϕ⟩ X ′ ,X .
The key problem-specific measure of the distance between u hτ and u that we will use in this paper is given by

(3.3) ∥R(u hτ )∥ X ′ + ∥u 0 -u hτ (•, 0)∥ H -1 (Ω) .
It follows from (2.2) that the measure (3.3) is zero if and only if the function u hτ coincides with the exact solution u. As we shall see below in Section 5, it in fact controls the energy error between u and u hτ and β(u)a n dβ(u hτ ). Moreover, this quantity can be easily bounded in terms of error estimators based on H(div; Ω)conforming flux reconstructions for piecewise affine-in-time u hτ that we show next.

3.2. General assumptions. In order to proceed with the analysis further, without the necessity to specify at this point any details on how the approximate solution u hτ was obtained, we are lead to make the following assumption. It requires Z-and X-conformity and u hτ to be piecewise affine and continuous in time on the time mesh {I n } 1≤n≤N of Section 2.2.1:

Assumption 3.1 (Approximate solution).
The function u hτ is such that

u hτ ∈ Z, ∂ t u hτ ∈ L 2 (0,T; L 2 (Ω)),β (u hτ ) ∈ X, u hτ | I n is affine in time on I n ∀1 ≤ n ≤ N.
Note that, consequently, the function u hτ is uniquely determined by the N +1 functions

u n h := u hτ (•,t n ), 0 ≤ n ≤ N ,a n d∂ t u n hτ := ∂ t u hτ | I n =( u n h -u n-1 h
)/τ n , 1 ≤ n ≤ N .W ew i l la l s oe m p l o yt h ea b r i d g e dn o t a t i o nu n hτ for u hτ | I n . The second assumption that we make is the existence of a piecewise constantin-time H(div; Ω)-conforming flux reconstruction t hτ , locally conservative on the meshes K n of Section 2.2.2. Let us first denote by f the piecewise constant-in-time function given by the time-mean values of the source function f on the intervals I n ,1≤ n ≤ N . Assumption 3.2 (Equilibrated flux reconstruction). For al l 1 ≤ n ≤ N ,t h e r e exists a vector field t n h ∈ H(div; Ω) such that

(∇•t n h , 1) K =( f n , 1) K -(∂ t u n hτ , 1) K ∀K ∈K n .
We denote by t hτ the space-time function such that t hτ | I n := t n h for all 1 ≤ n ≤ N . In Section 6 below, we show how to construct an equilibrated flux reconstruction t hτ in the context of vertex-centered finite volume (finite element with mass lumping and quadrature) spatial discretization.

3.3.

Ab a s i cap o s t e r i o r ie r r o re s t i m a t e .We now give an a p osteriori error estimate in the general setting of Assumptions 3.1 and 3.2. Note that the regularization of Section 2.1.2 is not used at the present stage.

We will estimate the error measure (3.3) by the lo cal residual expressed with the flux t hτ and by the difference of t hτ and the temperature flux, in the spirit of [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] Ladevèze | Comparaison de modèles de milieux continus.P h . D .t h e s i s[END_REF][START_REF] Luce | a lap o s t e r i o r ie r r o re s t i m a t o rb a s e do ne q u i l ibrated fluxes[END_REF][START_REF] Prager | Approximations in elasticity based on the concept of function space[END_REF][START_REF] Repin | AP o s t e r i o r iE s t i m a t e sf o rP a r t i a lD i ff e r e n t i a lE q u a t i o n s ,R a d o nS e r i e so nC o mputational and Applied Mathematics[END_REF] and [START_REF] Dolejš´ı | Af r a m e w o r kf o rr o b u s tap o s t e r i o r ie r r o rc o n t r o li n unsteady nonlinear advection-diffusion problems[END_REF][START_REF] Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF][START_REF] Ern | Ap o s t e r i o r ie r r o re s t i m a t i o nb a s e do np o t e n t i a la n dfl u xr e c o nstruction for the heat equation[END_REF][START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF][START_REF] Vohralík | Ap o s t e r i o r ie r r o re s t i m a t i o ni nt h ec o n f o r m i n gfi n i t ee l e m e n tm e t h o d based on its local conservativity and using local minimization (English, with English and French summaries)[END_REF]. More specifically, for 1 ≤ n ≤ N , t n h as in Assumption 3.2, and K ∈K n ,w ed e fi n et h eresidual estimator η n R,K and the flux estimator η n F,K as follows:

η n R,K := C P,K h K ∥ f n -∂ t u n hτ -∇•t n h ∥ L 2 (K) , (3.4a) η n F,K (t) := ∥t n h + ∇β(u hτ (•,t))∥ L 2 (K) t ∈ I n . (3.4b)
Here, C P,K is the constant from the Poincaré inequality

(3.5) ∥ϕ -Π n 0 ϕ∥ L 2 (K) ≤ C P,K h K ∥∇ϕ∥ L 2 (K) ∀ϕ ∈ H 1 (K).
It holds that C P,K =1 /π as the simplices K are convex; see [START_REF] Bebendorf | An o t eo nt h eP o i n c a r éi n e q u a l i t yf o rc o n v e xd o m a i n s[END_REF][START_REF] Payne | An optimal Poincaréine qualityforc onvexdomains[END_REF]. Finally, we define the initial condition estimator by

(3.6) η IC := ∥u 0 -u hτ (•, 0)∥ H -1 (Ω) .
We then have:

Theorem 3.3 (A posteriori estimate for the error measure (3.3)). Let u be the weak solution given by (2.2) and let u hτ and t hτ fulfill Assumptions 3.1 and 3.2, respectively. Then, there holds

(3.7) ∥R(u hτ )∥ X ′ + ∥u 0 -u hτ (•, 0)∥ H -1 (Ω) ≤ η + η IC ,
where

(3.8) η := N n=1 I n K∈K n η n R,K + η n F,K (t) 2 dt 1 2 + ∥f -f ∥ X ′ .
Proof. Let ϕ ∈ X with ∥ϕ∥ X = 1 be given. Then, adding and subtracting (t hτ , ∇ϕ) and using Green's theorem, it holds that

⟨R(u hτ ),ϕ⟩ X ′ ,X = T 0 {(f -∂ t u hτ -∇•t hτ ,ϕ) -(t hτ + ∇β(u hτ ), ∇ϕ)} (s)ds = T 0 (f -f,ϕ)+( f -∂ t u hτ -∇•t hτ ,ϕ) -(t hτ + ∇β(u hτ ), ∇ϕ) (s)ds = : T 1 + T 2 + T 3 .
For the first term we infer

T 1 ≤∥ f -f ∥ X ′ ∥ϕ∥ X = ∥f -f ∥ X ′ .
T h es e c o n dt e r m can be rewritten as follows:

T 2 = N n=1 I n ( f n -∂ t u n hτ -∇•t n h ,ϕ)(s)ds.
For all 1 ≤ n ≤ N and t ∈ I n ,thereholds(thedependenceofϕ on the time variable is omitted for brevity), where we have used the regularity of the arguments, Assumption 3.2, the Cauchy-Schwarz inequality, and the Poincaré inequality (3.5). For the third term, an application of the Cauchy-Schwarz inequality yields

( f n -∂ t u n hτ -∇•t n h ,ϕ)= K∈K n ( f n -∂ t u n hτ -∇•t n h ,ϕ) K = K∈K n ( f n -∂ t u n hτ -∇•t n h ,ϕ-Π n 0 ϕ) K ≤ K∈K n ∥ f n -∂ t u n hτ -∇•t n h ∥ L 2 (K) ∥ϕ -Π n 0 ϕ∥ L 2 (K) ≤ K∈K n C P,K h K ∥ f n -∂ t u n hτ -∇•t n h ∥ L 2 (K) ∥∇ϕ∥ L 2 (K) = K∈K n η n R,K ∥∇ϕ∥ L 2 (K) ,
T 3 ≤ N n=1 I n K∈K n η n F,K ∥∇ϕ∥ L 2 (K) (s)ds.
Collecting the above estimates, using the definition (3.2) of the dual norm of the residual, and using the Cauchy-Schwarz inequality yields (3.7).

3.4. An a posteriori error estimate distinguishing the space, time, regularization, linearization, and quadrature errors. Our next goal is to distinguish the different error components. This is an instrumental step to design an adaptive algorithm where the time step, the space mesh, the regularization parameter, and the stopping criterion for the linearization iterations are chosen optimally. We start by localizing in time the error measure introduced in Section 3.1. For 1 ≤ n ≤ N , we let

X n := L 2 (I n ; H 1 0 (Ω)),Z n := H 1 (I n ; H -1 (Ω)).
We localize in time the dual norm of the residual (3.2) by setting

∥R(u hτ )∥ X ′ n (3.9) :=s u p ϕ∈X n , ∥ϕ∥ X n =1 I n {⟨∂ t (u -u hτ ),ϕ⟩ +(∇β(u) -∇β(u hτ ), ∇ϕ)} (s)ds.
Note that, consequently,

∥R(u hτ )∥ 2 X ′ = N n=1 ∥R(u hτ )∥ 2 X ′ n for any u hτ ∈ Z with β(u hτ ) ∈ X.
Suppose now that we are marching in time from time t n-1 to time t n with a given time step τ n ,s t a r t i n gf r o mt h ea p p r o x i m a t i o nu n-1 h .W ea l s os u p p o s et h a t the regularization of Section 2.1.2 has been used for a given value of the parameter ϵ,a n dt h a tw ea r eo nt h ek-th step of some iterative linearization algorithm. We denote by u n,ϵ,k h the approximation of the solution u at time t n and prescribe the space-time function u n,ϵ,k hτ by the value

u n-1 h at time t n-1 ,b yt h ev a l u eu n,ϵ,k h at time t n ,a n db ya ffi n eb e h a v i o ri nt i m eo nI n ,i . e . , (3.10) u n,ϵ,k hτ (•,t)=(1-ρ(t))u n-1 h + ρ(t)u n,ϵ,k h ,ρ (t) := t -t n-1
τ n . We summarize our general requirements in the following: Assumption 3.4 (Adaptive setting). For al l 1 ≤ n ≤ N , a regularization parameter ϵ ≥ 0, and a linearization step k ≥ 1:

(i) u n,ϵ,k hτ is the approximate solution given by (3.10), u n,ϵ,k hτ ∈ Z n with ∂ t u n,ϵ,k hτ ∈ L 2 (I n ; L 2 (Ω)) and β(u n,ϵ,k hτ ) ∈ X n ; (ii) there exists an equilibrated flux t n,ϵ,k h ∈ H(div; Ω) such that (3.11) (∇•t n,ϵ,k h , 1) K =( f n , 1) K -(∂ t u n,ϵ,k hτ ,
(iii) l n,ϵ,k h ∈ [L 2 (Ω)] d is the available approximation of the flux ∇β ϵ (u(•,t n )); (iv) Π n is

an operator used for interpolatory numerical integration.

An example of the approximate solution u n,ϵ,k hτ ,t h el i n e a r i z e dfl u xl n,ϵ,k h ,a n dth e operator Π n in the context of the implicit vertex-centered finite volume discretization and Newton linearization is provided in Section 6.3 below.

Proceeding as in Theorem 3.3, it is immediately inferred that

(3.12) ∥R(u n,ϵ,k hτ )∥ X ′ n ≤ I n K∈K n η n,ϵ,k R,K + η n,ϵ,k F,K (t) 2 dt 1 2 + ∥f -f ∥ X ′ n ,
where

η n,ϵ,k R,K := C P,K h K ∥ f n -∂ t u n,ϵ,k hτ -∇•t n,ϵ,k h ∥ L 2 (K) , η n,ϵ,k F,K (t) := ∥t n,ϵ,k h + ∇β(u n,ϵ,k hτ (•,t))∥ L 2 (K) ,t ∈ I n .
For all K ∈K n , we next define the local spatial, temporal, quadrature, regularization, and linearization estimators as follows:

η n,ϵ,k sp,K := η n,ϵ,k R,K + ∥l n,ϵ,k h + t n,ϵ,k h ∥ L 2 (K) , (3.13a) η n,ϵ,k tm,K (t) := ∥∇(Π n β(u n,ϵ,k hτ (•,t))) -∇(Π n β(u n,ϵ,k h ))∥ L 2 (K) ,t ∈ I n , (3.13b) η n,ϵ,k qd,K (t) := ∥∇(β(u n,ϵ,k hτ (•,t))) -∇(Π n β(u n,ϵ,k hτ (•,t)))∥ L 2 (K) ,t ∈ I n , (3.13c) η n,ϵ,k reg,K := ∥∇(Π n β(u n,ϵ,k h )) -∇(Π n β ϵ (u n,ϵ,k h ))∥ L 2 (K) , (3.13d) η n,ϵ,k lin,K := ∥∇(Π n β ϵ (u n,ϵ,k h )) -l n,ϵ,k h ∥ L 2 (K) . (3.13e)
Global versions of these estimators are given by (η n,ϵ,k sp

) 2 := τ n K∈K n η n,ϵ,k sp,K 2 , (3.14a) (η n,ϵ,k tm ) 2 := I n K∈K n η n,ϵ,k tm,K (t) 2 dt, (3.14b) (η n,ϵ,k qd ) 2 := I n K∈K n η n,ϵ,k qd,K (t) 2 dt, (3.14c) (η n,ϵ,k reg ) 2 := τ n K∈K n η n,ϵ,k reg,K 2 , (3.14d) (η n,ϵ,k lin ) 2 := τ n K∈K n η n,ϵ,k lin,K 2 . (3.14e)
Using the inequality (3.12) followed by the triangle inequality we obtain the following estimate: Corollary 3.5 (Distinguishing the space, time, quadrature, regularization, linearization, and data oscillation errors). Let u be the weak solution given by (2.2),

let 1 ≤ n ≤ N , ϵ ≥ 0,a n dk ≥ 1,a n dl e tu n,ϵ,k hτ , t n,ϵ,k h , l n,ϵ,k h
,a n dΠ n be as described in Assumption 3.

4.T h e nt h e r eh o l d s

∥R(u n,ϵ,k hτ )∥ X ′ n ≤ η n,ϵ,k sp + η n,ϵ,k tm + η n,ϵ,k qd + η n,ϵ,k reg + η n,ϵ,k lin + ∥f -f ∥ X ′ n .
Licensed License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use Remark 3.6 (Time oscillation of the source term). The error due to the time oscillation of the source term ∥f -f ∥ X ′ n is zero provided that the source function f is piecewise constant in time.

Balancing and stopping criteria, adaptive algorithm, and efficiency

The individual error component estimators of Corollary 3.5 are used in this section to define adaptive criteria to stop the iterative linearizations, to select the value of the regularization parameter ϵ, to locally adapt the quadrature rule, to adjust the time step, and to select the mesh elements to refine/derefine. These criteria are incorporated in a fully adaptive algorithm detailed in Section 4.2. Finally, in Section 4.3 we show the efficiency of our estimators when the adaptive balancing and stopping criteria are used.

4.1.

Balancing and stopping criteria. Following [START_REF] Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF][START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF][START_REF] Jiránek | Ap o s t er i o r ier r o re s t i m a t esi n c l ud i n ga l g e b r a i cer r o r and stopping criteria for iterative solvers[END_REF], this section introduces stopping criteria for the iterative algorithms based on the estimators of Corollary 3.5. The goal is to stop the iterations as soon as the corresponding error component no longer affects significantly the overall error. We assume in what follows that we are marching in time from time t n-1 to time t n .L e tt h r e eu s e r -g i v e n parameters Γ lin ,Γ reg ,Γ qd ∈ (0, 1) be given. The criteria are:

(i) Linearization. The linearization iteration is pursued until step k n such that

(4.1) η n,ϵ,k n lin ≤ Γ lin η n,ϵ,k n sp + η n,ϵ,k n tm + η n,ϵ,k n qd + η n,ϵ,k n reg .
(ii) Regularization. The regularization parameter ϵ is reduced until the value ϵ n such that

(4.2) η n,ϵ n ,k n reg ≤ Γ reg η n,ϵ n ,k n sp + η n,ϵ n ,k n tm + η n,ϵ n ,k n qd .
(iii) Quadrature. The quadrature rule is improved until

(4.3) η n,ϵ n ,k n qd ≤ Γ qd η n,ϵ n ,k n sp + η n,ϵ n ,k n tm .
Note that all the linearization, regularization, and quadrature errors may be classified as subsidiary as they can be made as small as desired by increasing the computational effort for fixed mesh and time step; it is thus reasonable to expect that the above criteria will be attained. Local, element by element, versions of the criteria (4.1)-( 4.3) can be formulated using the local estimators (3.13) (see [START_REF] Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF][START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF][START_REF] Jiránek | Ap o s t er i o r ier r o re s t i m a t esi n c l ud i n ga l g e b r a i cer r o r and stopping criteria for iterative solvers[END_REF]), and require that the inequalities hold for all K ∈K n ;c f .( 7 . 4 )f o ra ne x a m p l e .

In the spirit of [START_REF] Nochetto | i o r ie r r o re s t i m a t i o na n da d a p t i v i t yf o r degenerate parabolic problems[END_REF][START_REF] Picasso | Adaptive finite elements for a linear parabolic problem,C o m p u t[END_REF][START_REF] Verfürth | i o r ie r r o re s t i m a t e sf o rfi n i t ee l e m e n td i s c e t i z a t i o n so ft h eh e a t equation[END_REF] and [START_REF] Dolejš´ı | Af r a m e w o r kf o rr o b u s tap o s t e r i o r ie r r o rc o n t r o li n unsteady nonlinear advection-diffusion problems[END_REF][START_REF] Ern | Ap o s t e r i o r ie r r o re s t i m a t i o nb a s e do np o t e n t i a la n dfl u xr e c o nstruction for the heat equation[END_REF][START_REF] Hilhorst | Ap o s t e r i o r ie r r o re s t i m a t e sf o rc o m b i n e dfi n i t ev o l u m efinite element discretizations of reactive transport equations on nonmatching grids[END_REF], we also propose the usual space-time adaptivity: (iv) Space-time error balancing. The space and time error components should be equilibrated by selecting the time step τ n and adjusting the spatial meshes K n in such a way that (4.4)

γ tm η n,ϵ n ,k n sp ≤ η n,ϵ n ,k n tm ≤ Γ tm η n,ϵ n ,k n sp .
Above, Γ tm >γ tm > 0 are again user-given parameters, typically close to 1. (v) Adaptive mesh refinement. The error in space should be evenly distributed throughout the domain Ω by local adaptation (refinement, coarsening) of the space mesh K n in such a way that, for all

K 1 ,K 2 ∈K n , η n,ϵ n ,k n sp,K 1 ≈ η n,ϵ n ,k n sp,K 2 .
Licensed In contrast to (4.1)-(4.3), the goal is to make η n,ϵ n ,k n sp and η n,ϵ n ,k n tm of comparable size as these error components are substantial and cannot be made arbitrarily small for a given choice of the mesh and of the time step. 4.2. Adaptive algorithm. In this section we propose an adaptive algorithm that implements the balancing and stopping criteria of Section 4.1. Moreover, for a prescribed ζ>0, we aim at satisfying the relation

(4.5) N n=1 ∥R(u hτ )∥ 2 X ′ n N n=1 ∥l n,ϵ,k h ∥ 2 L 2 (I n ;L 2 (Ω)) ≤ ζ 2 ,
i.e., to bring the relative error under the user-given precision ζ.T o a c c o u n tf o r limited computing resources, we fix refinement thresholds h,τ > 0f o rb o t ht h e mesh size and the time step and require, for all 0 ≤ n ≤ N , (4.6) min K∈K n h K ≥ h,τ n ≥ τ . Note that, in particular, because of (4.6), the attainment of (4.5) is not guaranteed.

Recall that u n,ϵ,k h stands for the approximation of the solution u n h at discrete time t n obtained after k linearization iterations using a regularization parameter ϵ. At each linearization iteration k,t h en e wa p p r o x i m a t i o nu n,ϵ,k h is obtained solving the linear problem written schematically as

u n,ϵ,k h =Ψ ( u n,ϵ,k-1 h ,τ n , K n ).
For the sake of simplicity, in what follows we neglect the quadrature and data oscillation estimators. Our adaptive algorithm is the following:

Algorithm 4.1 (Adaptive algorithm).
Fix the fractions of cells to refine, c ref ,a n dt od e r e fi n e ,c deref Choose an initial mesh K 0 ,r e g u l a r i z a t i o np a r a m e t e rϵ 0 ,a n dat o l e r a n c eζ IC > 0 u 0 h ← Π 0 (β -1 ϵ 0 (β(u 0 ))) repeat {Initial mesh and regularization parameter adaptation} Compute η IC Refine the cells K ∈K 0 such that η IC,K ≥ c ref max L∈K 0 η IC,L in accordance with (4.6) and adjust the regularization parameter

ϵ 0 u 0 h ← Π 0 (β -1 ϵ 0 (β(u 0 ))) until η IC ≤ ζ IC ∥∇(β ϵ 0 (u 0 h ))∥ L 2 (Ω) Choose an initial time step τ 0 ϵ ← ϵ 0 , t 0 ← 0, n ← 0 while t n ≤ T do {Time loop} n ← n +1 K n ←K n-1 τ n ← τ n-1 u n,ϵ,0 h ← u n-1 h repeat {Space refinement} repeat {Space and time error balancing} repeat {Regularization} k ← 0 repeat {Nonlinear solver} k ← k +1 u n,ϵ,k h =Ψ(u n,ϵ,k-1 h ,τ n , K n ) Compute η n,ϵ,k sp , η n,ϵ,k tm , η n,ϵ,k reg , η n,
k n ← k if (4.2) does not hold then ϵ ← ϵ/2 end if until (4.2) is satisfied ϵ n ← ϵ if η n,ϵ n ,k n tm <γ tm η n,ϵ n ,k n sp then τ n ← 2τ n else if η n,ϵ n ,k n tm > Γ tm η n,ϵ n ,k n sp and τ n ≥ 2τ then τ n ← τ n /2 end if until (4.4) is satisfied or τ n = τ Refine the cells K ∈K n such that η n,ϵ n ,k n sp,K ≥ c ref max L∈K n η n,ϵ n ,k n sp,L in accordance with (4.6) until η n,ϵ n ,k n sp + η n,ϵ n ,k n tm + η n,ϵ n ,k n reg + η n,ϵ n ,k n lin ≤ ζ∥l n,ϵ n ,k n h ∥ L 2 (I n ;L 2 (Ω)) or (h K = h,f o rt h em a r k e dc e l l s ) Derefine the cells K ∈K n such that η n,ϵ n ,k n sp,K ≤ c deref max L∈K n η n,ϵ n ,k n sp,L u n h ← u n,ϵ n ,k n h t n ← t n-1 + τ n ϵ ← 2ϵ end while 4.3.
Efficiency of the a posteriori error estimate. In this section we investigate the global efficiency of the estimators of Corollary 3.5 under the stopping and balancing criteria of Section 4.1. Hence, the quantities at discrete time t n are those obtained after performing k n linearization iterations to meet the criterion (4.1), using a regularization parameter ϵ n and a quadrature rule such that, respectively, (4.2) and (4.3) are satisfied, and a time step ensuring the time and space error balance (4.4). As usual, in order to use the argument of equivalence of norms on finite-dimensional spaces, we need to assume here:

Assumption 4.2 (Polynomial approximations). For al l 1 ≤ n ≤ N ,t h ef u n c t i o n u n,ϵ n ,k n hτ
is affine in time on the time interval I n and piecewise polynomial of order m in space on the mesh K n-1,n ;t h ef u n c t i o n sl n,ϵ n ,k n h and t n,ϵ n ,k n h are piecewise polynomial of order m in space on K n-1,n .

For 1 ≤ n ≤ N , we introduce the standard residual-based a posteriori error estimators, cf. [START_REF] Verfürth | i o r ie r r o re s t i m a t e sf o rfi n i t ee l e m e n td i s c e t i z a t i o n so ft h eh e a t equation[END_REF]:

η n res,1 2 := τ n K∈K n-1,n h 2 K ∥ f n -∂ t u n,ϵ n ,k n hτ + ∇•l n,ϵ n ,k n h ∥ 2 L 2 (K) , (4.7a) η n res,2 2 := τ n F ∈F i,n-1,n h F ∥[[ l n,ϵ n ,k n h ]] •n F ∥ 2 L 2 (F ) . (4.7b)
Let C be a generic constant only depending on the shape regularity parameter κ K of the meshes K n-1,n ,1≤ n ≤ N ,t h es p a c ed i m e n s i o nd,a n dt h ep o l y n o m i a ld e g r e e m.I n o r d e r t o s t i l l p r o c e e d g e n e r a l l y , w i t h o u t t h e s p e c i fi c a t i o n o f a p a r t i c u l a r spatial discretization scheme, we will suppose the following: Assumption 4.3 (Approximation property). For al l 1 ≤ n ≤ N , there holds

(4.8) τ n K∈K n-1,n ∥l n,ϵ n ,k n h + t n,ϵ n ,k n h ∥ 2 L 2 (K) ≤ C η n res,1 2 + η n res,2 2 .
This property will be verified in Section 6 below for the vertex-centered finite volume spatial discretization and specific constructions of the fluxes t n,ϵ n ,k n h and l n,ϵ n ,k n h . Under these assumptions, we have the following result, showing the equivalence of the error ∥R(u n,ϵ n ,k n hτ )∥ X ′ n and the estimators of Corollary 3.5, up to data oscillation: Theorem 4.4 (Global efficiency). Let, for all 1 ≤ n ≤ N ,t h es t o p p i n gc r i t eria (4.1)-( 4.3) as well as the second inequality in the balancing criterion (4.4) be satisfied with the parameters Γ lin , Γ reg , Γ qd ,a n dΓ tm small enough. Let Assumptions 4.2 and 4.3 hold true. Then

η n,ϵ n ,k n sp +η n,ϵ n ,k n tm +η n,ϵ n ,k n qd +η n,ϵ n ,k n reg +η n,ϵ n ,k n lin ≤ C ∥R(u n,ϵ n ,k n hτ )∥ X ′ n + ∥f -f ∥ X ′ n .
The proof of this result follows the techniques of [START_REF] Verfürth | Robust a posteriori error estimates for stationary convection-diffusion equations[END_REF] and the approach of [START_REF] Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF]. It is given in Appendix A.1.

5. An a posteriori error estimate for the error in temperature and enthalpy

In the previous sections we have given a posteriori error estimators for the dual norm of the residual. In this section we prove that these same estimators also bound an error in temperature and enthalpy. We rely on a duality argument which proves to be simpler than using the dual partial differential equation as in [START_REF] Nochetto | i o r ie r r o re s t i m a t i o na n da d a p t i v i t yf o r degenerate parabolic problems[END_REF].

5.1. Bounding the error of the temperature and enthalpy by the dual norm of the residual. For brevity of notation, for t ∈ (0,T], we let

Q t := L 2 (0,t; L 2 (Ω)),X t := L 2 (0,t; H 1 0 (Ω)),X ′ t := L 2 (0,t; H -1 (Ω)
). It is convenient to stress that the result of this section applies to all functions u hτ ∈ Z such that β(u hτ ) ∈ X.W efi r s ts t a t et h ef o l l o w i n gb o u n d : Lemma 5.1 (Simple bounds for the temperature and enthalpy errors). Let u be the solution of (2.2) and let u hτ ∈ Z be such that β(u hτ ) ∈ X.T h e nt h e r eh o l d s

L β 2 ∥u -u hτ ∥ 2 X ′ + L β 2 ∥(u -u hτ )(•,T)∥ 2 H -1 (Ω) + ∥β(u) -β(u hτ )∥ 2 Q T ≤ L β 2 (2e T -1) ∥R(u hτ )∥ 2 X ′ + ∥u 0 -u hτ (•, 0)∥ 2 H -1 (Ω) and L β 2 ∥(u -u hτ )(•,T)∥ 2 H -1 (Ω) + ∥β(u) -β(u hτ )∥ 2 Q T ≤ L β 2 e T ∥R(u hτ )∥ 2 X ′ + ∥u 0 -u hτ (•, 0)∥ 2 H -1 (Ω) .
The results of Lemma 5.1 are classical; we obtain them as a byproduct in the proof of Theorem 5.2 in Section A.2 below. These results are, however, not sufficiently precise. In particular, the use of the Gronwall lemma in its proof implies the appearance of the term e T on the right-hand sides, which grows exponentially with the final time T . The purpose of the following theorem is to improve considerably this point. Indeed, note that, in Theorem 5.2, the term e T does not appear other than in the approximation of the initial condition ∥u 0u hτ (•, 0)∥ 2 H -1 (Ω) which can be made sufficiently small. Theorem 5.2 takes a more complicated form than Lemma 5.1 but the numerical results based on its use (see Section 7) prove to be excellent, which is not the case for the framework of Lemma 5.1: Theorem 5.2 (An improved bound for the temperature and enthalpy errors). Let u be the solution of (2.2) and let u hτ ∈ Z be such that β(u hτ ) ∈ X.T h e nt h e r e holds

L β 2 ∥u -u hτ ∥ 2 X ′ + L β 2 ∥(u -u hτ )(•,T)∥ 2 H -1 (Ω) + ∥β(u) -β(u hτ )∥ 2 Q T +2 T 0 ∥β(u) -β(u hτ )∥ 2 Q t + t 0 ∥β(u) -β(u hτ )∥ 2 Q s e t-s ds dt ≤ L β 2 (2e T -1)∥u 0 -u hτ (•, 0)∥ 2 H -1 (Ω) + ∥R(u hτ )∥ 2 X ′ +2 T 0 ∥R(u hτ )∥ 2 X ′ t + t 0 ∥R(u hτ )∥ 2 X ′ s e t-s ds dt .
The proof of this result is given in Section A.2.

5.2.

The a posteriori error estimate. The upper bound in Theorem 5.2 can be combined with the results of Section 3.3 to obtain an a posteriori estimate for the temperature and enthalpy errors.

Theorem 5.3 (A posteriori estimate for the temperature and enthalpy errors).

Let u be the solution of (2.2) and let u hτ and t hτ fulfill Assumptions 3.1 and 3.2, respectively. Then there holds

(5.1)

L β 2 ∥u -u hτ ∥ 2 X ′ + L β 2 ∥(u -u hτ )(•,T)∥ 2 H -1 (Ω) + ∥β(u) -β(u hτ )∥ 2 Q T +2 T 0 ∥β(u) -β(u hτ )∥ 2 Q t + t 0 ∥β(u) -β(u hτ )∥ 2 Q s e t-s ds dt ≤ L β 2 (2e T -1)η 2 IC + η 2 +2 N n=1 τ n n l=1 (η l ) 2 + N n=1 n l=1 J nl l i=1 (η i ) 2 ,
with η IC defined by (3.6), η defined by (3.8), η n , 1 ≤ n ≤ N ,d e fi n e db y

(5.2)

η n := I n K∈K n (η n R,K + η n F,K (t)) 2 dt 1 2 + ∥f -f ∥ X ′ n ,
and setting, for 1 ≤ n, l ≤ N , Proof. To prove the result, we rely on Theorem 5.2. Applying Theorem 3.3, it follows that ∥R(u hτ )∥ X ′ ≤ η,s ow ea r el e f tt oe s t i m a t et h ef o l l o w i n gr i g h t -h a n d side contributions in terms of the a posteriori error estimators:

J nl := I n I l e t-
T 1 := T 0 ∥R(u hτ )∥ 2 X ′ t dt, T 2 := T 0 t 0 ∥R(u hτ )∥ 2 X ′
s e t-s ds dt.

As in Theorem 3.3, it is readily inferred that ∥R(u hτ )∥ X ′ l ≤ η l for all 1 ≤ l ≤ N , so that

∥R(u hτ )∥ 2 X ′ t n = n l=1 ∥R(u hτ )∥ 2 X ′ l ≤ n l=1 (η l ) 2 .
Using the fact that ∥R(u hτ )∥ X ′ t is a nondecreasing function of the time t together with the above inequality yields for the first term

T 1 ≤ N n=1 I n ∥R(u hτ )∥ 2 X ′ t n dt ≤ N n=1 I n n l=1 (η l ) 2 dt = N n=1 τ n n l=1 (η l ) 2 .
Proceeding in a similar way, for the second term T 2 we obtain

T 2 ≤ N n=1 I n n l=1 I l ∥R(u hτ )∥ 2 X ′ t l e t-s dsdt ≤ N n=1 I n n l=1 I l l i=1 (η i ) 2 e t-s ds dt = N n=1 n l=1 I n I l e t-s dsdt × l i=1 (η i ) 2 = N n=1 n l=1 J nl l i=1 (η i ) 2 ,
whence the conclusion follows.

Remark 5.4 (Simplified versions of the a posteriori estimate). In the spirit of Lemma 5.1, the following simplified versions of the a posteriori estimate of Theorem 5.3 hold:

L β 2 ∥u -u hτ ∥ 2 X ′ + L β 2 ∥(u -u hτ )(•,T)∥ 2 H -1 (Ω) + ∥β(u) -β(u hτ )∥ 2 Q T ≤ L β 2 (2e T -1) η 2 + η 2 IC , L β 2 ∥(u -u hτ )(•,T)∥ 2 H -1 (Ω) + ∥β(u) -β(u hτ )∥ 2 Q T ≤ L β 2 e T η 2 + η 2 IC .
Remark 5.5 (An a posteriori error estimate distinguishing the different error components). While relying on Corollary 3.5 instead of Theorem 3.3, equivalents of Theorem 5.3 and of the bounds of Remark 5.4, distinguishing the different error components can immediately be obtained.

Application to a vertex-centered finite volume discretization

In this section, we consider the vertex-centered finite volume spatial and backward Euler temporal discretization of the Stefan problem (1.1). The regularization of Section 2.1.2 is considered and the Newton linearization is used. We show how to construct the equilibrated flux t n,ϵ,k h ,thelinearizedfluxl n,ϵ,k h , and the interpolation operator Π n of Assumption 3.4 (in generalization of Assumptions 3.1 and 3.2) and verify Assumptions 4.2 and 4.3. Thus, all the results of Sections 3-5 will apply.

T n D n D K D Figure 2
. Simplicial mesh T n and the associated vertex-centered dual mesh D n (left) and the fine simplicial mesh K D of D ∈D n (right) 6.1. Dual and tertial space meshes. The vertex-centered finite volume method is defined using a sequence of dual meshes {D n } 0≤n≤N of the space domain Ω. For a given family of matching simplicial primal meshes {T n } 0≤n≤N ,w ec o n s t r u c t {D n } 0≤n≤N as follows: for any 0 ≤ n ≤ N and with every vertex a of the mesh T n , we associate one dual volume D,c o n s t r u c t e db yc o n n e c t i n gt h eb a r y c e n t e r so ft h e simplices sharing a through edge (and face for d =3)barycen ters(seeFigure2,left) for d =2. W esplitev erysetD n into interior dual volumes D n,i and boundary dual volumes D n,b .T h es i m p l i c i a lm e s hK n appearing in Sections 2-5 is constructed by dividing each D ∈D n into a mesh K D as indicated in Figure 2, right, if d =2and similarly for d =3.

6.2. The vertex-centered finite volume scheme. Let, for 0 ≤ n ≤ N ,

V n h := ϕ h ∈ C 0 (Ω); ϕ h | K ∈ P 1 (K) ∀K ∈T n
and let (6.1)

Π n : C 0 (Ω) → V n h
be the Lagrange interpolation operator (cf. Ciarlet [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]), which to a function ϕ ∈ C 0 (Ω) associates a function ϕ h ∈ V n h by setting ϕ h (a) := ϕ(a)f o ra n yv e r t e xa of the mesh T n .

Let u 0 h ∈ V 0 h be a suitable approximation of the regularized initial enthalpy

β -1 ϵ (β(u 0 )); see Algorithm 4.1. Next, let 1 ≤ n ≤ N , u n-1 h ∈ V n-1 h
,a n dam e s hT n (and consequently D n )beg i v e n . T h ev e r t e x -c e n t e r e dfi n i t ev o l u m es c h e m ef o rt h e regularized Stefan problem (2.3) reads: find u n,ϵ h ∈ V n h such that β ϵ (u n,ϵ h )(a)=0 for all vertices a of T n on ∂Ωa n ds u c ht h a t (6.2)

1 τ n (u n,ϵ h -u n-1 h , 1) D -(∇Π n (β ϵ (u n,ϵ h ))•n D , 1) ∂D =( f n , 1) D ∀D ∈D n,i .
Then the continuous and piecewise affine-in-time function u hτ appearing in the previous sections is given by u hτ | I n := u n,ϵ hτ , (6.3)

u n,ϵ hτ (•,t)=(1-ρ(t))u n-1 h + ρ(t)u n,ϵ h ,ρ (t) := t -t n-1 τ n t ∈ I n .
Remark 6.1 (Regularization). It is also possible to consider the vertex-centered finite volume discretization without any regularization, i.e., use β in place of β ϵ in (6.2), with u 0 h ∈ V 0 h an approximation of the initial enthalpy u 0 . Remark 6.2 (Links to the discretizations of [START_REF] Baughman | Co-volume methods for degenerate parabolic problems[END_REF][START_REF] Eymard | Finite volumes and nonlinear diffusion equations (English, with English and French summaries)[END_REF][START_REF] Nochetto | i o r ie r r o re s t i m a t i o na n da d a p t i v i t yf o r degenerate parabolic problems[END_REF]). Let for simplicity the meshes T n (and consequently D n ) do not move in time and let f n be piecewise constant on License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use the scheme (6.2) coincide with that of [36, equation (4.4)], because of the links of the vertex-centered finite volumes and finite elements with mass lumping/quadrature for the source term. Similarly, in two space dimensions and when all the angles of T n are smaller than or equal to 90 • ,r e p l a c i n gt h et r i a n g l eb a r y c e n t e r sb yt h e triangle circumcenters in the construction of D n ,t h es e c o n da n dt h i r dt e r m so ft h e scheme (6.2) coincide with that in the co-volume method of [START_REF] Baughman | Co-volume methods for degenerate parabolic problems[END_REF]. More generally, whenever T n is Delaunay and the mesh D n is its Voronoïd u a l ,t h es a m el i n kh o l d s true with the cell-centered finite volume scheme of [START_REF] Eymard | Finite volumes and nonlinear diffusion equations (English, with English and French summaries)[END_REF]; cf., e.g., [START_REF] Baughman | Co-volume methods for degenerate parabolic problems[END_REF][START_REF] Eymard | Finite volumes and nonlinear diffusion equations (English, with English and French summaries)[END_REF][START_REF] Nochetto | i o r ie r r o re s t i m a t i o na n da d a p t i v i t yf o r degenerate parabolic problems[END_REF], [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF][START_REF] Eymard | Finite volume methods[END_REF], or [51, Section 3]. Hence the only slight difference between (6.2) and these schemes is in the treatment of the time evolution term which is not mass-lumped herein. Remark 6.3 (Assumption 3.1). By the definition of u hτ by (6.2)-( 6.3) and by the fact that u hτ lies in a finite-dimensional space, u hτ ∈ Z and β(u hτ ) ∈ X,s ot h a t Assumption 3.1 is satisfied. A uniform bound could also be obtained by a priori stability analysis such as those in [START_REF] Baughman | Co-volume methods for degenerate parabolic problems[END_REF][START_REF] Eymard | Finite volumes and nonlinear diffusion equations (English, with English and French summaries)[END_REF][START_REF] Nochetto | i o r ie r r o re s t i m a t i o na n da d a p t i v i t yf o r degenerate parabolic problems[END_REF], but is not necessary in our setting. 6.3. Newton linearization. Let 1 ≤ n ≤ N and the mesh T n (and D n )befi x e d . Let the vector G n-1 be given by its components associated with the dual volumes

D n .

C o n s i d e rt h ec a s ew i t h o u tr e g u l a r i z a t i o n .T h e nt h es e c o n da n dt h i r dt e r m so f

D ∈D n,i , G n-1 D :=(u n-1 h , 1) D ,andsimilarlyforthevectorF n , F n D :=( f n , 1) D .L e t u b,n,ϵ h ∈ V n
h take the values β -1 ϵ (0) (0.5 for the example of Figure 1) at the boundary vertices of T n and the value zero at the other vertices of T n .T h el a s tv e c t o rt h a t we need is H n,ϵ , H n,ϵ D :=( u b,n,ϵ h , 1) D .L e t , f o r a g i v e n d u a l v o l u m e E ∈D n,i , φ E stand for the hat basis function of the space V n h associated with E; this is a function that takes the value 1 in the vertex associated with E and the value 0 at all other vertices of T n .W ea l s od e fi n et w om a t r i c e s ,w i t ht h ec o m p o n e n t so nt h e line associated with the dual volume D ∈D n,i and on the column associated with the dual volume E ∈D n,i given by M n D,E :=( φ E , 1) D , K n D,E :=( ∇φ E •n D , 1) ∂D . All the vectors are of size R |D n,i | and the matrices of size R |D n,i |×|D n,i | ,w i t h|D n,i | the number of dual volumes in D n,i (equal to the number of interior vertices of T n ). The equation (6.2) can be written in matrix form as follows: find the vector U n,ϵ such that (6.4)

M n U n,ϵ -τ n K n β ϵ (U n,ϵ )=τ n F n + G n-1 -H n,ϵ , where (β ϵ (U n,ϵ )) D := β ϵ (U n,ϵ D ). We have u n,ϵ h = E∈D n,i U n,ϵ E φ E + u b,n,ϵ
h . The algebraic system (6.4) is nonlinear. Its solution is approximated using the Newton linearization. Let U n,ϵ,0 be fixed; typically, U n,ϵ,0 := U n-1 .T h e n , f o r k ≥ 1, we approximate (6.5)

β ϵ (U n,ϵ,k ) ≈ β ϵ (U n,ϵ,k-1 )+β ′ ϵ (U n,ϵ,k-1 ) U n,ϵ,k -U n,ϵ,k-1 .
Since the regularized enthalpy-temperature function β ϵ is continuously differentiable, the Newton linearization (6.5) is well defined. At every Newton iteration k, we are thus lead to solve the following system of linear algebraic equations: find the vector U n,ϵ,k such that

M n -τ n K n β ′ ϵ (U n,ϵ,k-1 ) U n,ϵ,k = τ n F n + G n-1 -H n,ϵ (6.6) -τ n K n β ′ ϵ (U n,ϵ,k-1 )U n,ϵ,k-1 -β ϵ (U n,ϵ,k-1
) . At each linearization step k,w es e t (6.7)

u n,ϵ,k h := E∈D n,i U n,ϵ,k E φ E + u b,
:= ∇ E∈D n,i β ϵ (U n,ϵ,k-1 E )+β ′ ϵ (U n,ϵ,k-1 E ) U n,ϵ,k E -U n,ϵ,k-1 E φ E .
We p erform the Newton iterations until we meet the convergence criterion discussed in Section 4.1.

6.4. Flux reconstruction. Let a time step 1 ≤ n ≤ N ,aregularizationparameter ϵ>0, and a Newton linearization step k be fixed. We now show how to construct the flux t n,ϵ,k h of Assumption 3.4. For this purpose, we will solve a local Neumann problem by mixed finite elements on every dual volume, following [START_REF] Ern | Flux reconstruction and a posteriori error estimation for discontinuous Galerkin methods on general nonmatching grids (English, with English and French summaries)[END_REF][START_REF] Luce | a lap o s t e r i o r ie r r o re s t i m a t o rb a s e do ne q u i l ibrated fluxes[END_REF][START_REF] Vohralík | Guaranteed and fully robust a posteriorie r r o re s t i m a t e sf o rc o n f o r m i n gd i scretizations of diffusion problems with discontinuous coefficients[END_REF]. For a given D ∈D n , we introduce the spaces

RTN (K D ) := {v h ∈ H(div; D); v h | K ∈ [P 0 (K)] d + x P 0 (K) ∀K ∈K D }, RTN N (K D ) := {v h ∈ RTN (K D ); v h •n F = -l n,ϵ,k h •n F ∀F ∈ ∂K i D }, RTN N,0 (K D ) := {v h ∈ RTN (K D ); v h •n F =0 ∀F ∈ ∂K i D }
, where ∂K i D stands for all the faces of the submesh K D which are on the boundary of the dual volume D but not on the boundary of Ω. We will also need the space P * 0 (K D )w h i c hc o n s i s t so fp i e c e w i s ec o n s t a n t sf u n c t i o n so nK D ;w h e nD ∈D n,i ,w e additionally impose a zero mean value over D. The local problem consists in finding

t n,ϵ,k h ∈ RTN N (K D )a n dq h ∈ P * 0 (K D ),
(t n,ϵ,k h + l n,ϵ,k h , v h ) D -(q h , ∇•v h ) D =0 ∀v h ∈ RTN N,0 (K D ), (6.9a) (∇•t n,ϵ,k h ,φ h ) D -( f n -∂ t u n,ϵ,k hτ ,φ h ) D =0 ∀φ h ∈ P * 0 (K D ). (6.9b)
Note that the problem (6.9) is well-posed and one can take all φ h ∈ P 0 (K D )a st h e test functions in (6.9b). Indeed, it follows from (6.6) and (6.8) (compare to (6.2)) that (6.10)

1 τ n (u n,ϵ,k h -u n-1 h , 1) D -(l n,ϵ,k h •n D , 1) ∂D =( f n , 1) D ∀D ∈D n,i .
From (6.10), we see that the Neumann boundary condition encoded in RTN N (K D ) is in equilibrium with the boundary datum f n -∂ t u n,ϵ,k hτ of (6.9). We have the following key result: Lemma 6.4 (Assumptions 3.4, 4.2, and 4.3). Let 1 ≤ n ≤ N , ϵ>0,a n dk ≥ 1 be fixed. Let u n,ϵ,k h be given by (6.6)-(6.7), l n,ϵ,k h by (6.8), t n,ϵ,k h by (6.9),a n dΠ n by (6.

1).T h e nA s s u m p t i o n

Numerical experiments

We illustrate in this section our theoretical results on a series of numerical exp eriments for the vertex-centered finite volume discretization approach of Section 6. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 7.1. Setting. We consider the two-dimensional test case prop osed by No chetto et al. [START_REF] Nochetto | An adaptive finite element method for two-phase Stefan problems in two space dimensions. I. Stability and error estimates[END_REF][START_REF] Nochetto | An adaptive finite element method for twophase Stefan problems in two space dimensions. II. Implementation and numerical experiments[END_REF] on the space-time domain Ω × (0,T)w i t hΩ=( 0 , 5) 2 and T = π/1.25. The function β(•) is given by β(u)=u ⊖ +(u -1) ⊕ .T h ee x a c tt e m p e r a t u r eh a s the following expression:

(7.1) (β(u))(x, y, t)= 0.75(r 2 -1), if r<1, 1.5 -ρ ′ (t) y-ρ(t) r (r -1), if r ≥ 1,
where r 2 := x 2 +(yρ(t)) 2 and ρ(t) :=0.5 + sin(1.25t). The exact interface I(t) is ac i r c l ew i t hc e n t e r( 0 ,ρ(t)) and radius 1. The motion of the interface is governed by the Stefan law which prescribes that the normal velocity v satisfies

(∇β(u) + -∇β(u) -)•n = v on I(t),
where ∇β(u) + and ∇β(u) -denote the values of the temperature gradient on each side of the interface, while n is the unit normal to the interface with suitable orientation. The enthalpy u on Ω \ I(t)c a nb eo b t a i n e df r o mt h ee x p r e s s i o n( 7 . 1 ) of β(u). The homogeneous Neumann condition ∇β(u)•n =0i se n f o r c e da tx =0 , whereas Dirichlet boundary conditions on the temperature are prescribed at y =0, y =5 ,a n dx =5u s i n g( 7 . 1 ) . T h ei n i t i a le n t h a l p yu 0 and the source term f are likewise imposed using (7.1). The vertex-centered finite volume discretization of Section 6 is considered. No adaptation of the quadrature rule is performed; this is areasonablesimplificationsincetheenthalpy-temperaturefunctionβ(•) is piecewise affine.

7.2.

Computing approximately the negative norms. In practice we cannot compute the negative norms as the initial data indicator η IC ,t h ed a t ao s c i l l a t i o n ∥f -f ∥ X ′ ,andthedualnormoftheresidual∥R(u hτ )∥ X ′ , even if the exact solution u is known. For numerical experiments below, the dual norms are approximated by solving auxiliary problems. More specifically, for a function v ∈ X ′ to compute the negative norm ∥v∥ X ′ we consider for a.e. t ∈ (0,T)t h ep r o b l e m : fi n dψ(•,t) ∈ H 1 0 (Ω) such that

(7.2) (∇ψ(•,t), ∇ϕ)=⟨v(•,t),ϕ⟩, ∀ϕ ∈ H 1 0 (Ω).
Then

∥v∥ 2 X ′ = T 0 sup ϕ∈H 1 0 (Ω), ∥∇ϕ∥ L 2 (Ω) =1 ⟨v(•,t),ϕ⟩ 2 dt = T 0 sup ϕ∈H 1 0 (Ω), ∥∇ϕ∥ L 2 (Ω) =1 (∇ψ(•,t), ∇ϕ) 2 dt = T 0 ∥∇ψ∥ 2 L 2 (Ω) (t)dt = ∥∇ψ∥ 2 L 2 (0,T ;L 2 (Ω)) .
We obtain an approximation of the function ψ by solving the problem (7.2) numerically by the vertex-centered finite volume scheme on a refined spatial mesh and on discrete times which refine the given temporal mesh. We suppose that the ensuing discretization error is small and can be ignored. The computation of η IC is easier as it only involves the initial time t =0. As expected, the linearization error steadily decreases, while the other components stagnate starting from the second iteration. The stopping criterion (4.1) with Γ lin = 10 -2 allows us to profit from this behavior by stopping the Newton algorithm after the second iteration, while a classical criterion based on a fixed threshold,

(7.3) η n,ϵ,k lin ≤ ζ lin ,ζ lin =10 -7 ,
would require 10 iterations to converge. The overall gain for an entire simulation in terms of linearization iterations can be appreciated considering the results in Figure 4, left. We use the adaptive Algorithm 4.1 with different choices for the linearization stopping criterion: the classical criterion (7.3) then the stopping criterion (4.1) with Γ lin =0 .01 and Γ lin =0 .1. The other parameters used in the Algorithm are:

Γ reg =0 .1, ζ =1 ,ζ IC =1 ,h K 0 =0 .1, τ 0 =0 .1, ϵ 0 =0 .25, h =10 -2 , τ =10 -2 , c ref =0.7, c deref =0.2, γ tm =0.7
, and Γ tm =1.3. For the sake of completeness we also add a comparison with the local version of the stopping criterion (4.1), namely

(7.4) η n,ϵ,k n lin,K ≤ Γ lin,loc η n,ϵ,k n sp,K + η n,ϵ,k n tm,K + η n,ϵ,k n qd,K + η n,ϵ,k n reg,K ∀K ∈K n .
Even with this more stringent criterion, Figure 4, right, shows that a considerable gain in terms of number of linearization iterations can be achieved, whereas the precision on each time step (expressed by our error estimator η n (5.2)) is basically unchanged. In Figure 6, bottom left, we started by an initial mesh K 0 with h K 0 =0 .4a n da time step τ 0 =0 .1. Then the time step is adapted in order to satisfy (4.4), with γ tm =0.7, Γ tm =1.3. As a result, the spatial (3.14a) and temporal (3.14b) error estimators stay equilibrated during the whole simulation. Figure 6, top, on the other hand, shows two possible disequilibrated patterns corresponding to space and time over-refinement. In the top left we started by an initial mesh K 0 with h K 0 =0 .2 and a time step τ 0 =0.2, we fixed also γ tm =2andΓ tm =3,whileinthetoprigh t we started by an initial mesh K 0 with h K 0 =0 .5 and time step τ 0 =0 .05 and we fixed γ tm = 1 3 and Γ tm = 1 2 .F i n a l l y ,F i g u r e6 ,b o t t o mr i g h ts h o w st h ee ff e c to ft h i s violating of the balancing criterion (4.4) on the total error. These results make it apparent that the performance of an adaptive code may be considerably reduced when time and space errors are not balanced, and advocate the use of (4.4).

Next, we compare in Figure 7 the actual and predicted error distributions using the adaptive Algorithm 4.1 with

Γ lin =Γ reg =0 .1, ζ =1 ,ζ IC =1 ,h K 0 =0 .25,
τ 0 =0 .05, ϵ 0 =0 .25, h =1 0 -2 , τ =1 0 -2 , c ref =0 .7, c deref =0 .2, γ tm =0 .7,
and Γ tm =1 .3. We present the results at time t =0 .1. We see that the actual and predicted error distributions match very nicely. The corresponding exact and discrete enthalpies are depicted in Figure 8. 7.5. Overall performance. In this section we assess the overall performance of the adaptive algorithm of Section 4.2 in terms of precision vs. the number of unknowns.

In Figure 9, left, we depict the error and estimates as a function of the total number of space-time unknowns in the fully adaptive case and in the uniform case. In the adaptive case, we use Algorithm 4.1 with the parameters detailed in Section 7.4. In the uniform case, the temporal and spatial meshes as well as the regularization parameter are fixed during the simulation, and linearization is stopped when (7.3) is satisfied. The error is measured in the dual norm (3. 9. In both cases the adaptive strategy yields much better results than the uniform one, as expected. The right part of Figure 9 displays the corresponding effectivity indices, given by the ratio of the estimates over the error. These are remarkably close to the optimal value of one for the dual norm (3.3), even for the , where e denotes the chosen error measure while N st and M st are the total number of space-time unknowns corresponding to two subsequent levels of refinement.

The results for the uniform and adaptive cases are collected in Tables 1a and1b, respectively. We evaluate the dual norm of the residual (3.2), the L 2 (0,T; L 2 (Ω)) error in the temperature, and the estimator η of (3.8). We observe roughly twice faster convergence in the adaptive case in comparison with the uniform one.

Appendix A. Proofs

In this appendix, we collect the more involved proofs of some theorems of the paper.

A.1. Proof of Theorem 4.4. In this section, we will use the notation a b for the inequality a ≤ Cb with a generic constant C only depending on the shape regularity parameter κ K of the meshes K n-1,n and on the maximal level of coarsening between K n-1 and K n ,1≤ n ≤ N ,thespacedimensiond,andthepolynomialdegreem. Fix 1 ≤ n ≤ N .W es t a r tb yo b s e r v i n gt h a t ,o w i n gt ot h es t o p p i n gc r i t e r i a( 4 . 1 ) -( 4 . 3 ) and to the second inequality in the balancing criterion (4.4),

(A.1) η n,ϵ n ,k n sp + η n,ϵ n ,k n tm + η n,ϵ n ,k n qd + η n,ϵ n ,k n reg + η n,ϵ n ,k n lin η n,ϵ n ,k n sp .
Recall that we have supposed in Section 2.2.2 that the mesh K n ,1≤ n ≤ N , is obtained from K n-1 by limited refinement/coarsening and that the common refinements K n-1,n are uniformly shape regular. Thus, for K ∈K n , using the triangle inequality, Assumption 4.2, and the inverse inequality, cf. [44, Proposition 6.3.2], the first term of (3.13a) can be bounded by

η n,ϵ n ,k n R,K = C P,K h K ∥ f n -∂ t u n,ϵ n ,k n hτ -∇•t n,ϵ n ,k n h ∥ L 2 (K) ≤ C P,K h K ∥ f n -∂ t u n,ϵ n ,k n hτ + ∇•l n,ϵ n ,k n h ∥ L 2 (K) + C P,K h K ∥∇•(l n,ϵ n ,k n h + t n,ϵ n ,k n h )∥ L 2 (K) ⎧ ⎨ ⎩ K ′ ∈K n-1,n ,K ′ ⊂K h 2 K ′ ∥ f n -∂ t u n,ϵ n ,k n hτ + ∇•l n,ϵ n ,k n h ∥ 2 L 2 (K ′ ) ⎫ ⎬ ⎭ 1 2 + ⎧ ⎨ ⎩ K ′ ∈K n-1,n ,K ′ ⊂K ∥l n,ϵ n ,k n h + t n,ϵ n ,k n h ∥ 2 L 2 (K ′ ) ⎫ ⎬ ⎭ 1 2
.

Consequently, employing Assumption 4.3,

(A.2) η n,ϵ n ,k n sp η n res,1 + η n res,2 .
Proving the efficiency of the estimators introduced in Section 3.4 thus amounts to proving the efficiency of the residual estimators η n res,1 and η n res,2 . Henceforth, to simplify, we will use the shorthand notation

u n hτ = u n,ϵ n ,k n hτ , l n h = l n,ϵ n ,k n h and denote η n LRQT 2 := I n K∈K n-1,n ∥∇β(u n hτ (•,t)) -l n h ∥ 2 L 2 (K) dt. (A.3)
We have:

∥R(u n hτ )∥ X ′ n + η n LRQT + ∥f -f ∥ X ′ n .
Proof. For all K ∈K n-1,n , we let v

K :=( f n -∂ t u n hτ +∇•l n h )| K .B yA s s u m p t i o n4 . 2 , v K is polynomial in K.
W ed e n o t eb yψ K the usual bubble function on K,i . e . ,t h e product of the (d+1) hat basis functions (barycentric coordinates) ψ a associated with the vertices a of the element K,s e tλ K := h

2 K ψ K v K for all K ∈K n-1,n ,a n d let λ := K∈K n-1,n λ K .C l e a r l y ,λ ∈ H 1 0 (Ω) and λ| K ∈ H 1 0 (K)f o ra l lK ∈K n-1,n .
Using the equivalence of norms on finite-dimensional spaces, integrating by parts in space, the weak form (2.2c), and (3.9) together with the Cauchy-Schwarz inequality; we infer that (cf. [START_REF] Verfürth | Robust a posteriori error estimates for stationary convection-diffusion equations[END_REF])

(A.5) η n res,1 2
I n K∈K n-1,n h 2 K (v K ,ψ K v K ) K ds = I n ⟨∂ t (u -u n hτ ),λ⟩ +(∇β(u) -∇β(u n hτ ), ∇λ) +(∇β(u n hτ ) -l n h , ∇λ)+( f n -f, λ) ds ≤ ∥R(u n hτ )∥ X ′ n + η n LRQT + ∥f -f ∥ X ′ n ∥λ∥ X n .
By the shape regularity of the mesh K n-1,n and the inverse inequality (cf. [44, Proposition 6.3.2]) we have, for any K ∈K n-1,n ,

∥∇λ∥ L 2 (K) = h 2 K ∥∇(ψ K v K )∥ L 2 (K) h K ∥ψ K v K ∥ L 2 (K) ≤ h K ∥v K ∥ L 2 (K) .
An immediate consequence is that ∥λ∥ X n η n res,1 and (A.4) follows.

Lemma A.2 (Estimate of η n res,2 ). Under the assumptions of Theorem 4.4,t h e r e holds

(A.6) η n res,2 ∥R(u n hτ )∥ X ′ n + η n LRQT + ∥f -f ∥ X ′ n .
Proof. Let F ∈F i,n-1,n .W ed e n o t eb yK F the simplices K ∈K n-1,n that share the face F .L e tv F :=[ [ l n h ]] •n F and keep the same notation for the constant extension of v F into K F along the vectors face barycenter-opposite vertex. Owing to Assumption 4.2, v F is a polynomial on K F .L e tψ F be the usual face bubble function supported on K F ,i.e.,theproductofthed hat basis functions (barycentric coordinates) ψ a associated with the vertices a of the face F .F o ra l lF ∈F i,n-1,n , set λ F := h F ψ F v F and let λ := F ∈F i,n-1,n λ F .N o t e t h a t λ ∈ H 1 0 (Ω) and λ| K F ∈ H 1 0 (K F )f o ra l lF ∈F i,n-1,n . Using the equivalence of norms in finitedimensional spaces, integrating by parts in space, using the weak form (2.2c), and (3.9) together with the Cauchy-Schwarz inequality, it is inferred that (cf. [START_REF] Verfürth | Robust a posteriori error estimates for stationary convection-diffusion equations[END_REF]) (A.7)

η n res,2 2
I n F ∈F i,n-1,n h F (v F ,ψ F v F ) F ds = I n F ∈F i,n-
+ η n LRQT + ∥f -f ∥ X ′ n ∥λ∥ X n + η n res,1 τ n K∈K n-1,n h -2 K ∥λ∥ 2 L 2 (K) 1 2
.

Using the fact that, for all F ∈F i,n-1,n and 3) and Γ tm in the balancing criterion (4.4). Thus, choosing these parameters small enough, the term η n LRQT can be made small enough to be discarded from the right-hand side of (A.8); cf. [START_REF] Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF]Theorem 4.4] and the assertion of Theorem 4.4 follows from (A.1). A.2. Proof of Theorem 5.2. We start by proving the following intermediate result. (A.9) ≤∥u 0u hτ (•, 0)∥ 2 H -1 (Ω) + ∥R(u hτ )∥ 2

K ∈K n-1,n F , ∥ψ F v F ∥ L 2 (K) h 1 2 F ∥v F ∥ L 2 (F ) , it is inferred that τ n K∈K n-1,n h -2 K ∥λ∥ 2
X ′ t + ∥u -u hτ ∥ 2 X ′ t .
Proof. For a.e. t ∈ (0,T), we denote by W (•,t) ∈ H 1 0 (Ω) the solution to (A.10) (∇W (•,t), ∇ψ)=((uu hτ )(•,t),ψ) ∀ψ ∈ H 1 0 (Ω). The existence and uniqueness of W (•,t)f o l l o wf r o mt h eR i e s zr e p r e s e n t a t i o nt h e orem. Moreover, since u, u hτ ∈ Z, there holds W ∈ X. Using (A.10), it is inferred that This duality technique is rather standard; see [START_REF] Cancès | An a posteriori error estimate for vertex-centered finite volume discretizations of immiscible incompressible two-phase flow[END_REF] and the references therein. Its origins can be traced back at least to the elliptic projection of Wheeler [START_REF] Wheeler | Ap r i o r iL 2 error estimates for Galerkin approximations to parabolic partial differential equations[END_REF]. In some aspects, it is close to the elliptic reconstruction of Makridakis and Nochetto [START_REF] Makridakis | Elliptic reconstruction and a posteriori error estimates for parabolic problems,S I A MJ[END_REF]; however, in [START_REF] Makridakis | Elliptic reconstruction and a posteriori error estimates for parabolic problems,S I A MJ[END_REF] it is used to restore optimal order of the a posteriori estimate in L ∞ (0,T; L 2 (Ω)), whereas here we employ it to obtain a bound on an energy-like norm.

Taking ϕ = W 1 (0,t) with 1 (0,t) the characteristic function of the interval (0,t) in definition (3.1) and using (A.11) and the Young inequality, it is inferred that (A.12) Recalling (A.10), and since uu hτ ∈ H 1 (0,T; H -1 (Ω)), there holds ∂ t W ∈ X and, for a.e. s ∈ (0,T), ∂ t W (•,s) satisfies in a weak sense

⟨R(u hτ ),W⟩ X ′ t ,X t ≤∥R(u hτ )∥ X ′ t ∥u -u hτ ∥ X ′ t ≤ 1 2 ∥R(u hτ )∥ 2 X ′ t + 1 2 ∥u -u hτ ∥ 2
-∇•(∇∂ t W (•,s)) = ∂ t (u -u hτ )(•,s)
in Ω,

∂ t W (•,s)=0 on∂Ω.
Thus, it follows from the definition (A.10) of W and from the norm characterization (A.11) that (A.14)

R 1 = t 0 (∂ t ∇W, ∇W )(s)ds = 1 2 ∥∇W (•,t)∥ 2 L 2 (Ω) -∥∇W (•, 0)∥ 2 L 2 (Ω) = 1 2 ∥(u -u hτ )(•,t)∥ 2 H -1 (Ω) -∥u 0 -u hτ (•, 0)∥ 2 H -1 (Ω) .
Invoking again the definition (A.10) and using the fact that β is nondecreasing and L β -Lipschitz continuous, there holds The conclusion follows using the inequalities (A.12), (A.14), and (A.15) in equation (A.13).

(A.
Corollary A.4 (Application of the Gronwall lemma). Under the assumptions of Lemma A. The assertion follows by integrating over the interval (0,T).

We are now ready to prove Theorem 5.2:

Proof of Theorem 5.2. Using (A.9) with t = T and adding ∥uu hτ ∥ 2 X ′ to both sides we infer

L := 2 L β ∥β(u) -β(u hτ )∥ 2 Q T + ∥u -u hτ ∥ 2 X ′ + ∥(u -u hτ )(•,T)∥ 2 H -1 (Ω)
≤∥u 0u hτ (•, 0)∥ 2 H -1 (Ω) + ∥R(u hτ )∥ 2 X ′ +2∥uu hτ ∥ 2 X ′ . Using Corollary A.4 to estimate the last term in the right-hand side we obtain L ≤(2e T -1)∥u 0u hτ (•, 0)∥ 2 H -1 (Ω) + ∥R(u hτ )∥ 2

X ′ +2 T 0 ∥R(u hτ )∥ 2 X ′ t + t 0 ∥R(u hτ )∥ 2 X ′ s e t-s ds dt - 4 L β T 0 ∥β(u) -β(u hτ )∥ 2 Q t + t 0 ∥β(u) -β(u hτ )∥ 2
Q s e t-s ds dt.

The conclusion follows multiplying both sides by L β /2a n dr e a r r a n g i n gt h et e r m s .

  the mixed finite element approximations of local Neumann problems on D ∈D n,i and local Neumann/Dirichlet problems on D ∈D n,b :

  s3.4, 4.2,a n d 4.3 hold true. Proof. The equilibrium property (3.11) follows immediately from (6.9b), so that Assumption 3.4 is easily satisfied. Whereas Assumption 4.2 is trivial, Assumption 4.3 is obtained by proceeding exactly as in [51, proof of Theorem 5.5] or [15, proof of Lemma 5.3].

Figure 3 . 3 .

 33 Figure 3. Evolution of the spatial, temporal, regularization, and linearization error estimators (3.14) as a function of Newton iterations for a fixed mesh, time step, and regularization parameter

Figure 4 .Figure 5 .Figure 5 4 8

 45548 Figure 4. Error estimator η n (5.2) as a function of the cumulated Newton iterations at each time step (time steps are identified by markers). Global stopping criterion (4.1) or (7.3) (left), local stopping criterion (7.4) (right)

Figure 6 .

 6 Figure 6. Effect of the time step adaptation strategy on the global error estimator (3.8). Violations of the balancing criterion (4.4) by space over-refinement (top left) and time overrefinement (top right). Time step refinement honoring (4.4) (bottom left). Overall comparison (bottom right)

Figure 7 .

 7 Figure 7. Actual (left)a n de s t i m a t e d( right)e r r o rd i s t r i b u t i o n for Γ lin =Γ reg =0.1, adaptive Algorithm 4.1, entire domain (top), interface zoom (bottom)

  3) and estimated by Theorem 3.3 in the top part of Figure 9, whereas the energylike norm (5.1) and the estimate of Theorem 5.3 are used in the bottom part of

Figure 8 .

 8 Figure 8. Exact (left)a n da p p r o x i m a t e( right)e n t h a l p yc o r r esponding to the results of Figure 7

Table 1 .

 1 Comparison of the experimental orders of convergence (e.o.c.) in the uniform and fully adaptive cases. The total number of space-time unknowns is denoted by N st .T h ea c t u a le r r o r ∥R(u hτ )∥ X ′ and the estimated error η are defined by (3.2) and (3.8) respectively. (a) Uniform case N st ∥β(u)β(u hτ )∥ Q T e.o.c

L 2

 2 (K) (η n res,2 ) 2 ,w h e n c e by the inverse inequality, ∥λ∥ X n η n res,2 . Using this fact in (A.7) in conjunction with (A.4), (A.6) follows. Proof of Theorem 4.4. It follows from Lemmas A.1 and A.2 and from (A.2) that(A.8) η n,ϵ n ,k n sp ∥R(u n,ϵ n ,k n hτ )∥ X ′ n + η n LRQT + ∥f -f ∥ X ′ n .In order to bound the term η n LRQT ,w ep r o c e e da sf o l l o w s . T h et r i a n g l ei n e q u a l i t y and the definitions (3.13) and(3.14) giveη n LRQT ≤ η n,ϵ n ,k n lin + η n,ϵ n ,k n reg + η n,ϵ n ,k n qd + η n,ϵ n ,k n tm .Thus, proceeding as for the bound (A.1),η n LRQT ≤ Cη n,ϵ n ,k n sp ,where the constant C only depends on the parameters Γ lin ,Γ reg ,a n dΓ qd in the stopping criteria (4.1)-(4.

Lemma A. 3 ( 2 H

 32 Duality bound). Let u be the solution to (2.2) and let u hτ ∈ Z be such that β(u hτ ) ∈ X.T h e n ,t h e r eh o l d s ,f o ra . e .t ∈ (0,T), 2L β ∥β(u)β(u hτ )∥ 2 Q t + ∥(uu hτ )(•,t)∥ -1 (Ω)

(A. 11 )

 11 ∥∇W (•,t)∥ L 2 (Ω) =s u p ψ∈H 1 0 (Ω), ∥∇ψ∥ L 2 (Ω) =1(∇W (•,t), ∇ψ)=s u p ψ∈H 1 0 (Ω), ∥∇ψ∥ L 2 (Ω) =1 ((uu hτ )(•,t),ψ) = ∥(uu hτ )(•,t)∥ H -1 (Ω) .

0 (

 0 hτ ),W⟩ X ′ t ,X t = t 0 ⟨∂ t (uu hτ ),W⟩(s)ds (A.13) + t ∇β(u) -∇β(u hτ ), ∇W )(s)ds = : R 1 + R 2 .

15 ) 1 L β t 0 (

 1510 u hτ ,β(u)β(u hτ ))(s)ds ≥ β(u)β(u hτ ),β(u)β(u hτ ))(s)ds = 1 L β ∥β(u)β(u hτ )∥ 2 Q t .

  Licensed to Ecole Natl des Ponts et Chaussees. Prepared on Thu Jun 19 04:33:39 EDT 2014 for download from IP 195.221.193.15/195.221.193.114.

  . ∥R(u hτ )∥ X ′ e.o.c. hτ )∥ Q T e.o.c. ∥R(u hτ )∥ X ′ e.o.c.

						η	e.o.c.
	7020	7.13e-02	-	3.75e-01	-	1.22e-00	-
	66906	6.02e-02	0.224	3.30e-01	0.172 8.65e-01 0.455
	915840	5.07e-02	0.197	2.48e-01	0.364 6.50e-01 0.392
	1.12963e+07	2.19e-02	0.221	1.60e-01	0.115 2.40e-01 0.261
			(b) Adaptive case
	N st	∥β(u) -β(u η	e.o.c.
	9360	6.55e-02	-	3.51e-01	-	1.51e-00	-
	35370	5.28e-02	0.486	3.07e-01	0.303 1.08e-00 0.751
	224082	4.06e-02	0.427	2.19e-01	0.546 6.32e-01 0.868
	1.53329e+06	1.10e-02	0.392	1.18e-01	0.186 2.23e-01 0.312
	Figure

  Licensed to Ecole Natl des Ponts et Chaussees. Prepared on Thu Jun 19 04:33:39 EDT 2014 for download from IP 195.221.193.15/195.221.193.114.

			err. unif.	est. unif.			effectivity unif.
	Error/error estimate	10 0	err. ad.	est. ad.	Effectivity index	2 3 4	effectivity ad.
		10 -1
		10 4	10 5	10 6		10 7		10 4	10 5	10 6	10 7
		Total numb er of space-time unknowns		Total numb er of space-time unknowns
			err. unif.	est. unif.			effectivity unif.
			err. ad.	est. ad.		12	effectivity ad.
	Error/error estimate	10 -1 10 0					Effectivity index	8 10
		10 4	10 5	10 6		10 7		10 4	10 5	10 6	10 7
		Total numb er of space-time unknowns		Total numb er of space-time unknowns
		Figure 9. Comparison between adaptive and uniform refinement.
		Dual norm (3.3) (top), energy-like norm (5.1) (bottom). Error and
		estimators (left), effectivity indices (right)
	present time-dependent, degenerate problem with a moving free boundary. We
	regard the effectivity indices corresponding to Theorem 5.3 as likewise excellent;
	there are, in particular, several orders of magnitude smaller than the effectivity
	indices corresponding to the setting of Remark 5.4 that we have also assessed (not
	presented).
	A quantitative evaluation of the performance in terms of precision vs. the number
	of unknowns can be obtained by computing the experimental order of convergence
	(e.o.c.), defined as follows:
			e.o.c :=	log(e N st ) -log(e M st ) -1 3 (log N

License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use stlog M st )

  Licensed to Ecole Natl des Ponts et Chaussees. Prepared on Thu Jun 19 04:33:39 EDT 2014 for download from IP 195.221.193.15/195.221.193.114. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use Lemma A.1 (Estimate of η n res,1 ). Under the assumptions of Theorem 4.4,t h e r e holds

	(A.4)	η n res,1

  1,n K∈K F {(∇•l n h ,λ F ) K +(l n h , ∇λ F ) K } ds =

		{(∇•l n h ,λ)+(l n h , ∇λ)} ds
	I n
	=	⟨∂ t (u n hτ -u),λ⟩ +(∇β(u n hτ ) -∇β(u), ∇λ)
	I n
	+( f n -∂ t u n hτ + ∇•l n h ,λ)
	+(l n h -∇β(u n hτ ), ∇λ)+(f -f n ,λ) ds
	∥R(u n hτ )∥ X ′ n

  [START_REF] Andreianov | Discrete duality finite volume schemes for doubly nonlinear degenerate hyperbolic-parabolic equations[END_REF], there holds∥uu hτ ∥ 2 X ′ ≤ (e T -1)∥u 0u hτ (•, 0)∥ 2 )β(u hτ )∥ 2 Q s e t-s ds dt. ∥(u-u hτ )(•,t)∥ 2 H -1 (Ω) and α(t) := ∥u 0 -u hτ (•, 0)∥ 2 H -1 (Ω) +∥R(u hτ )∥ 2 L β ∥β(u)β(u hτ )∥ 2 Q t , it is inferred, for a.e. t ∈ (0,T), ∥(uu hτ )(•,t)∥ 2 H -1 (Ω) ≤ e t ∥u 0u hτ (•, 0)∥ 2 H -1 (Ω) + ∥R(u hτ )∥ 2

						H -1 (Ω)
			T		t
		+	0	∥R(u hτ )∥ 2 X ′ t +	0	∥R(u hτ )∥ 2 X ′ s e t-s ds dt
	-∥β(uProof. Using (A.9) followed by the Gronwall lemma T t 2 L β 0 ∥β(u) -β(u hτ )∥ 2 Q t + 0
		ξ(t) ≤ α(t)+
	2					X ′ t	-
						X ′ t
	t
	+	∥R(u hτ )∥ 2 X ′
	0

t 0 ξ(s)ds =⇒ ξ(t) ≤ α(t)+ t 0 α(s)e t-s ds, with ξ(t) := s e t-s ds -2 L β ∥β(u)β(u hτ )∥ 2 Q t + t 0 ∥β(u)β(u hτ )∥ 2

Q s e t-s ds .

This work was supported by the ERT project "Enhanced oil recovery and geological sequestration of CO ra d v a n c e d t e c h n i q u e s " (LJLL/IFPEN).

The third author is the corresponding author.