A need for numerical simulations

Power generation plays an important role in global warming (Audoly et al. 2018) and although the use of nuclear power in the energetic mix can be seen as sustainable (Brook et al. 2014; Knapp and Pevec 2018), nuclear energy is debated in many countries (Meserve 2004; Cici et al. 2012). Over 50 years of nuclear energy and the use of radioactive material in nuclear research and in industrial, medical and other applications have left a legacy of different kinds of nuclear waste awaiting final disposal worldwide. To ensure very long-term isolation to protect the environment and ensure the safety of the future generations (Linsley and Fattah 1994; Hummel and Schneider 2005), geological disposal facilities (‘repositories’) are considered to be the most suitable solution. When the activity in waste is relatively low and half-lives are less than about 30 a, i.e. for intermediate- and low-level waste, near-surface disposal is often considered to be adequate. For waste with a higher activity and/or longer life, i.e. for high-level long-lived waste (HLLW), the concept of geological disposal of radioactive waste emerged in the late 1950s (Hess 1957; De Marsily et al. 1977; Apted and Ahn 2010; Chapman and Hooper 2012). Since then, many international research and development programs have been launched to study deep argillaceous formations, granitic rocks and salt formations as potential host rocks for radioactive waste disposal (Landais and Aranyossy 2011)(Fig. 1).

Based on their outcomes, there is a broad technical consensus that geologic disposal will meet the safety requirement to minimize safety and environmental impacts, now and far into the future (Grambow and Bretesché 2014), meaning a hundred thousand to million years into the future. Demonstrating safety over time comparable to geological timescales relies on a rigorous, complex and iterative scientific approach referred to as “long-term behavior science” that is based upon three pillars: experiments, modeling, and natural/archaeological analogs (Poinssot and Gin 2012; Dillmann et al. 2014; Alexander et al. 2015; Martin et al. 2016). Today, numerical models are ubiquitous tools that can be used either for long-term predictive evaluation of solute transport or design optimization of nuclear waste geologic repositories (Tsang et al. 1994; Johnson et al. 2017). They are also used to make predictive multi-physical assessments within a timeframe and space scale larger than experiments can cover (Bredehoeft 2003; Bildstein and Claret 2015). These numerical simulations require integrating, in a consistent framework, an increasing amount of scientific knowledge (Geckeis and Rabung 2008) acquired for each of the individual components of such repositories. This implies considering couplings of different non-linear processes, applied to a wide range of materials with contrasting properties as a function of time and space in ever-larger systems. Reactive transport modeling (Yeh and Tripathi 1989) has established itself as a powerful, versatile and essential tool for tackling the complexity of the multi-space and temporal-scale issues engineers and scientists face. This includes describing the evolution of the components constituting the so-called “multiple-barrier system” between the waste matrix and the biosphere (Apted and Ahn 2010; Chapman and Hooper 2012). Predicting how (i) the waste matrices (e.g., glass, bitumen, cement), (ii) waste overpacks (e.g., metal canisters, concrete), (iii) engineered barriers such as bentonite (Sellin and Leupin 2013), (iv) and natural geological barriers will evolve with time in response to physical and chemical perturbations is of prime importance for performance and safety evaluations of the repository concepts. This chapter summarizes recent improvements and discusses future challenges in the application of reactive transport to deep geological nuclear waste disposal focusing on disposal in clay formations.

A need to consider coupled processes

The long-term safety of geological disposal is based on the multi-barrier concept consisting of a combination of natural (host rocks) and engineered (waste form, package, backfill and seal materials) barriers that will have mutual interactions during the lifetime of the repository. Predicting the effects of these interactions entails understanding, evaluating and prioritizing the pertinent thermal, hydraulic, mechanical, chemical and radiological (THMCR) processes (Fig. 2).

Thermal effects arise principally from heat generated by waste that is both dissipated in the over-packs and the geological medium, and evacuated by air ventilation (Benet et al. 2014a). The hydraulic stage and gas-related effects are a combination of repository resaturation and of gas generation. (e.g., H2 production by anaerobic corrosion of metals or through water radiolysis). The hydraulic transient is also impacted upon by the ventilation systems used during the operation period to ensure adequate gallery aeration for the workers and for the infrastructure. It will create heat and vapor exchanges with walls and storage packages, pressure variations and evaporation, and sometimes vapor condensation (Benet et al. 2014b). Mechanical effects can be induced during construction and excavation of underground drifts that will cause both tunnel convergence (Lisjak et al. 2015) and damage to the rock in the vicinity of the opening, with the formation of an associated excavated disturbed zone (EDZ) (Armand et al. 2014). There are also numerous chemical effects as foreign materials like borosilicate glasses (Poinssot et al. 2010; Gin et al. 2015), metallic canisters (King and Shoesmith 2010), and concrete (Alonso et al. 2010) that are introduced into the repository will induce chemical gradients across the repository components (Nagra 2002; Andra 2005, 2009). Because of these chemical gradients, perturbations such as pH and redox changes may alter the performance of the barriers over time (Bildstein and Claret 2015). Last but not least, radiological effects also exist in association with the radiological inventory of waste independently from the technological solution, retreatment or direct disposal of irradiated fuels (Odorowski et al. 2017) or after reprocessing (Gong et al. 1999; Berner et al. 2013). Moreover, an additional challenge arises because all of the THMC phenomena described above are coupled. The coupling of physical and chemical phenomena can be either weak or strong and may vary in time and space scales. As an illustration, excavation induces a mechanical rock effect producing an EDZ. The change in rock microstructure in the EDZ compared to pristine rock might impact radionuclide (RN) transport parameters if this damage zone acts as a preferential pathway for fluid flow and RN transport. However, issues are more complicated as the EDZ can heal and seal itself depending on hydraulic conditions (Thatcher et al. 2016). Another example of a coupled process is the evolution of redox conditions, a key parameter for RN solubility (Duro et al. 2014) which can be modified in connection with hydraulic and mechanical process. The oxidation event that will also occur in the repository near field (Matray et al. 2007; Craen et al. 2008; Vinsot et al. 2013; De Windt et al. 2014; Vinsot et al. 2014) as result of excavation and ventilation will change redox conditions. These conditions could also been modified if radiological effects are considered. The radiolysis effects of water under alpha irradiation simultaneously producing oxidizing species (e.g., hydrogen peroxide), and reducing species like hydrogen might play a role when competing with redox active species from the ambient environment (Odorowski et al. 2017). Many other examples of features and processes and their relevance in performance assessment and safety cases are given by (Bernier et al. 2017), including microbiological effects. These combined processes may be detrimental or favorable to the overall performance of a repository system over time. Since the combinations will occur on time-scales that are not accessible to experimentation, we must develop modeling approaches that can help to predict how barriers will evolve in time and space (Claret et al. 2018b). Reactive transport modeling is probably one of the most efficient techniques used to account for and quantify the complexity of coupling processes over a long period of time. Many codes exist (see Steefel et al. 2015b for a recent review) and to our knowledge, none of them are able to integrate all coupled THMCR effects. However, since the 34th volume of Reviews in Mineralogy (dedicated to reactive transport in porous media) was published about 20 years ago, numerical codes have continuously improved and are capable of representing more and more complex situations. Active topics of research dealing with RTM include the development of pore scale and hybrid, or multiple continua, models to capture the scale dependence of coupled reactive transport processes (Steefel et al. 2005). In this chapter, most of the reviewed studies that are described are based on continuum models (Lichtner 1996), although other kinds of approaches will be addressed at the end.

Governing equations

The reactive transport constitutive equations rely on the description of the porous medium of interest (e.g., one of the component of the multi-barrier system) at the continuum scale, with respect to its macroscopically measurable properties such as permeability, dispersivity and diffusability (Steefel et al. 2014, 2015b). In this framework, a generic equation for advective/dispersive transport in the liquid phase coupled to bio-geochemical reactions can be written:


where ϕ is the porosity (m3void·m−3medium), SL is the liquid saturation (unitless), Øi (mol·m−3water) is the total concentration term that integrates the distribution between primary species (with concentration Ci) and secondary species (Σj=1NxvjiCj),Di* is the diffusion coefficient for species i in the porous media (m2·s−1), q is the volumetric (or Darcy flux) of water (m3water·m−2medium·s−1), and Rr, Rm, Rs, and Rl are the aqueous phase, mineral, surface, and gas reactions (mol·m3medium·s−1) respectively. vji is the number of moles of component i in one mole of secondary species j. vik with k = r, m or l is the number of moles of component i in one mole of phase k.

