In the post-closure period of a geological disposal facility for radioactive waste, leaching of cement components is likely to give rise to an alkaline plume which will be in chemical disequilibrium with the host rock (which is clay in some concepts) and other engineered barrier system materials used in the facility, such as bentonite. An industrial analogue for cement-clay interaction can be found at Tournemire, southern France, where boreholes filled with concrete and cement remained in contact with the natural mudstone for 15–20 years. The boreholes have been overcored, extracted and mineralogical characterization has been performed. In this study, a reactive-transport model of the Tournemire system has been set up using the general-purpose modelling tool QPAC. Previous modelling work has been built upon by using the most up-to-date data and modelling techniques, and by adding both ion exchange and surface complexation processes in the mudstone. The main features observed at Tournemire were replicated by the model, including porosity variations and precipitation of carbonates, K-feldspar, ettringite and calcite. It was found that ion exchange needed to be included in order for C-S-H minerals to precipitate in the mudstone, providing a better match with the mineralogical characterization. The additional inclusion of surface complexation, however, led to limited calcite growth at the concrete-mudstone interface unlike samples taken from the Tournemire site that have a visible line of crusty carbonates along the interface.

The role of Ordinary Portland Cement (OPC)-based materials in geological disposal facilities (GDFs) for radioactive wastes is often a crucial one. Such materials are not only used for mechanical support of the GDF structure, and backfilling of cavities, but also contribute to the multi-barrier concept e.g. immobilising radionuclides within wasteforms and grouting rock fractures. Cement-based materials are therefore likely to be present within, or nearby, the majority of planned GDF systems (e.g. SKB, 2011; ANDRA, 2005; NDA, 2010).

Once the facility resaturates with natural groundwaters in the post closure period, an alkaline plume will be created due to the leaching of the cement-based materials. This plume will be in chemical disequilibrium with both the other materials used in the construction of the GDF and the natural host rock formation, with potentially deleterious effects on the capability of those materials to function in their intended manner as safety barriers (e.g. Wilson et al., 2011).

Many GDF concepts also include clay-based materials in some form, whether as a swelling buffer and/or backfill (e.g. SKB, 2011), or as a host rock (e.g. ANDRA, 2005). Chemical interactions between such materials and the alkaline plume could lead to reduction in the swelling capacity, alteration of the porosity/permeability and reduction of sorption capacity of the clay. Thus cement-clay interaction is an area of concern for the majority of waste management organisations across the world and has been, and continues to be, the focus of a large number of scientific studies.

The expected qualitative mineralogical changes that will occur when Ca-rich alkaline water (typically pH ⩾ 12) contacts clay-based materials, such as the formation of secondary C-(A)-S-H phases, are well understood (e.g. Michau, 2005; Gaucher & Blanc, 2006). However, the spatial extent of the interaction and the consequences for groundwater movement remain uncertain (NDA, 2011), and this is an area that laboratory-based experimental studies struggle to address.

Large-scale experimental studies in underground research laboratories such as the ones at Grimsel, Mont Terri (both in Switzerland) and Äspö (Sweden) go some way to providing answers to these questions but timescales are necessarily limited to a few years. One may look to natural cement-clay analogue systems that have been in situ for many tens, hundreds or thousands of years (e.g. Gaucher & Blanc, 2006; Savage, 2010), which are more relevant on both temporal and spatial scales applied to the GDF. However the complex interactions, heterogeneities and histories of natural analogue systems can make their interpretation difficult and limit their use in validating numerical models. Industrial or ‘man-made’ analogues sit somewhere between large-scale experimental studies and natural analogues. They may have been in place for decades or even centuries (relatively long compared to experimental studies, but short compared to natural analogues), and the initial and boundary conditions are often better understood than those for natural analogue systems. They therefore offer an opportunity to inform mathematical models which are in turn used to understand the evolution of the GDF over the very long time scales considered by safety cases.

The Tournemire analogue

An example of a cement-clay industrial analogue system can be found at Tournemire in southern France. The Institute for Radioprotection and Nuclear Safety (IRSN) and the National Centre for Scientific Research (CNRS) have developed a research programme based around a tunnel built at this site between 1882 and 1886. The tunnel passes through Jurassic marls and mudstones, which are ∼250 m thick and comprise mostly silicates (clays, quartz, feldspars, micas) with some carbonates (calcite, dolomite) and small amounts of sulphides and organic matter (Tinseau et al., 2006). In the 1990s, exploration boreholes were drilled into the basement of the tunnel and subsequently filled with cement and then concrete. This concrete was in contact with the mudstone for 15–20 years before it was identified as a useful industrial analogue, overcored and extracted for mineralogical characterization (Tinseau et al., 2006; Techer et al., 2012).

