The hydrological team at LSCE (http://www.lsce.ipsl.fr/) has been active in the field of hydrological & hydrogeological modeling with applications in heterogeneous porous media, fractured media, coupled transfers. See some references below.

The Cast3M code (http://www-cast3m.cea.fr/)

Cast3M is a multi-physics code dealing with various applications, initially developed with a finite element scheme for nuclear reactor applications, including solid and structure mechanics, fluid and heat transfers. It consists of various elemental bricks called procedures that can be organized together for the resolution of more complex problems or equations. The resolution of equations dealing with transfers in porous media has been developed from the 90s on. Cast3M now provides tools to resolve a large range of problems including saturated flows, unsaturated flows (Richard’s equation and multi-phase flow), Eulerian and Lagrangian transport by means of Finite Volume and Mixed Hybrid Finite Element schemes. The latter proved accurate and efficient for nuclear waste storage applications (flow and transport) within an inter-comparison exercise (Bernard-Michel, 2004) making Cast3M a reference code for applications in the transfer in porous media. Several extensions including these procedures were developed in the past taking advantage of the modular properties of the code, for instance to simulate couplings of hydrological processes with mechanics and/or heat transfers, or to solve a set of equations coupling surface and sub-surface transfers (Weil et al, 2009), or specialize the simulation tools for applications dealing with fractured media (e.g. Grenier et al, 2009).

Numerical approach for the coupled set of Thermo-Hydrological equations within the benchmark

The procedure solving coupled TH transfers with phase change was developed based on the core procedure TRANGEOL solving time step Eulerian transport (diffusive-dispersive & advective transport). TRANGEOL provides a variety of options: a VF (Le Potier 2004, 2005, 2010) or a MHFE (Mosé et al. 1994; Dabbene et al. 1998) numerical scheme, theta-methods (from implicit to explicit) for diffusion and advection, as well as various options for system matrix inversion and conditioning (www-cast3m.cea.fr). It was largely tested for nuclear waste storage applications showing the practical pros and cons of most of these options and approaches.

To resolve the TH benchmark cases the TRANGEOL procedure was used with the implementation of new development, 1°) adding the phase change term and 2°) solving the set of non-linear coupled equations by means of a Picard iteration scheme. Some points are detailed below.

Phase change resolution

The phase change term is treated as a storage term controlling the velocity of the propagation. Phase change has indeed a strong local influence on the propagation of a cooling front for instance by slowing it down as compared with the same temperature front case without phase change. This is due to the energy required to bring the elementary volume of porous medium to the temperature of phase change and then provide the energy for phase change. This leads to steep fronts. Special effort was put to stabilize the resolution with an under-relaxation scheme. The basic idea is to compute some terms of the equation which are function of the temperature (unknown for the heat equation) by their estimation based on the linear combination of temperature at present calculated time and at the previous. This approach stabilizes the convergence of the iteration but under-relaxation is generally less required after a certain amount of iteration and should be stopped to accelerate convergence. So, to optimize this approach we used the procedure by (Durbin & Delemos 2007). The basic idea is there that the under-relaxation factor depends on the progress of the iterations: the further from convergence, the larger the weight on the former time computation of the variable.

Resolution of conduction & advection terms

The procedure developed for the complete set of equations relies upon the elementary time step resolution of the diffusive & advective equation provided in TRANGEOL. It offers MHFE and VF schemes, a theta method (from implicit to explicit) for diffusion and advection. These numerical schemes were originally designed for the resolution of the conductive heat transfer (equivalently diffusive transport). The extensions to advection were made expressing the advection as a source term for MHFE (Dabbene et al 1998) and using a modified version of the method by () for the FV (Le Potier 2004). The stability and accuracy of the resolution has to be controlled with the evaluation of Fourier and CFL numbers respectively. For high Peclet numbers (strong advection), some options are provided in TRANGEOL to obtain an automatic unconditional stability of the resolution. The procedure behind it is an upwind scheme, formally equivalent to the introduction of a dispersion coefficient equal to mesh size. This leads to the CFL value of 1 in the corresponding meshes.

Resolution of the coupled set of equations

The resolution of the coupled set of equations is achieved iteratively, with a Picard scheme: the properties for both equations are computed from the values of the unknown at former iteration until convergence of both unknowns. The convergence criterion is expressed as the maximum value of the relative discrepancy between the two last computed pressure and temperature fields. For each time step, the convergence of the iteration scheme was considered achieved when the value of former criterion fell below a threshold. A threshold of 10-4 is commonly considered.

Some references

Bernard-Michel, G., Le Potier, C., Beccantini, A., Gounand, S., Chraibi M. : The Andra Couplex 1 test case: comparisons between finite-element, mixed hybrid finite element and finite volume element discretizations. Comput Geosci 8:187–201, 2004

Dabbene, F., Paillere, H., Magnaud, J.P.: Mixed hybrid finite elements for transport of pollutants by undergound water. Proc. 10th Conference on Finite Elements in Fluids, Tucson, AZ, January 1998

Durbin, T., Delemos, D.: Adaptative underrrelaxation of Picard iterations in ground water models. Ground Water, 45(5):648-651, 2007

Grenier, C., Bernard-Michel, G., Benabderrahmane, H.: Evaluation of retention properties of a semi-synthetic fractured block from modelling at performance assessment time scales (Äspö Hard Rock Laboratory, Sweden). Hydrogeology Journal 17: 1051–1066, 2009

Grenier, C., Régnier, D., Mouche, E., Benabderrahmane, H., Costard, F., Davy, P.: Impact of permafrost development on underground flow patterns: a numerical study considering freezing cycles on a two dimensional vertical cut through a generic river-plain system. Hydrogeology Journal, 2013, Volume 21, Issue 1, pp 257-270

Le Potier, C.: A nonlinear correction and maximum principle for diffusion operators discretized using cell-centered finite volume schemes. C. R. Acad. Sci. Paris, Ser. I, 348 (2010) 691-695, 2010.

Le Potier, C. : Schémas volumes finis pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés [Finite volume schemes for highly anisotropic diffusion operators on non-structured triangular meshes]. C R Acad Sci 341(12):787–792, 2005.

Le Potier, C.: Finite volume scheme in two or three dimensions for a diffusion-convection equation applied to porous media with CASTEM2000. Dev. Water Sci. 55(2), 1015-102,6 2004.

Mosé, R., Siegel, P., Ackerer, P., Chavent, G.: Application of the mixed hybrid finite element approximation in a groundwater flow model: luxury or necessity? Water Resour Res 30(11):3001–3012, 1994

Régnier D.: Modélisation physique et numérique de la dynamique d’un pergélisol au cours d’un cycle climatique. Implications pour le site Meuse / Haute-Marne. Mémoire de thèse, Université Rennes I, 2012.

Weill, S., Mouche, E., Patin, J.: A generalized Richards equation for surface/subsurface flow modeling. Journal of Hydrology 366, pp. 9-20, 2009