These equations can be expanded to consider (i) transport in phases other than the liquid (e.g., gases), (ii) a dispersion term, or (iii) the Nernst–Planck Equation, instead of Fick's Equation, to model diffusional transport. A comprehensive description of these equations can be found in Steefel et al. (2015b). Recently, the significance of the parameters in Equations (1) and (2), their relationship in case of mutual dependence (e.g., diffusion and saturation permeability and saturation) and the way they can be approximated using different empirical or based on theory laws have been reviewed (Claret et al. 2018b). Therefore, it will not be detailed again here. Equation (1) also highlights the fact that a porosity has to be defined for the considered materials. However, some of them (e.g., steel and other metals, glass) are not porous materials, and, unfortunately, reactive transport models are, by definition and according to Equation (1), not suited to include non-porous media (Bildstein et al. 2007; Claret et al. 2018b). To unravel this problem, the volume clearance that exists at the interface between the canister or the glass and the surrounding material, also known as the “technological gap”, is often included in the numerical cells that represent these materials.

The critical need for thermodynamic databases

Inherently, to model reactivity and couple it with transport, geochemical databases are needed. To illustrate this point, let us consider the transport of contaminants in a repository and its relationship with geochemistry. Speciation of the considered element will be a key parameter controlling the species distribution among the different phases and their solubility limits. While the solute species can diffuse with the porewater through the barriers, their precipitation represents an appreciable retardation mechanism (Wanner 2007). Speciation also controls sorption processes that can induce retardation mechanism for transport (Grenthe 1991; Grambow 2008). Nonetheless, it is not only a question of speciation. Transport is also governed by porosity change, and this is why thermodynamic databases should also account for minerals that constitute the barriers and secondary minerals that are able to precipitate when different barriers are in contact. In the field of nuclear energy, database comprehensiveness and consistency are key issues that have been initiated by Nuclear Energy Agencies and later prioritized at the national level (see Ragoussi and Costa 2019 for a review). In addition, the pivotal importance of databases for RTM calculation has recently been the focus of a special issue of Applied Geochemistry (Kulik et al. 2015). Within the framework of national radioactive waste disposal projects, database development has focused on consistent data for actinides, lanthanides, and chemotoxics, among other potential contaminants. When modeling disposal conditions, there is also a need to describe the evolution of the materials constituting the barriers. This research effort combines new thermodynamic data acquisition, especially for clay minerals (Gailhanou et al. 2007, 2009, 2012, 2013, 2017; Blanc et al. 2015a) and cement phases (Roosz et al. 2018), critical data selection for repository materials such as cement phases (Blanc et al. 2010; Walker et al. 2016; Lothenbach et al. 2019) and how to evaluate their capacity to describe interactions between materials (Blanc et al. 2015a). In addition, kinetic databases with kinetic rate parameters are being built (Marty et al. 2015a), as well as a sorption database (Brendler et al. 2003). An important constraint for kinetic rate parameters is consistency with the associated thermodynamic database for stoichiometry and solubility values, so that we can calculate how rates depend on saturation states.

Modeling materials and reactive interfaces

Studies on deep geological storage show that most of the physical and chemical reactivity will be concentrated close to the interface between the different materials. A first challenge therefore consists of predicting the extent of the alteration and the distance at which the perturbations will progress into the barriers and, in case of failure, into the waste matrices. Another challenge is to predict the evolution of the material properties in these barriers resulting from their alteration, and to assess the impact on the overall performance and safety of the system.

A non-exhaustive list of interfaces is given in Table 1 with a description of their occurrence in different disposal concepts including clay formations (e.g., in Belgium, France, Spain, and Switzerland) and crystalline host rocks (note that bentonite is used as a barrier between the containers and granite in this case, e.g., in Canada, Spain, and Switzerland). Another degree of complexity is added in HLLW repository concepts from Sweden and Finland where the use of copper containers is envisaged, as well as in France where a bentonitic cement may be emplaced at the extrados of the steel liner to avoid early corrosion in acidic conditions.

Concerning waste matrices, only the specific behavior of vitrified waste (glass) will be detailed in this chapter. For the alteration of spent fuel (SF), the authors refer the reader to the recent review of these aspects by (Ewing 2015). Waste matrices for ILLW per se are not detailed in this chapter, but most of the conclusions will be applicable to cementitious and metallic waste types. For bituminous and organic waste, the authors refer the reader to the modeling work of Sercombe et al. (2006), von Schenck and Källström (2014), and De Windt et al. (2015).

The direct interaction of materials with groundwater is not treated specifically in this chapter, even though the authors acknowledge that many interesting studies related to waste disposal in fractured hard rocks, where this issue is particularly relevant, have been conducted. In the disposal concepts that rely upon a chemical buffer of high pH (mainly ILLW), the issue of the interaction of groundwater with cement/concrete has been addressed (see (Wilson et al. 2018) for a recent review). RTM studies of cement-groundwater interactions have been reported with results concerning pH evolution, and chemical and mechanical degradation of concrete barriers (Höglund 2001, 2013; Bamforth et al. 2012; Cronstrand 2014; Wilson et al. 2018). In HLW disposal concepts, the bentonite buffer surrounding waste packages will evolve with time due to progressive reaction with saturating ambient groundwater, including ion exchange, dissolution-precipitation of accessory (gypsum, halite quartz, calcite) and clay minerals potentially affecting properties such as the hydraulic conductivity and the swelling pressure. RTM studies regarding these aspects have evolved from simple exchange models to more complex 1- and 2-D problem setup, including the effect of advection in a fracture intersecting the disposal cell (e.g., Arcos et al. 2003, 2008; Sena et al. 2010; Benbow et al. 2019).

Perturbations from atmospheric O2 and CO2 during the construction and ventilation stages of the repository life constitute a particular type of interface, a physical contact with air, where the aggressive agent is a gaseous component instead of interstitial water and minerals.

The modeling of the behavior of each type of material and the relevant interfaces (Table 1) will be detailed in this chapter with a particular emphasis on the evolution of the RTM approaches and achievements over the past two decades.

This “interface” is characterized by the direct physical contact between materials and air containing reactive gases (essentially O2 and CO2). In this situation, porous media are generally partially saturated with water and therefore fast diffusional transport of gas (four orders of magnitude faster than in water) can possibly occur into the pore network. The water saturation conditions will change with time resulting from the resaturation process with water coming from the surrounding host-rock and invading the tunnels, shafts, and wells. Resaturation will be initially impeded by ventilation during the operation phase and then by the generation of gas that may accumulate in the disposal cells (mainly H2) after closure. These two stages will produce very different geochemical environments, with contrasting reactivity, during which oxic vs. anoxic redox conditions will prevail. Long timespans with oxic conditions are to be avoided to limit the corrosion rate and to make sure reduced conditions are settled when radionuclides are released (because their oxidized forms tend to be more mobile). Calculating the duration of this stage is therefore important, but it is not thought to be very long once the disposal cell is closed (tens to thousands of years) due to the numerous oxygen consuming reactions that will occur in the vicinity of the disposal cell such as iron corrosion, dissolution of reducing mineral, and microbial respiration (see following sections).

The duration of the resaturation stage is difficult to predict because it strongly depends on the relative importance of the two counteracting processes (Marschall et al. 2005). It has been studied by conducting numerical simulations at disposal cell scale where the focus was to predict the transport and fate of hydrogen. The resaturation time calculated in the near field for different repository concepts varies in a very broad range from hundreds to hundreds of thousand years. It was found to be strongly linked to the diffusivity of dissolved hydrogen and the H2 production rate in the HLLW cells (Xu et al. 2008; Poller et al. 2011; Senger et al. 2011; Enssle et al. 2014; Brommundt et al. 2014; Sedighi et al. 2015) and was even longer in the ILLW cells where the degradation of waste matrices is also susceptible to produce H2 (Talandier et al. 2006; Poller et al. 2011; Avis et al. 2013). Note that flow and reactive transport in unsaturated conditions have also been extensively studied and modeled, though in a very different geological context and repository concept, in the framework of the Yucca Mountain project (USA); in this concept, the HLLW repository was envisaged in a fractured volcanic tuff formation at a depth approximately 300 m above the water table (e.g., Glassley et al. 2003; Spycher et al. 2003).

Corrosion of carbon steel components in oxic conditions

In repository conditions, metallic iron is not thermodynamically stable in the presence of water, and the corrosion of metallic iron occurs both in oxic (operation phase) and reduced conditions (post-closure period). In the first case, corrosion leads to the release of ferric iron and hydroxide in solution (Reaction 3), whereas in the second case corrosion produces aqueous ferrous iron, hydrogen and hydroxide (Reaction 4).

4 Fe(s)+6 H2O+3 O2(aq)4 Fe3++12 OH

This reaction leads either to a pH increase with an increase in redox potential or to a pH increase with a decrease in redox potential. While anoxic corrosion is often taken into account in RTM studies (see dedicated section below), oxic corrosion is seldom tackled, despite a much higher corrosion rate (Féron et al. 2008; Johnson and King 2008) and the development of more aggressive mechanisms such as pitting or crevice corrosion (Barnichon et al. 2018). This phenomenon is also observed in experiments with sequential aerobic–anaerobic conditions (Sherar et al. 2011; El Hajj et al. 2013).