The samples of interest for the present study come from the middle part of the borehole (at a depth of around 1.55 m from the tunnel floor), which was water-saturated and isolated from the atmosphere once sealed. This area is below the engineered damaged zone (EDZ) of the main tunnel. The diameter of the cores extracted were approximately 20 cm, with reaction fronts clearly visible in cross sections cut from the samples (see Tinseau et al., 2006; Techer et al., 2012; and De Windt et al., 2008). Changes in the texture and mineralogy were observed at distances of approximately 1 cm from the concrete contact in the case of the undisturbed mudstone, and 1–2 cm when a fracture was present (caused by damage during the original drilling process; such fractures are typically less than 1 mm wide and extend a few cm). Beyond this the mudstone matrix appeared largely undisturbed. Detailed analysis revealed three distinct front zones within the mudstone, summarized in Fig. 1. (1) The first few millimetres around the concrete/mudstone interface were characterized by the precipitation of carbonate polymorphs (vaterite, aragonite and calcite), gypsum and ettringite, and the dissolution of the clay and quartz contents of the mudstone. There was a corresponding decrease in porosity in this section of the mudstone. (2) A second cm-scale zone is characterized by calcite precipitation and honeycombed neoformed clay-like phases with the presence of calcium aluminosilicate hydrate (C-A-S-H) and calcium silicate hydrate (C-S-H) phases with high Si content. (3) In a third cm-scale zone, more calcite precipitation is observed, with possible K-feldspar overgrowths. The first stage of mineralogical study (Tinseau et al., 2006) suggested Na-zeolites were also present in this third zone, but later studies (Techer et al., 2012) have not confirmed this observation. Within the concrete an increase in porosity was observed near the interface, with a corresponding decrease in the portlandite content. There was also an increase in the amount of C-S-H gel, calcite and ettringite present.

Modelling the Tournemire system helps to build confidence in conceptual and numerical models of cement-clay systems (which will ultimately be used to inform safety cases for GDFs for radioactive wastes) and identify the key processes involved in the alteration of the clay by a high pH plume. The remainder of this paper describes how such a model was set up, and the results it produced.

Reactive Transport Model

Governing equations

The model is conceptually based on that originally devized by De Windt et al. (2008), which demonstrated that a kinetic approach to mineral dissolution and precipitation is required to accurately represent the mineralogical changes in the mudstone. In general, a good match was obtained by De Windt et al., with the exception that muscovite was predicted instead of K-feldspar overgrowth, and zeolite growth had to be somewhat suppressed. This work aims to build upon that of De Windt et al., by using (where possible) more recent input data and the most up-to-date computational techniques.

The current model is implemented in QPAC (Quintessa, 2012), a general-purpose modelling tool with a dedicated reactive transport module. A finite-volume scheme is used to discretize the system and the software is capable of solving fully-coupled nonlinear problems.

In geochemical modelling, a set of aqueous species is adopted as a ‘basis’. One basis species exists for each element in the periodic table, e.g. SiO2(aq) is often used as the Si basis, and so on. The basis species provide building blocks from which other chemical species (‘aqueous complexes’ and minerals) can be formed. Consider a system containing NB basis aqueous species, NC complex aqueous species and NM minerals, with NIE ion exchangers (e.g. Na-montmorillonite, Ca-montmorillonite etc.) and NSC surface complexation species (e.g. S1O, S1OH, etc.). Assuming a constant water density, the transport and reaction of basis species i in a saturated porous medium is described by:  
where ρw is the density of water (kg m–3), ϕ is the porosity (-), bi is the concentration of the basis species in solution (mol kg–1), and De is the effective diffusion coefficient (m2 s–1). The summation terms represent the consumption (or release) of basis species by aqueous complexes, mineral reactions, ion exchanges and surface complexation respectively. The rates of these reactions in the bulk volume are denoted by RC, RM, RIE and RSC respectively (mol m–3 s–1), whilst νij, μik, ξil and κim represent the stoichiometry of basis species i in the reaction for complex species j, mineral k, ion exchange reaction l or surface complexation reaction m.
Similarly the concentration of the aqueous complex species, cj (mol kg–1), satisfies:  
It is assumed that local equilibrium is always maintained between the complex and basis species, thus the concentration of the complex species is determined from the mass action equation:  
where γ denotes the activity coefficient for the basis or complex species (kg mol–1) and Kj is the equilibrium constant for the reaction (-) which can be determined experimentally and is available from thermodynamic databases that are prepared for use with geochemical modelling codes. The activity model used in the simulations presented here is the Davies Equation (Davies, 1962).

