Adaptive Mesh Refinement for a Finite Volume Method for Flow and Transport of Radionuclides in Heterogeneous Porous Media

— Adaptation de maillage pour un sche´ma volumes ﬁnis pour la simulation d’e´coulement et de transport de radionucle´ides en milieux poreux he´te´roge`nes — Cet article traite de l’adaptation dynamique de maillages pour la simulation nume´rique d’e´coulements incompressibles et miscibles en milieux poreux. Le proble`me est mode´lise´ par un syste`me couple´ entre une e´quation elliptique (pression-vitesse de Darcy) et une e´quation de diffusion-convection (concentration). Le syste`me est discre´tise´ par une me´thode volumes ﬁnis centre´s sur les sommets. On utilise un sche´ma de Godunov pour approcher le terme de convection et une approximation e´le´ment ﬁni P1 pour le terme de diffusion. Nous de´veloppons des estimateurs d’erreur a posteriori de type re´siduel. Nous introduisons deux sortes d’indicateurs. Abstract — Adaptive Mesh Reﬁnement for a Finite Volume Method for Flow and Transport of Radionuclides in Heterogeneous Porous Media — In this paper, we consider adaptive numerical simulation ofmiscible displacement problems in porous media, whichare modeled bysinglephase ﬂow equa-tions.Avertex-centredﬁnitevolumemethodisemployedtodiscretizethecoupledsystem: theDarcyﬂow equationandthediffusion-convectionconcentrationequation.Theconvectiontermisapproximatedwith a Godunov scheme over the dual ﬁnite volume mesh, whereas the diffusion-dispersion term is discretized by piecewise linear conforming ﬁnite elements. We introduce two kinds of indicators, both of them of residual type. The ﬁrst one is related to time discretization and is local with respect to the time discret-ization:thus, ateachtime, itprovidesanappropriateinformationforthechoiceofthenexttimestep.The second is related to space discretization and is local with respect to both evaluates where additional reﬁnement is needed and grid generation procedures dynamically create or remove ﬁne-grid patches as resolution requirements change. The method was implemented in the software MELODIE, developed by the French Institute for Radiological Protection and Nuclear Safety (IRSN, Institut de Radioprotection et de Suˆrete´ Nucle´aire ). The algorithm is then used to simulate the evolution of radionuclide migration from the waste packages through a heterogeneous disposal, dem-onstrating its capability to capture complex behavior of the resulting ﬂow.


Uncertainties
Analyse de sensibilité et optimisation sous incertitudes de procédés EOR de type surfactant-polymère   Abstract -Adaptive Mesh Refinement for a Finite Volume Method for Flow and Transport of Radionuclides in Heterogeneous Porous Media -In this paper, we consider adaptive numerical simulation of miscible displacement problems in porous media, which are modeled by single phase flow equations. A vertex-centred finite volume method is employed to discretize the coupled system: the Darcy flow equation and the diffusion-convection concentration equation. The convection term is approximated with a Godunov scheme over the dual finite volume mesh, whereas the diffusion-dispersion term is discretized by piecewise linear conforming finite elements. We introduce two kinds of indicators, both of them of residual type. The first one is related to time discretization and is local with respect to the time discretization: thus, at each time, it provides an appropriate information for the choice of the next time step. The second is related to space discretization and is local with respect to both the time and space variable and the idea is that at each time it is an efficient tool for mesh adaptivity. An error estimation procedure