There are few modeling studies dealing with corrosion during the oxic phase in the context of a disposal cell. The first attempts were dedicated to the assessment of the evolution of O2 and H2 content during the ventilation period using a 2D domain containing the container-EDZ-argillites system (Hoch and Wendling 2011). They used a TH(C) model where the only reaction is iron corrosion, and the model uses a water saturation-dependent formulation for the corrosion rate in oxic (consuming O2) and anoxic (producing H2) conditions. The results showed that oxygen diffusing from the tunnel is rapidly consumed by corrosion (with a rate ~100 µm/y) and the propagation front was limited to ~10 m (in a 40 m disposal cell). With the production of H2 on the other side (corrosion rate of 10 µm/y), the mole fraction increased to 0.3 at the end of the cell, but the fraction of both gases reached only ~0.1 where they mixed.

In their study, De Windt et al. (2014) added the kinetics of pyrite oxidative dissolution as an oxygen consuming reaction. They first checked the robustness of their RTM approach to reproduce an in situ experiment conducted in a borehole in the Tournemire URL (France). The predicted mineralogical evolution (goethite/magnetite) qualitatively matched the in situ experiment. Applying the model to the repository scale, they show that the gas diffusion coefficient in the partially saturated zones plays a major role in the location of the oxidizing/reducing front inside the engineered barriers with oxic corrosion of steel on one side and anoxic steel corrosion accompanied by hydrogen production on the other side. A similar study was performed by (Bond et al. 2013), matching the sulfate concentration from pyrite dissolution in an in situ experiment in the Mont-Terri underground research laboratory (Switzerland).

Recent studies tend to integrate more phenomenological models to determine the corrosion rate in oxic and anoxic conditions. The oxic corrosion model takes into account the thickness of the water film on the iron surface, and the diffusion of O2 in a goethite/lepidocrocite layer as the rate-limiting step (Hoerlé et al. 2004). The onset of corrosion happens when a critical amount of water condenses on the iron surface. The results of calculations at the disposal cell scale show in particular that corrosion concentrates in locations where the relative humidity is highest, especially where water “leaks” from fractures, with corrosion rates up to 200 µm/y locally.

Oxic transient: impact on claystone

Argillaceous rocks react with meteoric carbon dioxide (CO2) and oxygen (O2) from the atmosphere and nested chemical reaction fronts form in the subsurface in response to acid–base and redox reactions (Brantley et al. 2013). In a geological repository, similar phenomena will occur in the anaerobic host rock during construction (excavation, drilling operations) and operations (gallery ventilation). Under these conditions, the prevailing reducing conditions will be perturbed and redox-sensitive minerals (e.g., Fe-bearing minerals) will react. The expected chemical changes have been recently reviewed (Bildstein and Claret 2015) and the main weathering effects to be expected are recalled here. In claystone, oxidation mainly affects pyrite, which is ubiquitously found in the mineral assemblage. As it oxidizes with exposure to oxygen and water, pyrite releases sulfate and protons. The overall reaction under oxic conditions is expressed as:

FeS2+0.5 H2O+3.75 O2(aq)Fe3++2 SO42+H+

While the release of ferric ions causes the precipitation of iron (oxy)hydroxides, carbonates will play a major role in the buffering of proton release concomitantly affecting the population of the clay exchanger. At the laboratory scale, flow-through oxidation experiments conducted on COx tailings can be reproduced using RTM studies (Claret et al. 2018a). Recently an in situ experiment (OXITRAN) has been set up in the Tournemire underground research laboratory in which the time evolution of oxygen partial pressure in a measurement chamber isolated from the atmosphere has been recorded (Barnichon et al. 2018). Although pyrite is present and can buffer an oxygen plume while CO2 is produced (Vinsot et al. 2017), oxygen was never completely depleted in the test chambers. The key controlling parameters are the thickness of the Excavated Disturbed Zone and the ratio between the effective diffusion coefficient and the oxygen consumption first-order effective rate. In contrast, the ‘Full-Scale Emplacement’ (FE) experiment that was initiated at the Mont Terri rock laboratory aims to simulate, as realistically as possible, the construction, waste emplacement, backfilling and early post-closure evolution of a spent fuel/vitrified high-level waste disposal tunnel according to the Swiss repository concept. First results indicate rapid oxygen consumption at locations not affected by O2 inflow from the access tunnel (Müller et al. 2018). In addition to oxygen diffusion, drilling is also associated with desaturation and water evaporation, leading to increased salt concentration in the pore water (Zheng et al. 2008; Lerouge et al. 2013). A multiphase flow and reactive transport model of a ventilation experiment performed on Opalinus Clay indicates that changes in the clay mineral porosity caused by oxidation and the associated mineral dissolution/precipitation may seem weak (Zheng et al. 2008). The pore water seeping into the drifts and the gallery (e.g., after closure of the repository) will however interact with the oxidation products and the salts inherited by water evaporation. This more saline water will interact first with the repository materials.

Atmospheric concrete carbonation

This interface is relevant for concrete waste packages and cell structures (in most ILLW and the Belgian HLLW concepts) and for concrete components found in parts of other HLLW repositories (tunnels). Studies concerning this interface have benefited from the work that has been conducted for the atmospheric carbonation of buildings and civil engineering structures (see Ashraf 2016, for a review). In the carbonation process, the main chemical reaction is the acid attack caused by CO2 dissolving in water and triggering the dissolution of portlandite and the initial calcium-silicate-hydrate (C-S-H) phases to form calcium carbonate minerals (calcite, aragonite, vaterite) (Reaction 5), and new C-S-H phases with progressively decreasing Ca/Si ratio (amorphous silica being the ultimate alteration product, Reaction 6). Sulfate-bearing minerals are also affected by carbonation with the initial ettringite dissolving to form basanite (CaSO4·0.5H2O) or gypsum (CaSO4·H2O), depending on the relative humidity.


The models taking into account the coupling of CO2 diffusion in the gas phase with the drying process are of particular interest for atmospheric carbonation at the disposal cell scale. During the ventilation period, air will be forced into the tunnels at a constant, low relative humidity (~40%) at a temperature close to 25 °C (or may be slightly heated, up to 40 °C, in preliminary concepts from Andra 2005). The drying process is mainly controlled by water flow out of the concrete, the driving force being water evaporation at the surface. It can be described via full multiphase RT models or alternatively with Richards' Equation (depending on the concrete properties).

Most of the experiments conducted up to 2005 focused on carbonation “kinetics,” looking for the location of the carbonation front (identifying the pH front at ~9 with the phenolphthalein test) without providing quantitative mineralogical information. These experiments were modeled with simplified chemical models (e.g., Bary and Sellier 2004; Burkan Isgor and Razaqpur 2004). Only recently have experimental studies looked more closely at mineralogical changes, starting with portlandite and calcite, then adding aragonite and vaterite (Drouet 2010; Drouet et al. 2018), and ettringite, gypsum, and basanite (Auroy et al. 2018). The first modeling work actually using RTM was proposed by (Bary and Mügler 2006) with the objective of matching the portlandite and calcite profiles observed in experiments. They used a full multiphase flow code and developed a shrinking core model for the dissolution of portlandite and precipitation of calcite with a diffusion-limited rate. They also took into account the changes in porosity and the effect on permeability (through Kozeny–Carman relationship). The model only solves for the mass conservation of water, calcium and carbonate without full chemical treatment. Modeling results reproduce the observed carbonation front and a decrease in porosity for a series of experiments for three cement types at three different liquid saturations (0.65, 0.80, and 0.95). This model was then extended to account for the complete chemistry by Leterrier and Bary (2011), but was only calibrated on the carbonation front for the same experiment. A similar modeling approach was used by (Park 2008), although without the shrinking core model.

RTM studies of concrete carbonation at the disposal cell scale are scarce, especially those covering typical ventilation times (~100 years) and coupling carbonation with concrete drying (Trotignon et al. 2011; Thouvenot et al. 2013). In these two studies, the geometries of the waste package and the tunnel structure were simplified and handled as a 1D problem. The evolution of the complete mineralogy of concrete is investigated (see dedicated section below) and the problem is treated with full multiphase flow and transport using TOUGHREACT. It includes the feedback of porosity and phase saturation on diffusion using the Millington-Quirk relationship (with specific parameters for concrete) and a slipping factor to account for the larger gas-intrinsic permeability compared to water (Thiery et al. 2007; Zhang et al. 2015). Results show that the carbonation of concrete develops over 1 to 10 cm/100 years depending on the concrete properties and that, in the carbonated zone, the primary minerals dissolve all the way to the “ultimate” secondary minerals (calcite, amorphous silica, gypsum, gibbsite, ferric hydroxide), with no intermediate minerals (e.g., C-S-H). In contrast with previous results, porosity changes in these simulations are not significant (a few percent at most). It is noteworthy that these results are sensitive to the way the atmospheric cell in contact with concrete is implemented in codes as it has a direct impact on the flux of gaseous CO2 entering into concrete: size of the cell, diffusive properties in the gas and the aqueous phase.