Mineral kinetics

Mineral reactions are modelled using a kinetic assumption, with a standard rate equation of the form (e.g. Lasaga, 1998):  
where mk is the concentration (mol m–3) of mineral k per unit volume, Ak is the reactive surface area per unit volume (m2 m–3) of mineral k, k0k is the kinetic rate constant (mol m–2 s–1) and aH+ is the activity of H+ in solution (-). The reaction quotient, Q/K, is a measure of how near to equilibrium the solution is with respect to the mineral; Q is the ion activity product, whilst K is the equilibrium constant for the reaction. If Q/K > 1 the mineral precipitates; if Q/K < 1 the mineral dissolves. Note that this leads to an asymmetry between precipitation and dissolution (the range for the former is unbounded whilst dissolution can only occur for Q/K in the range [–1,0]), but a lack of data and experimental evidence usually results in the same form and kinetic rate being used in both cases.

Ion exchange

Ion exchange is assumed to be an instantaneous process, and is governed by a mass action equation in a similar manner to aqueous complexation:  

Here the Gaines Thomas approach (Gaines & Thomas, 1953) is adopted; thus the mass action equation involves the molalities, b (mol kg–1), of the basis species and the equivalent fractions E of the ion exchangers. The stoichiometries of these species are denoted by ξil and ξml respectively, for basis i and ion exchanger m in the ion exchange reaction for ion exchanger l. The equilibrium constant for the reaction, KlIE, is usually provided from experimental data.

The rate of change for an ion-exchange mineral is thus modified from equation 4 to  

Only one ion-exchanging mineral in any group (the “master” species) appears in the reactions of other ion exchangers, so the summation is reduced to a single term for all exchangers apart from the master.

The total number of exchange sites in a group is conserved unless the exchange minerals precipitate or dissolve. Thus the final equation needed to describe the ion exchange process is  

Surface complexation

The Bradbury and Baeyens model for surface complexation reactions on montmorillonite is used here (Bradbury & Baeyens, 1997). This does not take into account the corrective electrostatic term used by other models (see e.g. Bethke, 2008), and it has been argued by Davis et al. (1998) that the non-electrostatic approach may be the most appropriate for complex environmental applications since the surface charging behaviour of non-ideal natural mineral phases is not well known.

Like ion exchange a mass action equation governs the process of surface complexation:  

Here sn is the molality of surface species n (mol kg–1) and τnm is the stoichiometry of that species in the reaction for surface species m. The equilibrium constant for the reaction, KmSC, is determined experimentally.

The rate equation for surface species m is as follows:  
where mHOST is the concentration of the host mineral (mol m–3) and RMHOST is the associated precipitation/dissolution rate of that mineral (mol m–3 s–1).
As with ion exchange, the total number of complexation sites is conserved unless the host mineral dissolves or precipitates; thus the final equation in the suite is  


The system was discretized using a 1D radial geometry, with both the concrete and mudstone regions included explicitly. The actual system has an inner cement ring, surrounded by concrete, and then the mudstone; however for simplification both the cement and concrete were modelled as a concrete. Since the alteration is confined to a few centimetres around the concrete/mudstone interface, and the system is diffusion-limited, this approximation was deemed sufficient.

The radial geometry used in the model was centred on the borehole as shown in Fig. 2. A total of 46 cells were used radially in the concrete section, with resolution increasing towards the concrete/mudstone interface (cell widths of 2.5 mm to 0.2 mm). Similarly, 145 cells were used in the mudstone section, again with greater resolution used adjacent to the interface (cell widths ranging from 0.2 mm to 6 mm). No angular discretization was used due to the symmetry of the system.

Boundary conditions

A constant concentration boundary condition is applied to the outer mudstone boundary (at r = 50 cm), using the natural mudstone porewater. Symmetry arguments allow all other boundaries to be closed.

Input Data

Mudstone composition