INTRODUCTION
Numerical modeling of flow and transport in porous media is significant for many petroleum and environmental engineering problems. Most recently, modeling multiphase flow has received increasing attention in connection with the disposal of radioactive waste and for CO 2 storage sequestration in geological formations.
The long-term safety of the disposal of nuclear waste is an important issue in all countries with a significant nuclear program. One of the solutions envisaged for managing waste produced by nuclear industry is to dispose of the radioactive waste in deep geological formations chosen for their ability to delay and to attenuate possible releases of radionuclides in the biosphere. Repositories for the disposal of high-level and long-lived radioactive waste generally rely on a multi-barrier system to isolate the waste from the biosphere. The multibarrier system typically comprises the natural geological barrier provided by the repository host rock and its surroundings and an engineered barrier system, i.e. engineered materials placed within a repository, including the waste form, waste canisters, buffer materials, backfill and seals, for more details see for instance [1]. An important task of the safety assessment process is the handling of heterogeneities of the geological formation.
In this paper, we focus our attention on the numerical simulations of a single phase flow of an incompressible fluid with a dissolved radioactive solute: miscible flow, in connection with questions of safety of a nuclear waste repository. For more details on the formulation of such problems see, e.g., [2,3] and the references therein. Numerical simulation of flow and transport in porous media in petroleum and environmental applications has been a problem of interest for many years and many methods have been developed. There is an extensive literature on this subject. We will not attempt a literature review here, but merely mention a few references. We refer to the books [4][5][6] and the references therein.
In the safety assessment of deep geological repositories, a thoughtful consideration must be given to the mechanisms and possible pathways of migration of released radionuclides. However, when assessing confinement capabilities of disposal facilities and of the geological formations, the investigated domain covered by the numerical simulations includes a strong variability in the domain properties as well as in the geometrical scales (from cm to km) and presents significant computational challenges. Simulation models, if they are intended to provide realistic predictions, must accurately account for these effects. Furthermore, such flows are often characterized by localized phenomena such as steep concentration fronts. Accurately resolving these types of phenomena requires high resolution in regions where the solution is changing rapidly. For this reason, the development of some type of dynamic griding capability has long been of interest in the porous media community.
The adaptive mesh refinement has been a problem of interest for many years and many methods have been developed. There is an extensive literature on this subject for finite element approximations. Adaptive discontinuous Galerkin methods have tremendously developed in the last decade. We will not attempt a literature review here, but merely mention a few references. Here, we restrict ourself to adaptive methods for finite volume discretization in the context of flow and transport in porous media. Closely related to our work, we refer for instance to [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] and the references therein.
The French Institute for Radiological Protection and Nuclear Safety (IRSN) developed the code MELODIE, see for instance [26] and [27], and is constantly upgrading it, for numerical modeling of the physico-chemical phenomena involved in the release and in the migration of radionuclides from waste packages to the geosphere outlets (in deep geological formations, and more recently at surface or sub-surface locations). The mathematical formulation of these type of flow leads to a coupled system of partial differential equations, which includes an elliptic pressure-velocity equation and a diffusion-convection concentration equation. A vertex-centred finite volume method is employed to discretize the coupled system: the Darcy flow equation and the diffusion-convection concentration equation. The convection term is approximated with a Godunov scheme over the dual finite volume mesh, whereas the diffusion-dispersion term is discretized by piecewise linear conforming finite elements. The motivation for applying such finite volume procedure for modeling flow in porous media arises from the fact that this scheme is mass conservative element by element and satisfies a discrete maximum principle. IRSN is constantly upgrading MELODIE to improve the resolution time, the accuracy and the confidence in the results of the simulations.
The purpose of this paper is to build a posteriori error estimators with a strategy of adaptive mesh refinement, implement the method in the code MELODIE and present numerical simulations. The estimators are expected to allow managing of the process of local refinement of the mesh in order to minimize the discretization error at an optimal computational cost. We introduce two kinds of indicators, both of them of residual type. The first one is related to time discretization and is local with respect to the time discretization: thus, at each time, it provides an appropriate information for the choice of the next time step. The second is related to space discretization and is local with respect to both the time and space variable and the idea is that at each time it is an efficient tool for mesh adaptivity. An error estimation procedure evaluates where additional refinement is needed and grid generation procedures dynamically create or remove fine-grid patches as resolution requirements change. The aim is to concentrate the computational work near the regions of interest in the flow, such as in regions where the solution is changing rapidly and around wells. This can significantly reduce the computational effort required to obtain a desired level of accuracy in the simulation.
The reminder of the paper is organized as follows. In the next section, we give a short description of the mathematical and physical model used in this study. The third section describes their discretization using a vertexcentred finite volume method and a semi-implicit Euler approach is used for time discretization. In the fourth section, we present the construction of the residual error estimators and the adaptive algorithm. Results from numerical simulations are presented in the fifth section. To validate the efficiency and the accuracy of the method, two tests are investigated for 2D miscible flow problems. In the first test, a domain made of two porous media separated by an interface with one source is considered. The second test addresses the evolution of radionuclide migration from the waste packages through a complex heterogeneous disposal made of six different layers, where large permeability variations are allowed, with three source terms. Finally, we give some conclusions and remarks on this work, and discuss some of our future research.