The modeling results do not integrate some of the features observed in experiments. For instance, C-S-H carbonation may overrun CH carbonation (Groves et al. 1990), which is not predicted in numerical simulations. Some portlandite also remains, even after a long time in some experiments (Drouet 2010; Auroy et al. 2018; Drouet et al. 2018). Further tests with the shrinking core model may succeed in reproducing these observations. In combination with this phenomenon, the decrease of reactivity when liquid saturation is below 0.3–0.4 has been observed (Thiery 2005; Thiery et al. 2013) and implemented by Thouvenot et al. (2013) but to calibrate such a function in codes requires more data.

Another noteworthy feature of this type of simulation is that the computing time can be quite long. For instance, due to the lack of an implicit scheme for the coupling of transport in the gas phase and RT in the aqueous phase in TOUGHREACT, the CPU time can be one to several months for 100 years of physical time (Trotignon et al. 2011; Thouvenot et al. 2013). Note that a recent development proposes an interesting “look-up table approach” for the two-phase reactive transport, which could partly help solve the simulation CPU time (Huang et al. 2018).

Finally, some important effects of carbonation on the properties of concrete that have been largely observed are not actually taken into account in modeling (Czarnecki and Woyciechowski 2015; Ashraf 2016; Savija and Lukovic 2016): changes in porosity, transport properties, micro-and macro-mechanical properties, shrinkage and cracking. Experimental data has been recently acquired on changes in water retention properties and permeability (Auroy et al. 2013, 2015, 2018). The impact of carbonation is not the same for the different types of concrete (Auroy et al. 2015), which would also constitute an interesting challenge for RTM.

The iron–clay interface is studied in almost all countries where deep geological disposal of HLLW is envisaged. Significant amounts of metallic components (drift liner, structure, and containers) are used to ensure mechanical resistance of the waste package and to retard contact of waste matrices with water during the initial thermal stage, whereas clay barriers are used to retard RN migration. Although a variety of metallic alloys are envisaged in repositories (King and Shoesmith 2010; King 2014), the corrosion of these components in anoxic conditions is often simplified in the modeling studies into iron corrosion and takes place according to Reaction (7).


First modeling studies at the scale of the disposal cell

The first modeling studies concerning this interface attempted to implement the knowledge acquired in experiments on iron–clay interactions in deep geological conditions initially conducted in batch systems with powdered materials starting in the late 1990's (see Bildstein and Claret 2015 for a recent review). These experiments focused on understanding the processes, identifying the corrosion products (magnetite, siderite) and the secondary minerals resulting from the alteration of clay minerals (Fe-serpentine, and zeolites and Fechlorite at higher temperature), and quantifying the corrosion rate (in the order of 1 µm/year). The calculations presented in the study of Montes-H et al. (2005) were performed at the scale of the bentonite barrier perturbed by a source of iron on one side and in situ groundwater on the other side; mineral dissolution used the kinetics law, precipitation occurred at equilibrium, looking in particular at the montmorillonite-to-chlorite conversion.

Modeling studies then turned to a geometry closer to that of a disposal cell by including the iron canister into the calculation domain, adding a full kinetics approach, and focusing on porosity change and feedback on diffusional properties (Bildstein et al. 2006). The authors obtained the first result on porosity clogging after 5 000 years at the interface with bentonite and 15,000 years in the iron grid cell at the interface with the Callovo–Oxfordian claystone.

Both studies showed converging results in terms of pH increase, Eh decrease, and mineralogical evolution, although the presence of iron in the calculation domain produced stronger changes in the second study: pH from an initial circum-neutral value to ~10.5 and Eh from an initial −200 mV down to −800 mV during corrosion.

First modeling attempts to fit batch experiment results

Because significant uncertainties in the parameters used in the previous studies were highlighted (surface area, kinetics, evolution of the corrosion rate), the first experimental study to include modeling results came early with the objective of calibrating these parameters, at least partly, by matching the mineral paragenesis and the chemical indicators (pH, Eh, iron aqueous concentration) observed in the experiments. De Combarieu et al. (2007) used a mean corrosion rate (1.4 µm/year) measured in his batch experiments at 90 °C using powdered materials and iron foils. Modeling results matched the experimental data (pH = 10.5, Eh = −700 mV, precipitation of magnetite, Fe-serpentine, and Fe-silicate) using reactive surface areas calculated based on the particle size for primary minerals and the local equilibrium assumption for secondary minerals.

In a study looking at the diffusion of iron in bentonite, Hunter et al. (2007) were the first to introduce ion exchange and surface complexation in the modeling approach, using the corrosion rate measured in their experiments at 30 °C. Simple 1D calculations were performed to model the profile of iron measured into the bentonite; matching the amount of iron required making assumptions on the diffusion coefficient, surface complexation and magnetite precipitation.

Finally, Pena et al. (2008) introduced the first attempt to model a variable corrosion rate using a semi-analytical approach, hypothesizing that corrosion is controlled by diffusion through a growing magnetite film. The authors were able to match the data reported for corrosion in bentonite in batch experiments by (Smart et al. 2006) on a time scale of 0.5 year using a solid phase diffusion coefficient of 10−20 m2/s for the Fe2+ and 10−19 m2/s for H2, H2O, and OH species.

A decade of RTM evolution of simulations at the scale of the disposal cell

Later, a long series of modeling studies came back to the scale of the disposal cell; see review in Claret et al. (2018a) for a detailed description of the modeling set-up of the different studies.

While the previous studies focused on corrosion rate and mineral paragenesis, the subsequent studies at the disposal cell put the focus on full ion exchange and surface complexation with proton and Fe2+ (Samper et al. 2008; Wersin et al. 2008). While a 1D grid is used by Wersin et al. (2008) for the iron-bentonite system, the first 2D-axisymetric calculations are proposed by Samper et al. (2008) with an application to the iron–bentonite–granite system. Note that in both studies, the mineralogical system is simplified (no reactivity for clay minerals) and a constant corrosion rate is used in both cases, the emphasis being put on the role of sorption. The results show the importance of the surface protonation of bentonite, which reduces the pH increase during corrosion (from 11 to 9 in the bentonite), and a decrease of porosity in bentonite. The mass balance of iron shows that the precipitation of magnetite, and to a lesser extent siderite, accounts for most of the iron immobilization (iron sorbed on complexation surface is anecdotal). Note that Samper et al. (2008) introduced an interesting “progressive”, cell-by-cell, corrosion process. This is also the first study to perform a thorough sensitivity analysis of uncertain parameters (corrosion rate, protonation sites, and surface complexation vs. exchange of Fe2+).

In a series of simulations with a constant corrosion rate, Savage et al. (2010) tried to tackle the tricky issue of the long-term evolution of the Fe-minerals in the mineralogical assemblage in the iron-bentonite system by using evidence for mineral parageneses from analogous natural systems. They allowed for faster minerals to precipitate initially and other more stable minerals to take over through nucleation and Oswald ripening, based on the early work of (Steefel and Vancappellen 1990). This is the only work where the authors were able to describe a magnetite → cronstedtite → berthierine → chlorite mineralogical sequence at the iron–clay interface over a period of 1 My (according to experimental results reviewed by Mosser-Ruck et al. 2010).