The mudstone composition used by De Windt et al. (2008) was adopted for the present study and is shown in Table 1. Clays are the most abundant minerals and are composed mostly of kaolinite and ‘illite’ (pure illite and interstratified illite/smectite with more than 70% of illite), while smectite represents a minor component. Detrital muscovite and K-feldspar have also been identified. In the present model (and also as described by De Windt et al.), montmorillonite represents the smectitic part of the mudstone, and muscovite is used a proxy for illite. Calcite is the predominant carbonate while dolomite is present in smaller proportion. Framboidal pyrite, which represents less than 1 wt.% of the mudstone, was not considered in the calculations of De Windt et al., and is also excluded from the present model.

Concrete composition

The concrete filling the borehole was modelled by De Windt et al. (2008) by a mixture of about 35 wt.% Portland cement phases (C-S-H, hydrotalcite, portlandite and sulfoaluminates) and 65 wt.% calcareous aggregates (calcite).

The main difference here is that the C-S-H phase is represented by an ideal solid solution, whereas De Windt et al. used a fixed composition. The ideal solid solution is based on the model developed by Kulik & Kersten (2001), which has since been redefined by Matschei et al. (2007) and Lothenbach et al. (2008), and uses tobermorite-like (C-S-HTob) and jennite-like (C-S-HJen) gels to represent Ca/Si ratios of 0.83 and 1.67, respectively.

Other phases in the concrete remain the same as those used by De Windt et al., except for monocarbonate replacing monosulfate since preliminary calculations indicated that the former would be rapidly converted to the latter in the environmental conditions assumed. The full concrete composition used is given in Table 2.

Porewater compositions

The initial porewater compositions of the concrete and mudstone regions are based on those described by Beaucaire et al. (2008), and are given in Table 3. The mudstone porewater is of Na-Cl-(HCO3) type, with a Cl content of approximately 10 mmol l–1 and a pH close to 8. The bulk chemical analysis of the concrete indicated the presence of both potassium and sodium, interpreted as the presence of a K-Na-OH pore fluid with pH 13.4 at 15°C (the average temperature within the tunnel). The main assumption associated with these compositions is that concentrations of some key components are controlled by the equilibrium solubility of minerals/solids present in the concrete and mudstone, rather than by ion exchange or other processes.

Thermodynamic data

The thermodynamic data used in the modelling were mainly chosen from the Lawrence Livermore (LLNL) database, ‘thermo.com.V8.R6.230’, and calculated at 15°C using the polynomial interpolation included within Geochemist’s Workbench (Bethke, 2008). This database was used mainly to take advantage of its large range of data for zeolite-type minerals; however it lacks comprehensive data for cement phases. Thus the database was supplemented by cement data extracted from the compilations of Matschei et al. (2007) and Lothenbach et al. (2008). In addition, data for vaterite were taken from Busenberg & Plummer (1982). The thermodynamic data used for all minerals considered in the study are given in Table 4, and for aqueous species in Table 5.

In reactive transport simulations it is not feasible to allow every mineral in the database to precipitate; rather a list of potential secondary minerals is suggested by the user. The selection of secondary minerals was based on previous similar modelling studies, supporting calculations in Geochemist’s Workbench (Bethke, 2008), and in part on the results of De Windt et al. (2008). This list includes crystalline C-S-H minerals jennite and tobermorite, which are only permitted to precipitate in the mudstone section (C-S-H is represented by a solid solution in the concrete section). The C-S-H materials in the mudstone at Tournemire investigated by TEM show strong crystalline morphology (e.g. figure 2e in Techer et al., 2012) as compared with the clearly amorphous C-S-H in the concrete samples (e.g. figure 4a in Techer et al., 2012). Previous incarnations of the present model (Watson et al., 2011, 2012) used a solid solution throughout both concrete and mudstone, but this representation failed to predict precipitation of C-S-H phases in the mudstone.

The molar volumes and weights of the minerals are read directly from the LLNL database, but as the cement minerals are not present, values for these minerals are given in Table 6.

Kinetic data

The reactive surface areas of the various minerals were modified from those presented in De Windt et al. (2008) by decreasing values for kaolinite, illite and montmorillonite to reflect ‘edge’ rather than ‘total’ surface areas (e.g. Kline & Fogler, 1981), as shown in Table 7. The reactive surface areas of the secondary minerals were calculated assuming a crystal size of 100 μm.