MODEL SYSTEM FOR MISCIBLE FLOW
We consider for simplicity a two-dimensional horizontal reservoir where the gravity effects are negligible.
The single-phase flow of an incompressible fluid with a dissolved solute in a horizontal porous reservoir X & R 2 over a time period 0; T f ½, is given by Darcy's law and mass conservation (e.g., [4][5][6]): where P andq are the pressure and Darcy velocity of the fluid mixture, x and K are the porosity and the permeability of the medium, l is the viscosity of the mixture, C is the concentration of the contaminant solute, and f is a given source term. D is the diffusion-dispersion tensor given by: with E ij ðqÞ ¼ q i q j jqj 2 , 1 i; j 2, jqj is the Euclidean norm ofq, d e is the effective diffusion coefficient, and a l and a t are the magnitudes of longitudinal and transverse dispersion respectively.
The system (1-3) is subjected to appropriate boundary conditions and an initial condition. In what follows we use standard assumptions for miscible flow in porous media. Namely, we will assume that the porosity x is bounded above and below by positive constants, the permeability K and the diffusion-dispersion D are uniformly positive definite matrices, the viscosity l is such that, 0 < l À lðcÞ l þ ; c 2 ½0; 1, and d e , a l and a t are positive constants such that a t a l .

FINITE VOLUME DISCRETIZATION
For the spatial discretization of the system (1-3), we use the vertex-centered finite volume method, also called the control volume finite element method. We consider here the two dimensional case, for more details see for instance [2,28] and the references therein. The method is based on two spatial grids: a primary grid, which is a conforming finite element grid, and a secondary grid composed of control volumes centered in the vertices of the primary grid. Control volumes are constructed around grid nodes by joining the midpoints of the edges of a triangle with the barycenter of the triangle (Fig. 1). All model data, permeability, porosity, sources, and dispersivity are constant element-wise on the primary grid. The material interfaces are therefore aligned with the Before introducing the scheme, let us give some definitions and notations for the meshes. We introduce a partition of the interval ½0; T f into sub-intervals ½t nÀ1 ; t n ; 1 n N such that 0 ¼ t 0 < t 1 < ::: < t N ¼ T f : We denote by s n the length t nþ1 À t n ; and by s the maximum of s n ; 1 n N : For each n; we consider T n h a regular triangulation of X by closed triangles. We denote by V n h the dual decomposition associated to T n h : The partition V n h is chosen as the set of N control volumes V i that constitute the dual of the triangulation T n h known as the Voronoı¨mesh and such that X ¼ [ i¼1;:::;N V i : The mesh is constructed by connecting the middle points of edges and circumcenters of each neighboring pair of triangles having a common edge with a straight-line segment. Furthermore, we denote by E n h the set of edges E of the triangulation T n h ; and C n h the set of edges c of the dual decomposition V n h : In the following, we denote respectively by r 1 and r 2 the edges gm ij and gm ik and byñ r 1 andñ r 2 their outward unit normals, and by R 1 et R 2 the sub-edges of T : im ij and im ik ; andñ R 1 andñ R 2 their outward unit normals.

Discretization of the Flow Equation
To efficiently solve the governing equations, we use a sequential approach in which the pressure and concentration equations are solved consecutively. For simplicity, we will consider a permanent regime and take l ¼ 1, so the time level is omitted to shorten the notation.
Integrating Equation (2) over a control volume V i and using the divergence theorem and (1), we get: where K T is an approximation of the permeability tensor in T and i a vertex belonging to T : Denote by V T ðiÞ ¼ fj neighbor of i ; j 2 T g, the neighbors of i belonging to T , and by: the elementary diffusive term, where N j is the P1 finite element base function on the triangle associated to the vertex x j . Using a P1 Galerkin expansion of P and using the fact that P j2V T ðiÞ N j ¼ 1 in T , we obtain: Observing that oV i \ T ¼ r 1 [ r 2 , we get: We can express the flux of K T rN j on oV i \ T using the outward normals at edges of the triangle T : To do so, we use the Gaussian theorem: where V is a surface with the boundary oV : We take V as a surface limited by R 1 ; R 2 ; r 1 and r 2 : Then, we get: Taking into account that: and using (8), we get: whereñ jk is the outward unit normal at the edge ðx j x k Þ of T . We recall that rN j ¼ À jx i x k j 2jT jñ R 2 , then we get:   Finally, we incorporate the boundary conditions in (7) to obtain the linear system for the flow equation which is solved by the preconditioned conjugate gradient method. The preconditioning step utilizes incomplete Gaussian elimination.
Finally, the Darcy velocityq, performed as a postprocessing step, is approximated by a piecewise constant in each element T , denoted byq T , and given by: For each vertex x j of T ; we have rN j ¼ À jx i x k j 2jTjñ R 2 and then we can expressq T ; using the opposite edge F j of the vertex x j in T andñ j the outward unit normal vector at F j ; as following:q

Discretization of the Transport Equation
Discretization of the concentration Equation (3) is performed by using a vertex-centred finite volume method, (e.g., [2]), with a semi-implicit time stepping, the time discretization is implicit for the diffusion term and it is explicit for the convective term. The diffusion term is discretized by piecewise linear conforming finite elements, whereas the convective term is approximated with the aid of a Godunov scheme. Let us mention that a fully implicit scheme is also available in the code MELODIE. Let C n i be an approximation of C at the point ðx i ; t n Þ. First, we integrate Equation (3) over the control volume V i , we apply the divergence theorem, a semi-implicit time discretization, the mass lumping in the accumulation term and a full upwind stabilization in the convective term, resulting in the following discrete concentration equation: where The discretization for the diffusive flux in (13) is performed using the same steps as above for the pressure equation. We have: where: with D T is an approximation of the diffusion-dispersion tensor D in T .
To approximate the convective flux, we use the Godunov scheme. We have: On the dual edge r 1 separating the vertices i and j, we approximate C nÀ1 jr 1q T :ñ r 1 as following: where for r 2 R; r þ ¼ maxf0; rg and r À ¼ minf0; rg: Note that C nÀ1 jr 1 is approximated taking care of the flow direction. Denote by: the elementary convective terms. Then (17) implies: Finally, the scheme could be written in the following form: Then, we incorporate the boundary conditions to obtain, at each time step, the linear system which is solved by the preconditioned conjugate gradient method. The preconditioning step utilizes incomplete Gaussian elimination.
Let us end this section by the following remark. To ensure the discret maximum principle and hence the stability of this scheme we define locally the CFL condition as: Taking into account (18), we get the following estimation for a conv ii : where jq T :ñj max is the maximum of jq T :ñj over all dual edges belonging to T : Furthermore we get: It follows thanks to (22) and (23) that (21) holds true if: Then we define the following CFL condition:

A POSTERIORI ESTIMATORS AND ADAPTIVE ALGORITHMS
In this section, we derive an adaptive numerical technique using the finite volume approximation described in the previous section. The proposed error estimators were suggested in [8] for vertex-centered finite volumes method. Here, we extend this strategy to finite volumes finite elements discretization and for more complex geometry. The estimator was obtained by dual volumes summation of the residuals, while in the current work it is obtained by summation over the underlying triangulation. The method expresses the error in terms of the residual of the approximate solution. Based on the error estimators, the mesh adaptivity relies on mesh generator.
A new mesh must be generated at each time step. After that, there are two alternatives: either restart the computation from scratch or project all the information from the old mesh to the new mesh. The two meshes may have very different topologies and number of elements and the number of degrees of freedom can change arbitrarily to meet a prescribed accuracy. The initial mesh does not drastically influence the adaptive process, because a new mesh is rebuilt at each step. Each triangulation T nþ1 h is derived from T n h by refining some elements of T n h into a few sub-elements or by derefining some elements into a new elements. We denote D n ¼ Dðt n Þ,q n ¼qðt n Þ and f n ¼ f ðt n Þ, and we introduce D n h (resp. f n h ) a piecewise approximation of D n (resp. f n ). In the following, C n h will denote the finite volume approximation of C n ; such that C n h ¼ P N i¼1 C n i N i ; andq n h the piecewise P 0 approximation ofq n expressed by the hydraulic head.
Following [30] and [8], we define the residuals and inter-element jumps of the approximation: where T stands for the current element T 2 T n h . Let: where ½; E denotes the difference between limits from either side of the edge E 2 E n h ; andñ E is the outward unit normal vector at the edge E: Let: where c 2 C n h is an edge of the dual decomposition V n h andñ c is the outward unit normal vector at c: We define local spatial error indicators: diffusive jumps: and convective jumps: where jjvjj 0;S ¼ R S v 2 dx denotes the L 2 norm for any function v and any subset S of X: The global spatial error indicator is given by: For each n 2 f1; :::; N g, we define the temporal error indicator as: With the family ðC n h Þ 0 n N , we associate the function C hs in ½0; T which is linear on each interval ½t nÀ1 ; t n ; 1 n N ; and equal to C n h at t n ; 0 n N : This function writes, for 1 n N : We define the norm jjC h;s jj ðt nÀ1 ;t n Þ as follows: We define the global error indicator as follows: Then for an approximate solution C h;s ; we define the relative error estimator g r by: Our aim is to control the relative error between the exact and approximated solution by a prescribed tolerance. In the following, we give an overview of the basis time-stepping routine and then describe in detail each time step of the algorithm. We denote by N T nÀ1 h the number of triangles belonging to T nÀ1 h : We present here an adaptive algorithm based on the above a posteriori error estimates which is designed to ensure that the relative energy error between the exact and approximate solutions will be below a prescribed tolerance. The spatial refinement-derefinement is performed using the maximum strategy, while the temporal adaption is based on the temporal error indicator and the CFL condition. Let e a predefined parameter. The aim is to build an adaptive mesh T n h ; n ¼ 1; :::; N for an adapted time step s n such that the error indicators ensure that the relative error estimator is below the prescribed parameter e; such that g r e 2 so as to equilibrate the space and the time estimators.
For a given time level t nÀ1 , we set: Tol ¼ e jju hs jj t nÀ1 ;t n ð Þ ffiffi ffi 2 p We denote by 0 < d ref < 1 and 0 < d deref < 1 two parameters. For practical implementation purposes and because computer limitations, we introduce maximal level parameters N sp ; Fr sp and Fr tm ; where: -N sp , is the limit number (level) of refinement, -Fr sp , the frequence of using the spatial refinementderefinement process, -Fr tm , the frequence of using the temporal adaptation process, -Cfl lim , a prescribed value constraint that the CFL condition limit does not exceed.

Adaptive algorithm
Let an initial mesh T 0 h , and an initial time step s 0 Set n = 1 and t 1 = t 0 + s 0 Time iterations: h to obtain the final mesh T n h such that the refinement'slevel is less than N sp While g n sp ! Tol Depending on the value fr tm If g n tm > Tol and CFL < Cfl lim Set t n ¼ t n À s nÀ1 and s nÀ1 ¼ s nÀ1 =2 Else Set t n ¼ t n þ s nÀ1 and s nÀ1 ¼ 1:1s nÀ1 Save the approximation C n h ; the mesh T n h and the temporal step size s n : Set n ¼ n þ 1: The mesh refinement algorithm is based on bisecting an element into three triangles. Suppose that a set of elements are scheduled for refinement, then those elements are bisected by creating a gravity center node, see Figure 2. The derefinement process consists in deleting every couple of triangles, based on error indication information. If so the mesh would be coarsened to the one shown in Figure 3 or Figure 4. Suppose that error indicators reveal that these elements may be coarsened. The strategy of collapsing the vertices of the current triangle onto his gravity center is done by deleting all edges connected to those vertices. After coarsening or refining we update the mesh calling a Delaunay correction program to ensure that the mesh respects the empty circle criterion (Fig. 5, 6). The shape of elements containing small or large angles that were created during refinement or coarsening may be improved by edge swapping. This procedure operates on pairs of triangles T 1 and T 2 ; which violate the empty circle criterion, that share a common edge E. The edge swapping occurs deleting this edge and the empty circle criterion is satisfied for the triangles K 1 and K 2 Figure 7. The process is continued until all irregular triangles are removed. As such, swapping will have to provide mesh quality.
Remark 2. In addition of the set of elements that are marked in this way, we also mark all their neighbors with a common edge. We moreover add to this set those elements that now have more than half of the neighbors marked (Fig. 3).  Refinement process of the current triangle.

Figure 3
Coarsening of a polygonal region. Remark 3. In order to approximate the flux at the fine grid in the boundaries we introduce additional (auxiliary) points, in which the solution is obtained by interpolation on the coarse grid (Fig. 8).
We end up this section by the following remark about the interpolation of the solution after a modification of the mesh.
Remark 4. Note that the solution interpolation is also a key point in the mesh adaptation algorithm. For that, we need an interpolation scheme to transfer the information of the solution from the current mesh to the newly adapted mesh. Before proceeding to mesh adaption, we save the history of the current mesh. The new generated vertices are located in the current mesh by identifying the elements containing them. Then an interpolation scheme is used to extract the information from the solution. More precisely, the mesh refinement algorithm is based on bisecting an element into three triangles by adding a gravity center node. Then we assign to this new node the average of the values of the solution to the vertices of the triangle containing it. The temporal error indicators are used to control the errors due to the solution interpolation stage. The software automatically selects time steps to satisfy a prescribed (local) temporal error tolerance and to maintain stability.

NUMERICAL RESULTS
In this section, we illustrate the adaptive methodology. We study the performance of the algorithm by looking at the spatial adaptation of the mesh and aim at obtaining a good distribution of nodes, in taking into account also the temporal adaption. The following numerical two examples are performed to demonstrate the implementation of the proposed algorithm. In the first example, we verify our underlying adaptive numerical algorithm. In the second one, we present a numerical example for a waste-disposal in a complex domain.

Example 1: Heterogeneous Domain with One Source
This scenario consists of a pulse of activity in a heterogeneous domain. This example helps to get the first conclusions about the efficiency of the adaptive method. We consider the domain X ¼0; 50½Â0; 20½ composed of two parts, i.e., the left medium X l ¼0; 20½Â0; 20½ and the right medium X r ¼20; 50½Â0; 20½ (Fig. 9). The values of physical parameters are presented in Table 1. The dispersivity is isotropic for the right medium and anisotropic for the left. The flow calculation is carried out with a condition of constant pressure head on the right vertical limit of a value of 30 meters and on the left vertical limit of a value of 0 meter. The initial concentration is 10 4 Bq at the center of the left medium and releases immediately at the first step of time of the simulation. The transport equation is solved for a short time interval (0-5 years) with adapted time step size.
For this domain, the adaptive strategy is applied. The error indicators are computed for each element T at time The edge swap mesh modification.   Heterogeneous domain: X l (left) and X r (right).  plume and wide in the right medium after the interface between the two media. The mesh follows the movement of the plume. This numerical test shows that the proposed adaptive strategy is effective for guiding mesh adaption.

Example 2: Waste-Case Model with Different Sources
One of goals of this work is to simulate the flow of water and transport for 10 000 years in realistic conditions. The model is a cross cutting geological of 5 000 meters long and 600 meters high with 3 different layers cutting by 3 faults cf. Figure 11. The different values of parameters are presented in Table 2. The groundwater flow is under steady state condition. The flow calculation is carried out with a condition of the pressure head equal to the topographic level on the top of the model by applying a first spatial adaptation of the a posteriori error method.
On the field of hydraulic head obtained, the simulation of the immediate release of the three sources of waste is modeled with a condition that the adapt mesh cannot be less precise than the first mesh on which the water flow calculation is made. The simulation time is 10 000 years and the calculation of the transport of the contaminant is made upon all the domains. Figure 12 displays the state of the concentration in the cross cutting for 3 dates of the simulation and the corresponding propagation of the mesh like for the example 1. In this case, the adaptive strategy also performs fairly well with mesh refinement which is made almost solely around sources points.

CONCLUSION
In this paper, we presented an adaptive mesh refinement algorithm for single phase flow in porous media. The processes modeled are miscible flow of an incompressible fluid with a dissolved radioactive solute to simulate the evolution of radionuclide migration in a nuclear waste repository. The resulting system is discretized by a vertex-centred finite volume method. The method was implemented in the software MELODIE. This leads to a considerable improvement in computational efficiency. The adaptive strategy has been illustrated by means of numerical examples in an academic scenario and a realistic scenario. The above results illustrate that the proposed adaptive mesh refinement is capable of tackling in a robust and accurate fashion various physical phenomena relevant to flow and transport of radionuclides in heterogeneous porous media. Our goal in future work will be to extend this approach to more realistic 3D problems.