It was only a couple of years after the first attempt to develop a model corrosion rate depending on the chemical environment by Peña et al. (2008), that non-constant iron corrosion rates were introduced in large-scale calculations, and to date only Wilson et al. (2015) have implemented a corrosion rate controlled by diffusion. In most other cases, a simpler model was adopted, in the form of a “standard” rate law, i.e. with a term considering the departure from equilibrium considering reaction (3). This assumption results in a progressive decrease of the corrosion rate, usually by an order of magnitude over a period of 100 000 years (Marty et al. 2010; Lu et al. 2011; Ngo et al. 2014). Note that a decrease of the corrosion rate is also achieved, to some extent, by considering that the reactive surface area depends on the amount of iron and on porosity (Bildstein et al. 2006; Savage et al. 2010; Wersin and Birgersson 2014). Interestingly, in most of these studies (except Lu et al. 2011 and Wilson et al. 2015), ion exchange and surface complexation were not included in the calculations. This choice was motivated either by the fact that the focus was put on other processes (mineralogical changes, effect of transport), or because explicit surface reactions were considered to create mass by double-counting the counter ions when the mineral-bearing surface sites dissolves away in significant amounts (see sensitivity analysis in (Bildstein et al. 2012). This problem was recently treated by creating specific thermodynamic data for surface-bearing minerals to reconcile ion exchange and dissolution/precipitation processes (Benbow et al. 2019). However, note that very few data exist for surface reactions at temperature higher than 25 °C, which undermines the robustness of the interpretation.

Non-isothermal calculations were introduced into simulations quite recently by imposing a variable temperature field calculated by a TH code (Bildstein et al. 2012) or using a full THC modeling tool (Samper et al. 2016; Mon et al. 2017). This additional feature usually tends to stiffen the calculations, because RT processes and temperature are not implicitly coupled and temperature changes affect both the reactive processes (creating a thermodynamic disequilibrium and modifying kinetics) and transport.

Table 2 provides a synoptic view of the results of these studies and shows that the most abundant corrosion product predicted by the models in the long term is magnetite, sometimes with Fe-carbonates (siderite), and Fe-silicates (greenalite), sometimes incorporating Al (berthierine, cronstedtite). Primary minerals in clay are often destabilized in favor of Fephyllosilicates or zeolites if they are allowed to precipitate. Numerical studies often differ on the precise nature of the secondary minerals. The transformation of clay minerals into Fe-chlorite, and the timing, very much depend on whether (i) it is included as a secondary mineral, in which case it is the most stable phase and precipitates from the beginning of the simulation (e.g., Marty et al. 2010, or (ii) it results from a ripening process, taking the place of a precursor mineral in Savage et al. 2010). One of the most sensitive parameter remains the corrosion rate. The extent of the perturbation is always predicted to be limited to a few centimeters, up to 20 centimeters into the clay barrier.

Concerning porosity clogging, a major assumption was made in some simulations: because the phenomenology is not known when porosity vanishes, porosity update was disabled in order to reach the end of the corrosion stage (~50,000 years) (Bildstein et al. 2012). A complete inhibition of the corrosion process has never been observed in experiments, even if a dense corrosion product layer is often identified. In addition, the ubiquity of magnetite in simulation results as the dominant corrosion product in the long term is questioned by many experimental results and archeological analogs. These issues, along with recent experimental work (Martin et al. 2008; Bourdelle et al. 2014, 2017), motivated the modelers to come back to modeling experimental results.

A recent return to modeling experimental results

Two of the most recent modeling studies revisited interactions at the experimental scale to refine the understanding and modeling of this interface. In the first work, Ngo et al. (2015) modeled iron-COx claystone interaction results obtained from the batch experiment with powdered material at 90 °C for 90 days (Bourdelle et al. 2014). Using a mean corrosion rate (determined from the experiment), standard and diffusion-controlled kinetics for mineral dissolution, and the local equilibrium assumption for precipitation, they were able to match the low pH (value ~7) and the evolution of aqueous species concentrations. They also matched the set of secondary minerals observed at the end the experiment: Ca-saponite, and greenalite. Interestingly, no magnetite was observed or simulated in these conditions. This particular feature, also observed by Bourdelle et al. (2017), was attributed to the high reactive surface area in the powdered system, producing a high precipitation rate for greenalite in the experiment (and matched by the modeling with equilibrium precipitation). Note that the simulation also predicted the transient precipitation of chukanovite, a Fe-hydroxyl-carbonate recently observed in experiments (Schlegel et al. 2010) and archeological analogs (Saheb et al. 2012), which disappears after the corrosion stage to form Ca-saponite.

In the second study, the modeling of iron corrosion in a Callovo–Oxfordian claystone block (machined from a core plug) at 90 °C for 2 years (Martin et al. 2008) was conducted by Bildstein et al. (2016). The objective was to reproduce the mineralogical paragenesis identified by (Schlegel et al. 2014) using the evolution of the corrosion rate measured in the experiment using electrochemical impedance spectroscopy (from 100 to 0.1 µm/year). The sequence of minerals observed from the contact with iron towards the claystone was as follows: magnetite, Fe-silicate, and Ca-siderite. This sequence could not be matched with the standard simulation setup (typical parameters for mineral dissolution and precipitation kinetics data, evolution of diffusion coefficient using Archie's law). The sensitivity analysis performed on the kinetic parameters (quartz dissolution, precipitation) was unsuccessful in generating the correct sequence. It was particularly difficult to maintain magnetite at the iron surface as the corrosion rate decreased with time. The observed mineral paragenesis could only be reproduced by using a cell-by-cell corrosion process and by attributing very slow diffusional transport properties to pre-corroded cells (i.e. by using a large value for the cementation factor). This layer “isolated” the iron surface from the claystone, allowing specific chemical conditions to develop at this location favoring magnetite precipitation.

The RTM simulations of iron corrosion will evolve in the short future to integrate more phenomenological corrosion models that calculate the evolution of the corrosion rate as a function of the geochemical conditions. Such electrochemical corrosion models have already been developed (Bataillon et al. 2010; King et al. 2014) and have to be coupled to RTM codes. This approach will be useful to simulate and interpret the existing laboratory experiments and also in situ results obtained recently (Necib et al. 2016; Schlegel et al. 2016, 2018).

A brief materials mineralogy overview and associated chemical gradient at the materials interface

Before dealing with the interaction between clay and concrete, it might be interesting to recall briefly the mineralogy of these two materials. In a deep geological disposal, clay materials can be the clay-rock itself where clay minerals are the main constituents of these rocks, but also the bentonites used for the construction of engineered barriers (Sellin and Leupin 2013). The mineralogy of clay-rocks is complex, meaning that its clay fraction not only contains pure clay mineral end-members, such as kaolinite, smectite and illite, but also mixed layer minerals (Claret et al. 2004). In addition, carbonate minerals with a range of chemical compositions and structures are present (see for example (Lerouge et al. 2013), and pyrite, along with organic matter and quartz (Gaucher et al. 2004a; Jenni et al. 2014; Zeelmaekers et al. 2015). Proportions of total phyllosilicates, carbonates, quartz, pyrite and organic matter can be found in the literature for various clay formations (e.g., Opalinus Clay, Boom Clay, Callovian-Oxfordian formation, and Boda Clay, in Altmann et al. (2012)) and can be very different from bentonites (e.g., Ufer et al. 2008; Savage and Cloet 2018).

The mineralogical composition of cementitious materials is also complex. Concrete is a composite material made of a porous matrix (the hydrated binder) filled with water, into which are embedded filler materials such as quartz and calcite, which act as a granular skeleton. As the hydration reaction proceeds due to the cement–water interaction, the anhydrous phases, calcium silicates (C3S or C2S) and calcium aluminates (C3A and C4AF) are converted into hydrates such as portlandite, C-S-H (calcium silicate hydrate), ettringite, monosulfate or monocarbonate. This decreases the bulk porosity, since the molar volume of the hydrates is much larger than that of the anhydrous phases (Van Damme et al. 2013; Gaboreau et al. 2017; Claret et al. 2018b). Concrete has a long history that began in the pre-Roman age, and its formulation, which was greatly improved at the beginning of nineteenth century with the invention of Portland cement, has recently become more and more sophisticated. Within the framework of nuclear waste disposal, such sophistication has been introduced with the development of low alkaline concrete (Codina et al. 2008). The idea behind this development was to improve concrete compatibility with the repository environment, yet it remains a high-strength concrete. On one hand, low alkaline concrete (often called “low pH” cement) has a lower alkali content than that of ordinary Portland-based cement material, which may reduce the pH gradient at the interface and therefore the changes in clay in contact with the concrete. On the other hand, it has a low-heat hydration temperature, which minimizes the microcracking that can have negative consequences on cement's long-term durability. One must also keep in mind that whatever the formulation, the hydration of cement-based material is kinetically driven, and thus the composition of its pore water evolves with its curing time.

While the measured pH in concrete is in the range 10.4 to 13.6 depending on time and formulation (Lothenbach and Wieland 2006; Luke and Lachowski 2008; García Calvo et al. 2010; Lothenbach 2010; Bach et al. 2012; Lothenbach et al. 2012, 2014), the pH of pore water in clay materials is generally in the range of 7 to 9 (Bradbury and Baeyens 2003, 2009; Wersin 2003; Bildstein and Claret 2015; Lerouge et al. 2018). Therefore, even for low-pH concrete material, the pH difference at the interface may be as much as two pH units. When the pH is above 13, the alkali (Na+, K+) concentrations in pore water are higher in cement material than in clay-rocks. The reverse is true for K+ when the cement material is made of low-pH cement. Still, in relation to the pH values, partial CO2 pressures are far lower in cement materials than in the pore water of clay-rocks, where partial pressures of CO2 that are higher than the atmospheric pressure (~10−3 to ~10−2 bars) are reported (Gailhanou et al. 2009; Lassin et al. 2016). Redox conditions are reducing in both claystone (Brendler et al. 2003) and in reinforced steel cement materials (Aréna et al. 2018). Nevertheless, in the latter, a strong Eh gradient exists from the steel surface outward to the bulk concrete. As described just above, a steep geochemical gradient exists at the clay–concrete interface. In the 90's, early calculations based on mass balance assumptions only, and reported in Gaucher and Blanc (2006) and Savage et al. (2007), led to the conclusion that 0.2–1 m3 of bentonite are needed to buffer the chemical perturbation created by 1 m3 of concrete. If true, this conclusion would have been problematic for the storage concepts that rely on the properties of unaltered clay materials, and this explains why so much effort has been put in reactive transport modeling studies of both experiments that mimic clay–concrete interfaces and the long-term evolution of clay–concrete interfaces.

Reactive transport modeling of laboratory and in situ scale experiments

Before their use for large-scale simulation and long-term evolution, reactive transport models and computer codes have to be tested against experimental data in order to test and improve their robustness. Experimental data can be obtained in the laboratory with experimental apparatus that mimic storage situations (e.g., column experiment), in situ to better reproduce field conditions (this is possible thanks to many underground research laboratories that are in operation) or gathering data from natural analogs for which hydrogeology and geochemical conditions reproduced phenomena that are expected in the repository. In Table 2, reactive transport studies that are compared against experimental data dealing with clay–concrete interaction are reported. The focus of our chapter is claystone, which is why experimental studies that deal with high-pH plume but in granite (e.g., (Soler 2003; Pfingsten et al. 2006; Soler and Mader 2007)) have not been reported here. In addition, when some geochemical modeling is made without coupling chemistry and transport (e.g., Lalan et al. 2016) these references are also not mentioned as the focus is on reactive transport. Among the thirteen references discussed in Table 2, three deal with column laboratory experiments (two focusing on the same experiment), four are in situ experiments and four focus on the same natural cement analog. The latter is the Maqarin natural analog in the north of Jordan (Khoury et al. 1985; Alexander et al. 1992; Khoury et al. 1992; Chaou et al. 2017). In this area, hyperalkaline groundwater compositions are the product of low-temperature leaching of an assemblage of natural cement minerals produced as a result of high-temperature/low-pressure metamorphism of marls (i.e. clay biomicrites) and limestones. On the hydraulic downstream of the cement zone, alkaline groundwaters circulated through fractures within the biomicrite clay. On the edge of the fractures, calcite, kaolinite, silica, low amounts of illite, albite and organic matter dissolved. Within the fractures, different opening/clogging stages may have occurred, leading to a complex mineralogical pathway. Two of the modeled cement claystone interfaces were sampled during the Cement–Opalinus Clay Interaction (CI) Experiment at the Mont Terri rock laboratory (Jenni et al. 2017). This is a long-term passive diffusion-reaction experiment between contrasting materials of relevance to engineered barrier systems/near-field for deep disposal of radioactive waste in claystone (Opalinus Clay). The sampled interfaces have been (and are still) extensively characterized from the mineralogical and petrophysical point of view (Jenni et al. 2014, 2017; Dauzeres et al. 2016; Lerouge et al. 2017). The two other in situ interfaces were sampled in the Tournemire underground research laboratory. Again, these interfaces have been extensively investigated (Tinseau et al. 2006; Gaboreau et al. 2011; Bartier et al. 2013). In one of the column laboratory experiments, bentonite was exposed on one face to a solution that mimic cementitious pore water, while the opposite face was a homo-ionization solution. In the second, concrete and bentonite plugs were juxtaposed. For the column laboratory experiments, heat was applied either in isothermal conditions (60 °C or 90 °C) or by applying a thermal gradient. These are the only studies in where temperature plays a role. Except for (i) the hydration concrete-bentonite column test, where non-saturated conditions and thermal gradients imply the use of a THMC coupled model, (ii) the natural cement analog and one in situ experiment where fractures play a role, a diffusive regime was the main dominant transport process in all other experiments. The column experiments and in situ experiments were sized in the 10 cm range, whereas for the Maqarin natural analog site a ~100 m fracture was considered in simulation, but matrix diffusion perpendicular to the fracture was in the 10 cm range. In most of the publications described in Table 2, the thermodynamic data used come from existing databases that are either merged together (i.e. a clay-oriented database with a cement-oriented database) or completed with some missing data (e.g. zeolite, M-S-H). Depending on the reactive transport code used, the C-S-H representation is made either by considering a discrete calcium to silica ratio (e.g. 1.6, 1.2 and 0.8), or by considering a jennite-tobermorite solid solution. The latter approach was used for the studies that used GEM (Wagner et al. 2012; Kulik et al. 2013) for the reactive part. The main findings and conclusions of the various studies are given in Table 2; therefore, they will not be detailed in the text. However, general conclusions can be drawn. Although some discrepancies between models and experiments could appear because new or completed datasets become available after modeling (Chaou et al. 2017), reactive transport modeling shows a great capability for reproducing the experiments (e.g., mineralogical transformation pathway and clogging processes). It is also a very useful tool for performing sensitivity analyses of input parameters and testing hypotheses. One of the most critical parameters that has been described in reported studies deals with kinetics and reactive surface areas that can play a large role in sequential minerals' appearance or disappearance, as well as the localization of porosity reduction.

Reactive transport modeling of the long-term evolution of clay–concrete interfaces

As previously stated, reactive transport modeling is a powerful tool for describing and reproducing phenomena that occur at both the time and space experimental scale. However, these scales are at least four or five orders of magnitude lower than the one we have to deal with for the long-term evolution of clay–concrete interfaces with geometries that are relevant for the repository gallery scale. In that case, reactive transport models are used to bridge the scale and make predictions about long-term evolution. Our literature review (which might not be fully exhaustive) found about twenty reactive transport studies addressing these issues during the past two decades (Table 3). In these studies, anionic exclusion in clay and multi-components diffusion were not considered. The complexity of the phenomena that were modeled lay more in considering an exhaustive mineralogy and steep chemical gradients that are not always easy to handle from a numerical point of view, even for modern reactive transport codes. The number of reactive transport codes used was also considerable: PRECIP (Savage et al. 2002), GIMRT (Lichtner 1996), HYTEC (van der Lee et al. 2003), PHREEQC (Parkhurst and Appelo 1999), TOUGHREACT (Xu et al. 2004; Burnol et al. 2006), ALLIANCES (PHREEDC/MT3D, (Montarnal et al. 2007), CORE2D V4 (Samper et al. 2003), CRUNCH (Bildstein and Claret 2015), OpenGeoSys-GEM (Kolditz et al. 2012), ORCHESTRA (Meeussen 2003), MIN3P-THCm (Mayer et al. 2002). A detailed description of most of the above-mentioned codes and their history can be found in Steefel et al. (2015b). Many databases have also been used: EQ3/6 (Wolery 1983), CEMDATA2007 (Lothenbach and Wieland 2006; Lothenbach and Winnefeld 2006), Nagra/PSI database (Hummel et al. 2002), THERMOCHIMIE (Blanc et al. 2015a,b), THERMODDEM (Bartier et al. 2013). The geometries (1D or 2D Cartesian, 1D radial, 1D or 2D axisymmetric and 2D with cylindrical coordinates), the mesh size (usually in the cm range), the transport parameters, the way mineral dissolution or precipitation was accounted for (either at local equilibrium and/or in kinetics) have been recently reviewed (Claret et al. 2018b). Otherwise, in the case of the mass balance calculations reported earlier in this chapter, all reported studies are consistent with very limited spatial extension, in terms of mineral precipitation and dissolution, of cement perturbation in clay materials and vice versa. Not all the reported studies consider exactly the same starting mineralogical assemblages, the same secondary potential mineral phases, the same transport properties, or the same geometries (see (Table 3). This is inherent to the starting conceptualization, but also to the quality and completeness of the databases that have been improved with time, as well as the code capability. For example, in earlier calculations, cement material was not represented as a porous medium but rather as an alkaline plume imposed as a boundary condition on one side of the clay domain. Therefore, in order to evaluate the impact of different modeling assumptions among available studies, Marty et al. (2014) carried out calculations with a consistent set of data and input parameters arranged with increasing order of complexity (e.g., the considered geometry, the representation of porous media, the choice of secondary minerals phases). This standardized approach allowed for a proper comparison of numerical results and showed that modeled reaction pathways were mostly independent of the modeling assumptions for simulations carried out in the presence of water-saturated conditions. To compare the results of the various simulations after 100,000 years, concrete and clay degradation indicators were selected. In the concrete zone, the simulations were compared by examining both the total portlandite-dissolution (accounting for both hydrolysis and carbonation) and (ii) the ettringite-precipitation (volumetric variations above twice the initial volume) front. The dolomite-dissolution front, the smectite-dissolution front (volumetric variations below half of the initial volume) and the pH plume in the clay barrier (pH > 9) were chosen to analyze clay-rock alterations. Dolomite is one of the most destabilized minerals in the claystone (Callovo–Oxfordian formation in this study) due to secondary phase formation (e.g., saponite) incorporating magnesium in their structural formulas, while montmorillonite (bearing hydroxyl groups) contributes to the buffer capacity. The criterion chosen for pH corresponds to a limit above which the dissolution rates of many aluminum silicate phases increase significantly. Regardless of the complexity of the simulation, after 100,000 years of interaction, not one of the above criteria indicates a perturbation higher than one meter. For the cases that account for full complexity, the perturbation extension is even 5 times lower. Still, in saturated conditions and besides the simulation conceptualization, the accuracy and numerical stability of the reactive transport codes have been successfully benchmarked (Blanc et al. 2015b). While simulations conducted in saturated conditions are numerous, simulations accounting for non-saturated and nonisothermal conditions are rare, making their results more difficult to evaluate. In these simulations, the main uncertainty lies in the period necessary to reach complete saturation inside the concrete. Considering an initial saturation of 30 %, this period is estimated to be approximately 2,000 years (Burnol et al. 2006), but it may depend on the grid size resolution and on the power exponent of the Millington relationship that describes dependency of pore aqueous and gas diffusion coefficients as a function of saturation and porosity (Trotignon et al. 2011). From the mineralogical point of view, the transformation pathways in concrete were found to be similar to those simulated in saturated conditions (e.g., portlandite degradation, carbonation, Ca to Si ratio decrease in C-S-H), except for sulfate/carbonate salts being deposited where drying takes place.

The glass–(iron)–clay interface

This interface is studied in countries where the spent fuel is retreated, in total or in part, by removing and recycling uranium and plutonium, e.g., in Belgium, France, Germany, Japan, Russia, UK, and the USA (Gin et al. 2013). In this case, borosilicate glasses are developed as the waste matrix for HLLW. Note that preliminary studies have also been conducted for the use of glass for ILLW (e.g., in France, UK and the USA). This interface is complex, since iron may be present if steel corrosion is still going on when water enters the canister. Once corrosion is finished, corrosion products will be present and have an influence on glass alteration. RN can be considered as traces and are not thought to have a great impact on the glass alteration process.

Numerical studies of glass alteration at the space scale of the disposal cell (50 m) and the corresponding time scale (100,000 years) are scarce (Bildstein et al. 2007, 2012). In these simulations, glass alteration follows a very simple “operational” model in which the alteration proceeds in two stages, after a time lag corresponding to the time before canister failure (700 years): a first initial phase where the alteration rate (r0) is high and a second stage where a residual, much lower, rate has been established (~r0/10,000). The feedback between glass alteration and the chemical environment is indirectly taken into account by the duration of the initial stage, which stops when silica derived from the glass has saturated the “sorption” capacity of the system close to glass (essentially the corrosion products). In these simulations, iron corrosion starts at the beginning of the simulations so that from 700 to 45,000 years, when corrosion is completed, it proceeds simultaneously with glass alteration favoring the precipitation of Fe-silicates.

The most noticeable evolution in the RTM approach to glass alteration in these two studies concerns the complexity of the geochemical and mineral system. Bildstein et al. (2007) took into account only four simple mineral phases as glass alteration products, including one zeolite. The second paper (Bildstein et al. 2012) includes 26 secondary minerals, with 14 glass alteration products, some of them determined from experimental studies (including Fe-silicates and Fe-aluminosilicates)(Figure 3). An explicit limitation of the H2 partial pressure (hydrogen is considered to form a gas phase when pH2 reaches 60 bar) was also taken into account as well as the presence of an excavation damaged zone with degraded transport properties in contact with the canister. Finally, to approach realistic conditions, the thermal gradient evolutions determined from THM calculations were included in the calculations, with temperature-dependent thermodynamic, kinetic and transport properties.

The relative simplicity of the simulations, in terms of glass behavior, can be explained by the fact that the alteration of glass in the presence of complex mineralogical environments has been only studied extensively in the last decade, allowing for specific models for glass alteration to be developed and calibrated.

Along with transport by diffusion, the dissolution of corrosion products and clay minerals and the precipitation of secondary crystallized minerals should be the main processes governing the long-term fate of glass. Dissolved cations such as Fe2+, Mg2+, Ni2+ have been shown to precipitate with silicon (if the pH is sufficiently high), which is the major element of the glass protective amorphous layer, and sustain glass dissolution (Aréna et al. 2016, 2017, 2018). The influence of Mg-rich Callovo–Oxfordian groundwater on glass dissolution has been modeled and fitted to the experimental results in Jollivet et al. (2012). The mechanism hypothesized by Rebiscoul et al. (2015) for modeling the long term alteration rate of glass in the presence of magnetite is slow magnetite dissolution followed by iron silicate precipitation. The ability of magnesium-rich minerals like hydromagnesite, dolomite and the argillite clay fraction to provide magnesium and sustain silicon consumption and glass alteration has been experimentally studied and modeled by (Debure et al. 2012, 2013, 2018). The effect of clay on glass alteration has often been described with a partition coefficient (Kd) of silicon in the clay, reflecting the complexity of the underlying dissolution/precipitation mechanisms (Godon et al. 1989; Gin et al. 2001; Pozo et al. 2007).

Current research in glass alteration modeling focuses on the description of the amorphous layer (formation, structure, composition, solubility) and on the capacity of the amorphous layer to lower the glass dissolution rate (Gin et al. 2016, 2018). Following an approach initiated by Bourcier et al. (1994), the hypothesis of a backward reaction of formation of the amorphous layer from the bulk fluid composition is made in order to build a thermodynamic description of the amorphous layer (McGrail and Chick 1984; Daux et al. 1997; Frugier et al. 2009; Rajmohan et al. 2010; Steefel et al. 2015a; Frugier and Godon 2018). The amorphous layer composition and solubility are currently described by concatenating database minerals or hypothetical phases that account for both the amorphous layer composition and the fluid composition at steady state. The GRAAL model (Frugier 2008; Frugier and Godon 2018) computes the local amount of protective layer to calculate the local thickness of the protective layer and deduce the glass alteration rate.

The remaining challenge is to couple the specific model for glass alteration (e.g., GRAAL) with large-scale RTM simulations including the materials present in the near field of the disposal cell, i.e. to ensure compatibility between the two modeling systems.

The glass–concrete–(clay) interface

This interface plays a role in deep geological repository, especially in the Belgian concept of concrete supercontainer. Once the overpack (carbon steel) and primary container (stainless steel) have failed, glass will be in direct physical contact with concrete. There are currently no RTM studies of this interface at the scale of the disposal cell, to best of the authors' knowledge. This is partly due to the challenge of understanding the complex behavior of glass in water and its interactions with materials such as concrete, which has only been studied in recent years.

The glass–concrete interface has been a growing field of interest since Belgium proposed the supercontainer concept for the disposal of high-level wastes (Kursten and Druyts 2008). A thick concrete buffer is meant to maintain high pH, to preserve the integrity of the carbon steel overpack at least during the waste's thermal phase. France is studying the option of filling the gap between the rock wall and the casing with bentonite/cement grout. The grout thickness is dimensioned to compensate for acidic fluids coming from the oxidation of iron and sulfur in the clay occurring during the facility's operational period. With high pH for long periods of time, the Belgian concept relies more on the overpack's lifetime, whereas the French concept takes more advantage of the confinement properties of the waste glass. However, both require understanding, deterministic models, and long-term quantification relative to glass alteration close to concrete.

Alteration of glass at a pH higher than 10 results in the precipitation of other secondary crystallized minerals, especially zeolites, as well as C-S-H phases. Zeolite precipitation has been proven to sustain glass corrosion and to change the composition of the protective amorphous layer at the glass surface (Ribet and Gin 2004; Ferrand et al. 2013, 2014; Mercado-Depierre et al. 2013, 2017).

The higher the pH, the faster the zeolite nucleation-growth process and the slower the induction time required for zeolite surface to be high enough to impact the glass dissolution rate (Fournier et al. 2017). Zeolite precipitation consumes Si and Al, two key constituents in the amorphous layer, and therefore sustains glass dissolution. However, zeolites also consume alkali, which in return lowers both the pH and the zeolite precipitation rate. Zeolite-seeded experiments have been a major tool for investigating those experimental conditions where a mineral's precipitation kinetics drive the glass dissolution rate (Fournier et al. 2017). The first time-dependent geochemical modeling of alteration resumption has been achieved with the GRAAL model (Frugier et al. 2017). Only a limited number of modeling studies of glass–cement interactions exist and none of them integrated the transport component. Liu et al. (2015) report simulations of glass alteration in hyperalkaline solution from their 300-day experiments at the laboratory scale, showing a diffusion-controlled glass alteration rate at 30 °C and a shift towards a surface reaction-controlled rate after 100 days at 70 °C. Baston et al. (2017) performed batch type scoping thermodynamic modeling on the interactions between two illustrative vitrified ILLW products and two cementitious backfill, predicting limited changes in mineralogy and properties of the two materials. Improving knowledge about i) amorphous layer composition and solubility at high pH and ii) zeolite precipitation kinetics is still required before applying the model to the complex chemistry and the long time scales of geological disposal.

The iron–concrete interface in anoxic conditions

This interface will be encountered in the Belgian concept with the concrete supercontainers (steel envelope and overpack) and also in all type of repository concepts where steel reinforced concrete is envisaged. Although experimental studies on iron corrosion in concrete have been performed in the context of waste disposal (e.g., L'Hostis et al. 2011; Kursten et al. 2017), no RTM studies have been found in the literature concerning this interface in anoxic conditions.

In the highly alkaline interstitial solution of Portland concrete, the corrosion of carbon steel is characterized by a low corrosion rate (passive mode) in the order of 0.1 µm/yr to 1 µm/yr (Smart et al. 2013; Chomat et al. 2017; Kursten et al. 2017). This corrosion mode is due to the formation of a thin protective layer, with corrosion products similar to those observed at the iron–clay interface: magnetite dominates, accompanied by small amounts of hematite. The stability of this layer depends highly on the solution chemistry: pH, buffering effect, carbonate and sulfate content (Chomat et al. 2017). Corrosion can therefore resume when the conditions are no longer favorable, i.e. when concrete does not buffer the pH at high values (typically upon degradation by carbonation).

Since corrosion models are at the heart of this problematic, significant advances will be achieved in the near future by coupling RTM codes with electrochemical corrosion models (Bataillon et al. 2010; Macdonald et al. 2011).

Over the last two decades, the use of RTM in simulations of the long-term behavior of materials in waste repositories has evolved to consider geochemical systems of increasing complexity, driven by the improvement of numerical codes (more functionalities), and their efficiency (solver, parallelization), as well as the accumulation of data acquired to feed thermodynamic and kinetics databases. Overall, these simulations provide a better understanding of physical and chemical change and how materials interact; they give more confidence in the prediction of the durability of materials by converging on the alteration extent that can be expected at the interfaces between different materials over long periods of time. Benchmarking RTM codes on problems related to waste repositories (e.g., Marty et al. 2015b) and continuous database development (Blanc et al. 2012, 2015b; Giffaut et al. 2014; Lothenbach et al. 2019) also improve the robustness of the predictions. This is not to say that we fully understand all the physical and chemical phenomena or that modeling is fully representative of all the processes occurring in such complex systems. These are two different aspects to be considered which call for different treatment: either acquiring more scientific knowledge through experiments, or new developments or approaches in numerical codes. Often, both issues have to be tackled together.

An example of the “lack of knowledge” issue is the question of how redox conditions are controlled in deep geological environments (e.g., see the results of redox measurements in European project RECOSY (Duro et al. 2014). Oxygen is a very reactive chemical compound, readily producing oxic conditions along its diffusive pathway during the operational stage (see dedicated section in this chapter). These oxic conditions will not remain long after the closure of the disposal cells and the onset of the corrosion phase. The re-establishment of reduced conditions in the host rock is expected within a few hundreds to thousands of years. Moreover, reduced conditions do not mean reducing conditions, i.e. the formation of reactive reductants: the reactivity of hydrogen is low, sulfate reduction is thought to be very slow and only catalyzed by microorganisms (Truche et al. 2009) and argillaceous rocks only contain small amounts of ferric compounds (e.g., in the Callovo–Oxfordian claystone). The question therefore remains about the buffering capacity of the host rock with respect to additional oxidative perturbations once the hydrogen has diffused away from the cell disposal. Such perturbations could come from the alteration of glass that will release significant amounts of ferric species. This is of course a crucial issue for the concomitant migration of RNs. It is not clear whether the host rock will contain enough reactive reductants (ferrous and/or sulfide compounds, such as pyrite) to maintain reduced conditions at this stage. Most RTM codes are ready to simulate this type of situation, since redox couples can usually be deactivated to define only the reactive ones (even H+/H2 in post-corrosion conditions).

Another recurrent issue in RTM is the value of reactive surface area that should be used for minerals in the simulation of these systems (see discussion in Li 2019, this volume). An example of assessment of the reactive surface of glass dissolving in water is given by Fournier et al. (2016), who demonstrated that the use of geometrical surface area provides a consistent approach whenever it can be easily calculated. When trying to apply this approach to porous media, and even more so to evolving porous media, the question is inevitably linked to the concepts of porosity network, accessibility to water, and physical and chemical heterogeneities (e.g., Landrot et al. 2012, Noiriel et al. 2012). Variable porosity is the result of the volume balance between mineral dissolution and precipitation (Seigneur et al. 2019, this volume); clogging phenomena are predicted at the interfaces between reactive materials in the repository because secondary minerals take up more volume than primary minerals. This is particularly true for metallic components, especially if technological gaps are rapidly filled by swelling clay upon rehydration (Wilson et al. 2015). This behavior is also predicted at the interface between concrete and clay. No consensus has been reached today about how the RT processes are affected by porosity clogging, i.e. how the macroscopic mineral reaction rate changes when porosity vanishes. This is a domain where RT simulations at the pore scale could help decipher the mechanisms controlling the dissolution and precipitation rate (Molins and Knabner 2019, this volume). A corollary aspect of this question brings us to the way we consider the evolution of transport properties in (strongly) altered zones. For instance, the diffusion properties in the corrosion layer, which controls the mineral sequence of corrosion products, had to be set to very low values to reproduce the observed paragenesis (Bildstein et al. 2016). In these situations, Kozeny–Carman and Archie's law fail to reproduce experimental observations and thermodynamic properties such as solubility are also modified if pores reach micrometric size (e.g. Emmanuel and Berkowitz 2007, Mürmann et al. 2013, Bergonzi et al. 2016, Rajyaguru et al. 2019). Although the precipitation of minerals, in general and in particular in this type of pores, should theoretically be described in RTM as a sequence of steps including nucleation, scavenging (or ripening), and growth, this approach has only been used by Savage et al. (2010) in the context of waste repositories. This type of model has been developed in the last decade and could be used more widely, providing that the necessary data are acquired (Fritz and Noguera 2009, Noguera et al. 2006, 2016). Combined with the thermodynamic approach in small pores, they may also shed new light on the processes controlling mineral reactivity during porosity clogging.

Recognizing that the combination of heat release in HLLW disposal cells and the production of hydrogen gas (also in the ILLW cells) will significantly delay the resaturation of the repository near field, more THC processes have been implemented to perform simulations with realistic temperature and water saturation changes during the lifetime of the repository. Multiphase RTM has been conducted, first looking at TH processes with a source of hydrogen (e.g., Xu et al. 2008, Zheng et al. 2008, Senger et al. 2011, Treille et al. 2012). It has progressively integrated more complex chemistry (e.g., pyrite oxidation in De Windt et al. 2014; bentonite illitization in Zheng et al. 2015), but full chemical treatment has been limited so far mainly to concrete carbonation (Trotignon et al. 2011, Thouvenot et al. 2013). The coupling between the fast diffusive transport of very reactive gas has proven to be a challenge for numerical codes. In this regard, the RTM of waste repositories could benefit from the experience return of modeling in other systems such as CO2 sequestration (Sin and Corvisier 2019, this volume). Another aspect of unsaturated porous media that is not yet widely treated in RTM relates to the dependence of chemical reactivity with regard to the water content, in particular looking at how the mineral and gas solubility and the aqueous speciation relate to capillary pressure at low water saturation (Lassin et al. 2005, 2011, 2016). This process may be enhanced by the fact that some chemical reactions, such as steel corrosion and glass alteration, consume water; this coupled effect is only starting to be fully integrated in RTM codes (Seigneur et al. 2018).

Finally, even though large-scale 3D RTM calculations started to be tractable very recently thanks to improvements in computer calculation capacity, studies looking at waste repositories in this way are scarce (Trinchero et al. 2017). The main reason remains the large scale of the problems considered: metric to kilometric spatial domain, time period of 100,000 to millions of years. This explains, at least partly, why numerical simulations in this context have not yet fully benefited from recent advancements in RTM such as the diffusion in charged porous media (Tournassat and Steefel 2019, this volume). This feature would be particularly appropriate for clay materials encountered in deep geological storage to account for different “porosities” accessible to anions and cations. The same statement is applicable for modeling at the pore scale, which could be combined with the macroscopic approach, for instance, for problems involving porosity clogging.

Open access: Article available to all readers online.