Kinetic rate data were selected from the literature and are shown in Table 8. Activation energy data are not available for all minerals, so data quoted here are for a temperature of 25°C. Cement minerals (including crystalline tobermorite and jennite) were assigned a uniform rate of 10–12 mol m–2 s–1 (e.g. Baur et al., 2004).

Ion exchange and surface complexation

A further improvement on the original De Windt et al. (2008) model considered here is the inclusion of ion exchange and clay surface protonation-deprotonation reactions. Both of these processes can have a large impact on the pore water composition and hence the relative saturation states of minerals, and thus variant cases were considered that included neither process, just ion exchange, and both processes together.

Data for the ion exchange reactions is given in Table 9. The initial fraction of each ion exchange species is determined internally by the software, assuming equilibrium with the porewater. Data for the clay edge site reactions are shown in Table 10. Two types of surface sites are considered, type 1 (S1) and type 2 (S2). These have different equilibrium constants but are assumed to have the same capacity, 4.0 × 10–2 mol kg–1. The initial fraction of each type of surface site is also determined by imposing equilibrium with the porewater initially in situ.

Transport properties

Within the matrix of the mudstone, porewater contents are very low (3–5% by weight) and solute transport is substantially driven by diffusion. The effective diffusion coefficient of all aqueous species in the mudstone was based on a value for tritiated water (HTO) parallel to the argillite bedding, selected from previous measurements as 2 × 10–11 m2 s–1 for a porosity of 9.5% (De Windt et al., 2008).

Assuming non-porous aggregates and typical OPC parameters, the effective diffusion coefficient of all aqueous species in the concrete was set to 3 × 10–12 m2 s–1 for a porosity of 13%.

The effects of changing porosity on diffusion are incorporated using Archie’s Law (e.g. Appelo & Postma, 2007):  
where De is the effective diffusion coefficient (m2 s–1), Dp is the pore diffusion coefficient (m2 s–1), ϕ is the porosity (-) and m is a constant, in this case taken to be 2. Note that De Windt et al. (2008) used a constant effective diffusion coefficient, which does not account for the coupling of mineral dissolution/precipitation with transport of the aqueous solutes.

Porosity considerations

A small ‘residual’ porosity (0.001%) is enforcedly retained in the model, which allows the simulation to continue beyond pore blocking events (but greatly reduces the effective diffusion coefficient in blocked cells, inhibiting transport). Whilst this is introduced primarily as a means to progress the numerical solution, it is not unreasonable to assume that minerals do not precipitate perfectly uniformly and that a small amount of pore space is still available in regions where precipitation is not possible. The question of whether this pore space is connected for transport is a difficult one, but here it is assumed that transport paths do still exist, although the use of Archie’s Law (as described above) impedes movement of aqueous species at low porosity.

The discretization imposed in a numerical model also plays a role here – in reality, there will always be a face of a blocked pathway that is exposed to porewater (and therefore reaction, which may act to dissolve the blockage at a later stage in the evolution of the system). In the numerical model presented here, porewater cannot react with minerals in adjacent cells where the porosity is zero – there must be some transport between the two cells to initiate the reaction. Thus the residual porosity employed here is also a means of representing reaction at blockage end-faces.

On the time scales of interest here (∼15 years) this very slow diffusion through “blocked” compartments is essentially equivalent to no transport. When considering models of systems where much larger timescales are involved (e.g. geological disposal facilities for radioactive waste) this modelling assumption may have to be revisited.

The other consideration that must be taken into account in models of this type is the relationship between the discretization employed and the time at which pore space clogging (and hence restriction of transport) occurs (e.g. Marty et al., 2009). Thus it is important to avoid large cells in areas where clogging is likely, such as around the concrete-mudstone interface in the present case. Similarly it is important to avoid cells that are too small in regions where clogging is likely. Total clogging of a cell, regardless of its size, will lead to an almost complete reduction in diffusive transport across the clogged region. Therefore cells should not be chosen to be smaller than a thickness of material that is genuinely believed to lead to a halt in transport in the real system (i.e. a thickness that would lead to a chemically and mechanically stable barrier that could realistically resist flows/transport). Choosing cells that are smaller than this size will lead both to an underestimation of the time required for clogging and possibly an underestimation of the width of the alteration zone (dependent on the relative magnitudes of diffusive transport and mineral reactions).

The approach taken here was to use finer discretization around the media interface, and to gradually increase resolution in the model until a convergent solution was obtained in terms of depth of alteration and approximate time of clogging events (to the residual porosity). In the final model, the minimum cell width was 0.2 mm.


Base model – excluding ion exchange and surface complexation

A volume fraction plot of a 1.5 cm section around the concrete/mudstone interface after 15 years is shown in Fig. 3. The initial mineralogy can be largely inferred from the concrete and mudstone sections lying furthest from the interface. As expected, the mineralogical alteration is concentrated around the interface where the high-pH water from the concrete mixes with the lower pH natural groundwater. The dissolution of portlandite releases Ca2+, leading to precipitation of calcite on both sides of the interface, particularly on the concrete side. Ettringite also precipitates near the interface on the concrete side, as Al3+ is released from dissolving aluminosilicates on the mudstone side and SO42 levels increase due to the diffusion of the natural mudstone porewater into the concrete.

The Ca/Si ratio in the concrete is decreased as the C-S-H gel solid solution shifts from a C-S-HJen to C-S-HTob end-member composition, caused by the inwards diffusion of SiO2(aq) from the dissolving mudstone minerals. The influx of silica and release of Ca2+ from the dissolving portlandite leads to an increase in the total amount of C-S-H within the concrete close to the interface, and also to the precipitation of tobermorite in the first mudstone cell (immediately adjacent to the concrete).

The profile of the porosity around the interface at a number of intervals between 0 and 15 years is shown in Fig. 4. In the bulk of the concrete there is a general increase in porosity, caused by the dissolution of portlandite (consistent with the observations summarized in Fig. 1), but a reduction in porosity occurs nearer the interface due to the precipitation of calcite, ettringite and C-S-H phases (though blockage is not to the minimum allowable porosity). In the mudstone, porosity reduction only occurs immediately adjacent to the concrete due to calcite and tobermorite precipitation; in the remainder of the mudstone there is a small increase in porosity, tailing off deeper into the mudstone, caused by dissolution of the primary minerals.

This model replicates many of the features observed by the characterization studies (summarized in Fig. 1). In the concrete, portlandite dissolves; there is an increase in C-S-H, calcite and ettringite towards the interface and porosity increases (away from the immediate vicinity of the interface). Around the interface itself the observation of a thin crust of carbonates was represented in the model by precipitation of calcite on both sides of the interface. The precipitation of aragonite and vaterite in addition to calcite in this region is a classic Ostwald step sequence but the metastable states cannot be captured by the present model, although it is possible to include this mechanism (e.g. Savage et al., 2010). In the mudstone, no ettringite precipitation in the first few mm was predicted by the model, but quartz dissolved over the first 1 cm or so, as did the clay minerals muscovite (illite) and montmorillonite. Calcite also precipitated over a range of about a centimetre, but C-S-H (in the form of 14Å-tobermorite) only appeared immediately adjacent to the interface and not through a centimetre-scale section of mudstone. C-A-S-H was not included in the model due to a current lack of thermodynamic data. Finally, there were small amounts of K-feldspar overgrowth up to ∼1.5 cm from the interface.

Zeolites were not included as possible secondary minerals in this case but previous simulations (Watson et al., 2011) indicated that they precipitate readily; this could represent an analogue for the missing C-A-S-H phases (in the absence of thermodynamic data for C-A-S-H solids, low Al/Si zeolites such as stilbite and mordenite are the nearest analogue phase).

Variant 1 – including ion exchange

For the base case all the exchange sites on the montmorillonite in the mudstone were assumed to be occupied by Ca ions. In this variant, ion exchange is included and the initial fraction of sites occupied by Ca, Na and Mg ions is determined internally by the software assuming equilibrium with the initial porewater. This results in about three-fifths of the sites being filled with Na ions, with one-fifth each of Ca and Mg ions (noting that one Ca2+ or Mg2+ ion occupies two sites). During the simulation the influx of calcium ions from the concrete causes Mg ions to be displaced from the exchange sites; there is also an increase in the number of sites occupied by Na ions a few mm from the interface due to release of sodium by dissolving clays.

Figure 5 indicates that the inclusion of ion exchange in the model reduces the spatial extent of alteration around the interface; calcite still precipitates but in smaller quantities, on both sides of the interface. This can be attributed to the preferential take-up of Ca ions by the exchange sites. The second effect is that C-S-H (in the form of 14Å-tobermorite, with small amounts of jennite) appears much deeper into the mudstone, appearing in significant quantities up to ∼0.5 cm from the interface. This is probably a later effect, which occurs when the Ca-bearing montmorillonite dissolves releasing Ca2+ ions that can be taken up by C-S-H phases.

Small amounts of ettringite also appear on the mudstone side, giving better agreement with the characterization study results.

The porosity profile (Fig. 6) now shows a decrease in the porosity of the mudstone across the whole range shown, consistent with the findings of the second characterization study (Techer et al., 2012). This porosity reduction in fact now extends to approximately 1.5–2 cm from the interface, consistent with the growth of the C-S-H phases.

The inclusion of ion exchange in the model also reduces the extent of the pH plume in the mudstone, as shown in Fig. 7. In the case with no ion exchange, pH > 9 extends to a depth of ∼20 cm after 15 years, whereas when ion exchange is included pH > 9 extends to less than 5 cm. Again, this can probably be attributed to the growth of C-S-H phases in the mudstone which releases H+ ions and lowers pH.

Variant 2 – including ion exchange and surface complexation

When protonation-deprotonation reactions on the montmorillonite surface is also included in the model, calcite growth at the interface is suppressed further (Fig. 8). This deviates from the observation of the “thin white crust” clearly visible on the Tournemire samples, indicating that either surface complexation is not an important process at Tournemire, or the combined effect of ion exchange and surface complexation has been overestimated in the model.

Ettringite precipitation acts to fill up the pore space on the concrete side of the interface, and this combined with the precipitation of 14Å-tobermorite in the first 0.5 cm or so of mudstone means that the porosity profile is rather similar to that of Variant 1 (Fig. 9). The pH profile is also little affected by the inclusion of this process, as shown in Fig. 7.

In complex systems such as this, small changes in the aqueous chemistry can send the solution in a different direction. Clearly the small change in pH around the interface caused by the inclusion of surface complexation processes leads to an environment which is more conducive to ettringite growth than calcite growth, compared to the other cases considered here. It is likely that the slightly higher pH in the first few centimetres of mudstone causes more dissolution of the clay minerals, releasing more Al ions which are consumed by the ettringite.


The Tournemire industrial analogue provides a good opportunity to study a large-scale, in situ, long-term cement/clay system. Acceptable simulation of this system helps to build confidence in the predictive capabilities of reactive transport models, which will ultimately be used to inform performance assessment models that are used to develop safety cases for the geological disposal of radioactive waste.

The modelling work of De Windt et al. (2008) has been built upon in this study by using the most up-to-date data and modelling techniques, and by adding both ion exchange and surface complexation processes in the mudstone. It was found that ion exchange plays an important role in the model, reducing the extent of the alkaline plume and providing favourable conditions for the precipitation of C-S-H in the first few centimetres of mudstone (represented in the model by the crystalline mineral 14Å-tobermorite). The inclusion of surface complexation in the model (in addition to ion exchange) overly suppressed the growth of calcite around the interface, contrary to the visible crust of carbonates seen in the samples taken from the site. This may indicate that the processes are over-represented in the model.

All the main features reported by the two main mineral characterization studies (Tinseau et al., 2006; Techer et al., 2012) have been replicated by the model (when ion exchange on montmorillonite surface sites is included). This includes the dissolution of portlandite, increase of C-S-H, calcite and ettringite in the concrete; the precipitation of calcite around the interface itself; dissolution of the clay minerals and quartz in the first few millimetres of the mudstone (with small amounts of ettringite precipitation); precipitation of calcite and C-S-H phases up to ∼1 cm into the mudstone; some K-feldspar overgrowth extending to ∼2 cm; and a general decrease in the porosity of the mudstone over a scale of a few centimetres.

The model was unable to capture the presence of vaterite and aragonite as well as calcite at the interface. The inclusion of Ostwald step sequencing may address this, allowing the presence of such metastable phases to be modelled. A lack of thermodynamic data for C-A-S-H phases also led to their exclusion from the model, despite the evidence from the mineralogical studies of their presence. It is hoped that ongoing work, partly funded by the Long-Term Cement Studies project (of which this work is also a part), will produce some useful data for C-A-S-H in the near future.

This work was carried out as part of the Long-Term Cement Studies (LCS) project and was funded by the Nuclear Decommissioning Authority Radioactive Waste Management Directorate (NDA RWMD) in the UK.

The authors would also like to thank the reviewers of this paper for their helpful comments, and colleagues at IRSN, in particular François Marsal, for their valued contributions.