Modeling cesium migration through Opalinus clay: a benchmark for single- and multi-species sorption-diffusion models

Sophisticated modeling of the migration of sorbing radionuclides in compacted claystones is needed for supporting the safety analysis of deep geological repositories for radioactive waste, which requires robust modeling tools/codes. Here, a benchmark related to a long term laboratory scale diffusion experiment of cesium, a moderately sorbing radionuclide, through Opalinus clay is presented. The benchmark was performed with the following codes: CORE2DV5, Flotran, COMSOL Multiphysics, OpenGeoSys-GEM, MCOTAC and PHREEQC v.3. The migration setup was solved with two different conceptual models, i) a single-species model by using a look-up table for a cesium sorption isotherm and ii) a multi-species diffusion model including a complex mechanistic cesium sorption model. The calculations were performed for three different cesium boundary concentrations (10−3, 10−5, 10−7 mol / L) to investigate the models/codes capabilities taking into account the nonlinear sorption behavior of cesium. Generally, good agreement for both single- and multi-species benchmark concepts could be achieved, however, some discrepancies have been identified, especially near the boundaries, where code specific spatial (and time) discretization had to be improved to achieve better agreement at the expense of longer computation times. In addition, the benchmark exercise yielded useful information on code performance, setup options, input and output data management, and post processing options. Finally, the comparison of single-species and multi-species model concepts showed that the single-species approach yielded generally earlier breakthrough, because this approach accounts neither for cation exchange of Cs+ with K+ and Na+, nor K+ and Na+ diffusion in the pore water.


Introduction
Clay formations are being considered as potential host rocks for the deep geological disposal of radioactive waste. The reasons for such decision are among others: a) their very low hydraulic conductivity and b) the large sorption capability for many radionuclides, which make such formations as an ideal containment for isolating radioactive waste for very long time periods. The performance assessment of such host rocks as geological barriers requires the understanding and quantification of the processes affecting the migration of radionuclides by combining laboratory experiments, field studies and numerical modeling. Due to the low hydraulic conductivity, diffusion is the main mechanism of mass transport in clays. Radionuclide diffusion and sorption properties are then key processes for the safety assessment of underground radioactive waste repositories.
Diffusion coefficients and sorption parameters for radionuclides are usually derived from small-cm-scale diffusion experiments and related modeling [1][2][3][4][5][6] or from in situ diffusion tests [7][8][9]. For example, several in-situ diffusion experiments in Opalinus clay have been performed or are still on-going at the Mont Terri Underground Research Laboratory (URL) in Switzerland to obtain diffusion and retention data for moderately and strongly sorbing tracers [10]. In particular, one moderately sorbing tracer that has been extensively studied in the last 15 years is Cs + [2,3,[11][12][13][14][15] which shows non-linear sorption behavior in clays by cation exchange [16].
On the other hand, selecting the appropriate sorption conceptual model in migration models is at least as important as estimating the best sorption parameters. In general, single tracer models by using a simple K d approach have been applied for analyzing such diffusion experiments in the past (e.g. [7][8][17][18][19][20][21][22][23]). Non-linear sorption was later included through Freundlich sorption isotherms [15,24,25]. More complex geochemical models have also been applied especially when single-species models failed [26][27][28][29][30]. This later conceptual model, called "multi-component approach" allowed the use of a more mechanistic cation-exchange model [11,[28][29][30][31][32][33] and can be considered as more appropriate from a scientific point of view, even though the application of the K d approach was enough for many experiments. Other approximations used in the multiscale approach, include dual porosity models which assume the presence of dead-end pores in the clay with extensions even larger than the whole sample [29]. Further sorption models describing cesium sorption on clay, have been proposed recently [34]. Their implementation into the multi-species transport codes used here would be similar to the sorption models used for benchmark, i.e. K d look-up table and complex ion-exchange model using three surface sites. The focus of this benchmark is not comparing all different cesium sorption/diffusion models. The strength of such a multi-component approach lies in the fact that the effect of the negatively charged clay surfaces can also be included. These models can consider a diffuse layer region, from which anions are (partly) excluded and where cations are enriched (and contribute to "surface diffusion") compared to bulk porewater [35]. Inevitably, such multi-porosity approaches introduce additional parameters such as the width of the diffuse layer region. These parameters are typically unknown, but there are attempts to determine them from independent data. However, it should be mentioned that multi-species approaches are computationally much more expensive than single-species approaches, and therefore not always applicable [36]. As described before, a large amount of diffusion experiments and their modeling has been published in the literature with different model concepts. However, because of the complexity of both the system itself and the setups, it is not always possible to reproduce the same results from one model/code to another. For example, sometimes it is not clear which experimental setup and processes were taken into account in the models or there are incompatibilities in the model parameters. Additionally, code comparison for cesium diffusion in clay has not been done so far.
The benchmark presented here considers several numerical models to simulate the migration of cesium at 3 different tracer concentrations (10 −3 , 10 −5 and 10 −7 mol /L) in 1 cm length of Opalinus clay during 10 years including: i) A single-species transport model by using CORE 2 D V5, COMSOL Multiphysics, MCOTAC and OpenGeoSys-GEM, and ii) A multi-species transport model by using CORE 2D V5, MCOTAC, Flotran and PHREEQC v. 3. In the first case, non-linear cesium sorption in Opalinus clay is simulated with a single-species transport model using a tabulated non-linear sorption isotherm measured in batch sorption experiments [16]. In the second case, the non-linear sorption behavior of cesium in Opalinus clay is simulated with a multicomponent reactive transport model with cation exchange on three surface sites of Opalinus clay as described by Bradbury and Baeyens [37]. The results are compared among the codes by using two different concepts and then between the two concepts. For avoiding incompatibility of model parameters, in this exercise, both geochemical sorption models -K d and multi-species approach with Bradbury and Baeyens cesium sorption modelare defined with the same set of experimental data so that a comparison of the results of the different codes for fixed transport and thermodynamic sorption models parameters is possible.
Finally, it should be mentioned that benchmarking of different codes, approaches and model setups is an unavoidable and important part of the scientific work to increase the confidence in complex reactive transport models and their predictions.

Benchmark setup
The numerical modeling performed to simulate the non-linear sorption of cesium through Opalinus clay is related to laboratory diffusion experiments carried out by Jakob et al. [11]. High-pressure cells were loaded with Opalinus clay samples originated from the Mont Terri Underground Research Laboratory to carry out the diffusion experiments. The diffusion experimental setup is described in detail in [2]. During the diffusion experiments a constant 134 Cs (half-life = 2.0648 y) concentration gradient has been established across the fully water saturated Opalinus clay sample, equilibrated with Opalinus clay pore water prior to the start of the cesium diffusion into/through the clay. Boundary conditions were kept constant in the inlet and outlet of the system to have "simple" boundary conditions in the numerical modeling. Some specific experimental conditions were ignored in setting up the benchmark of the cesium diffusion in Opalinus clay to facilitate the implementation of specific conditions in the codes. For example, the filter plates on both sides of the sample within the diffusion cell, slightly decreasing cesium concentration at the "high concentration" side of the sample, were not included. Additionally, the effect of the "sawtooth" cesium concentration profile at the "low concentration" side of the sample, which results from frequent fresh water exchange in the reservoir to keep "low cesium concentration" were also not taken into account. Time dependent boundary conditions were not implemented in numerical models to improve the reliability when the results obtained with the different codes are compared. Radioactive decay of 134 Cs is also not considered due to the short duration of the experiment compared with the half-life of this radioisotope.

Description of the different conceptual models
Two conceptual models -a single-species and a multi-species reactive transport models -have been used within this benchmark exercise to simulate the non-linear sorption behavior of cesium diffusing through Opalinus clay. Three different cesium boundary conditions, cesium concentrations at 10 −3 M, 10 −5 M, and 10 −7 M, are assumed to test especially the nonlinear sorption behavior of cesium on Opalinus clay and its implementation in the different codes.
The concept of multi-species reactive transport used in this work accounts for transport of all ions in solution (in the pore water) coupled with thermodynamic equilibrium calculations for chemical reactions in solution and on surfaces at the same time. Details regarding transport equations, coupling and geochemical equilibrium reactions can be found in [38][39][40]. Within this formalism a mechanistic sorption model for cesium on Opalinus clay and major ions in solution [37,41] has been considered.
Both conceptual model options are implemented in some of the codes used in this benchmark: i.e. CORE 2D V5 and MCOTAC. The equations implemented for both concepts are described below.
The multi-species transport equations describing diffusion in saturated porous media consists of a set of diffusion equations, one for each solute species i: is the coordinate in space and R i [mol m −3 s −1 ] denotes a source/sink term, which accounts for equilibrium reactions between the solution and the mineral surfaces according to: n aqueous and n sorbed denote the number of species in solution and sorbed on the mineral surfaces, respectively, ε [−] is the diffusion accessible porosity (volume of void space per total volume), ρ [kg m −3 ] is the solid density as mass of solid phase per unit volume of solid phase, S j [mol kg −1 ] is the amount of species sorbed per unit mass of solid phase and ν i, j [−] are stoichiometric coefficients of the reactions. A constant mean value for the diffusion coefficient for all ions in solution is assumed to maintain charge balance during transport calculations.
The mass action law or the assumed equilibrium chemistry in solution and between solution and minerals and mineral surfaces yields constraints that define the source/sink term for each species and represent the mechanistic understanding of surface and/or ion exchange reactions of the following form In Eq. (3) an ion c 1 with charge z c1 exchanges with an ion c 2 with charge z c2 on a surface/exchange site ≡s 1 satisfying the mass action law where [C i ] = γ C i , i = 1, 2 and K s C1C2 is the formation constant. [C i ] denotes activities which are the product of activity coefficient γ [−] evaluated with the Davies equation [42] for ionic strength corrections and concentrations C i . ≡s 1 c 1, 2 are immobile species representing a surface/exchange site occupied by an ion c 1,2 (in equivalent fractional occupancies defined as equivalents of c 1 (or c 2 ) sorbed per unit mass divided by the cation exchange capacity in equivalents per unit mass), where for immobile species an activity coefficient of 1 is assumed. Selectivity coefficients for the above Eq. (3) are defined via Eq. (4) according to the Gaines and Thomas convention given in [43]; further details for the Cs ion exchange model are given in [37,41].
The concept of single-species reactive transport for cesium in Opalinus clay can be implemented in the multi-species concept very easily by using only a single reactive transport equation, Eq. (1), for cesium and an additional equilibrium reaction between cesium in solution and cesium sorbed (Cs_sorbed) of the form of Eqs. (3) and (4) where K(Cs) is the non-linear sorption isotherm for cesium, represented by the concentrations of cesium in solution, Cs, and cesium sorbed (S). This was the implementation used in CORE 2D V5, MCOTAC and OpenGeoSys. Another implementation for a single-species reactive transport concept is used in the COMSOL Multiphysics software environment.
Here, the 1D diffusion equation for dissolved cesium taking into account non-linear sorption in terms of a tabulated isotherm according to S = f(C), where S [mol kg −1 ] is the amount of sorbed cesium on the solid phase; f(C) is a function that depends on the time-and space dependent cesium concentration in solution [mol l −1 ]. The resulting transport equation is given by [11]: or where D e [m 2 s −1 ] is the effective diffusion coefficient, which in turn is related to the pore-diffusion coefficient D p [m 2 s −1 ] by D e = ε D p ; and D p [m 2 s −1 ] and α [−] denotes the tracerconcentration-dependent rock capacity factor accounting for non-linear sorption processes, which is given by: For conservative tracers, i.e. d f/d C = K d = 0, α is equal to the porosity of the porous medium. In this representation, the non-linear sorption isotherm for cesium is represented by the derivative dS/dC of the concentrations of cesium sorbed (dS) and cesium in solution, dC. The tabulated cesium sorption isotherm, either in form of Eq. (6) or in the form of Eq.(9), used for the single-species model concept, was taken from sorption measurements carried out by Van Loon et al. [16] (see also Appendix 1).

Description of the different codes
Six different codes used in this benchmark exercise are described below. Table 1 summarizes their specific application in this benchmark.
COMSOL Multiphysics® is a general-purpose simulation software for modeling designs, devices, and processes in all fields of engineering, manufacturing, and scientific research. Version 5.3 is used here and includes an additional "Chemical Reaction Engineering" module (see https://www. comsol.com/). CORE 2D V5 is a code for transient saturated and unsaturated water flow, heat transport, and multicomponent reactive solute transport under both local chemical equilibrium and kinetic conditions in heterogeneous and anisotropic media. The flow and transport equations are solved with Galerkin finite elements and an Euler scheme for time discretization [44,45]. CORE 2D V5 is capable of simulating chemical reactions such as acid-base, redox, aqueous complexation, cation exchange, surface complexation, mineral dissolution/ precipitation and gas dissolution/exsolution. The chemical formulation is based on the ion association theory and uses an extended version of the Debye-Hückel (B-dot) equation for activity coefficients of aqueous species. CORE 2D V5 is based on the sequential iteration approach to solve chemical reactive solute transport. Iterations are repeated until prescribed convergence criteria are attained [46]. The code also takes into account the feedback effect of the porosity changes due to mineral dissolution/precipitation on flow, transport and chemical parameters [47,48]. CORE 2D V5 has been extensively verified against analytical solutions and other reactive transport codes [49]. In addition, CORE 2D V5 has been widely used to model laboratory and in situ experiments [50][51][52][53][54] to simulate the interactions of corrosion products and bentonite [55][56][57], and to evaluate the long-term geochemical evolution of repositories in granite and clay [58,59].
OpenGeoSys-GEM (OGS) is the coupling of two codes, where the fluid flow and mass transport equations are solved by OpenGeoSys version 5 based on a standard finite element formulation, and the chemical processes by the GEMS3K kernel code of GEM-Selektor V3 [60]. The coupling of these two codes is referred to as OpenGeoSys-GEM, and its capabilities are described in [61,62]. Mass transport and chemical reactions are solved in a sequential non-iterative approach (SNIA); i.e., the transport and reaction equations are solved separately in a sequential manner without the iteration between them. The GEM approach as implemented in GEMS3K consists of calculating the equilibrium state of a chemical system via minimization of its Gibbs free energy. The minimization is constrained by mass balance equations where the given total amounts of chemical elements are conserved. An additional Model concept used for benchmark single-species single-species multi-species single-species multi-species single-species multi-species multi-species charge balance equation is also imposed to enforce the electroneutrality condition of the system. The equilibrium state calculated by GEMS3K provides the mole amounts of every species in the system and the composition of all solid, liquid, or gaseous phases [63]. In addition, other chemical quantities such as species activities or saturation indices that are needed for the calculation of kinetic rates of mineral dissolution are provided.
Flotran is a reactive flow and transport code [64] capable of describing coupled thermal-hydrologic-chemical (THC) processes in heterogeneous porous media under variably saturated and non-isothermal conditions. Flotran can describe two-phase fluid flow (water and air). Multi-component reactive transport can be calculated including homogeneous aqueous speciation reactions, heterogeneous gas reactions, redox reactions, ion-exchange and sorption, and mineral precipitation and dissolution. Geochemical systems can be modeled in 1D, 2D and 3D using structured or unstructured grids as discretization scheme. Flotran can solve the reactive transport equations with different finite difference methods (explicit, implicit and operator-splitting). The Flotran version used here for modeling is based on version of 2007 including some slight modifications made at university of Bern.
MCOTAC (Modular Coupling Of Transport And Chemistry) calculates groundwater hydraulics, transport of solutes and their geochemical interaction with solids in one or two spatial dimensions using a multi-species random walk method. The coupling is done sequentially, i.e. using exchange terms between the different modules with and without iterations in between. Consequently, back-coupling between different processes is included. This allows the modeling of quite complex geochemical systems, as for example, mineral dissolution and its influence on solute transport of radionuclides: Mineral dissolution influences flow and transport properties due to changed porosity, related hydraulic conductivity and also diffusivity. Changed flow and transport conditions influence in turn dissolution processes. The result is a "fully coupled" hydro-geochemical reactive transport model (see e.g. [39][40][41] or http://www.polyql.ethz.ch/pmwiki.php?n= Main.MCOTAC). PHREEQCv.3.
[65] is a code written in C and C ++ for simulating chemical reactions and transport processes (1D) in water including reversible and irreversible reactions. The program is based on equilibrium and kinetic reactions of aqueous solutions interacting with pure minerals, gases, ideal and non-ideal binary solid solutions exchangers and surfaces. Distribution of redox elements among their valence states can be modeled as well. Rate equations are completely userspecifiable in the form of Basic statements. Several types of aqueous models: Davies, Debye Hückel, Pitzer specific-ioninteraction and SIT (Specific ion Interaction Theory) are implemented as well. Sorption and desorption can be modeled as surface complexation reactions or as (charge neutral) ion exchange reactions. Ion exchange can be modeled with the Gaines-Thomas, Gapon or Vanselow conventions. Changes on pressure and temperature, density, conductance and inverse geochemical calculations are also possible. An explicit finite difference algorithm is included for calculations of 1D advective-dispersive transport, and optionally, diffusion in stagnant zones, transport in dual porosity media or multicomponent diffusion. In the multicomponent approach, species can have individual, temperature-dependent diffusion coefficients, but ion fluxes are then modified to maintain charge balance during transport. Additionally, the diffusion coefficients can be coupled to porosity changes that may result from mineral dissolution and precipitation.

Description of the geometric and geochemical setup
The model geometry, transport parameters and boundary conditions for cesium diffusion in Opalinus clay used for the singlespecies transport model concepts/codes and the multi-species model concepts/codes are the same. However, different meshing/discretization strategies were used according to the different numerical methods implemented in each codes (see section 2.2 and Table 2) to solve the set of equations that defines the system (see section 2.1) and after a rigorous sensitivity analysis of different options and criteria: a) number of nodes/volumes, b) uniform or non-uniform grids, c) time stepping d) simulation times and e) convergence tolerance (see Appendix 3). Although the experimental setup of the cesium diffusion in Opalinus clay diffusion was one-dimensional, 2D meshes were also used, especially for the finite element transport codes CORE 2D V5 and COMSOL Multiphysics (see Table 1). The optimized mesh size or discretization varies from one code to the other from constant, uniform 2-D mesh of triangular finite elements to logarithmic increasing/decreasing meshing in direction of the cesium gradient. Therefore, discretization varies from the sub mm range up to 0.003 m in x-direction for the whole model domain. Examples are given in Fig. 1.
Related to the different discretization used by different codes, also the time stepping has been done differently by the codes (fixed or automatic time stepping). Some of the codes have an automatic time stepping scheme, however, in those codes where this was not the case, a time stepping scheme has been optimized (i.e. PHREEQC) as a compromise between the spatial resolution and the computation time for a given studied period.
For the geochemical setup of boundary and initial conditions, it is assumed to have an equilibrated clay pore water (Mont Terri PI-Water [66,67], with a "high" cesium concentration at x = 0 cm (y = 0 up to y = y max ) and an equilibrated clay pore water with a "low" cesium concentration for x > 0 up to x = L = 1 cm (y = 0 up to y = y max ) initially. "Low" means a cesium concentration level of 10 −10 M, while for "high" cesium concentration level, three different calculations are performed for 10 −3 M, 10 −5 M, and 10 −7 M to investigate the correct handling of the non-linear sorption behavior by the different codes. Hence, the initial and the boundary conditions are defined by the following equations: This type of boundary conditions were already implemented in all codes used in this benchmark, except in PHREEQC where additional implementation was needed (see Appendix 4).
Values for the transport parameters used in all the numerical (single-and multi-species) modeling are fixed as constant values, not necessarily the result of a specific experiment, and as follows: porosity ε = 0.15, pore diffusion coefficient D p = 10 −10 m 2 /s and effective diffusion coefficient D e = D p · ε = 1.5· 10 −11 m 2 /s. Transport calculationscesium diffusion into Table 2 Look-up tables for the non-linear cesium sorption isotherm used by CORE 2D V5, MCOTAC and OpenGeoSys-GEM (first two columns) and COMSOL (columns three and four) for the single-species modeling benchmark. The values are taken from similar measurements as described in [16]  Opalinus claywere performed for a time up to 10 years to achieve a reasonable extent of cesium migration into the clay sample. Comparison of the different codes has been done by comparing cesium breakthrough curves calculated at different locations in the clay sample (x = 1, 5, 7 and 9 mm).

Single-species transport model concept
Here, the non-linear sorption behavior of cesium in Opalinus clay is simulated in terms of a tabulated non-linear sorption isotherm measured in the laboratory [16] (see Appendix 1) and implemented in CORE 2D V5, OpenGeoSys, MCOTAC, and COMSOL. In a slightly different numerical setups (compare Eqs. (6-9)), the sorption can be modeled using the K d approach which assumes that there is a fast and instantaneous chemical equilibrium between dissolved and sorbed species. CORE 2D V5, MCOTAC and OpenGeoSys-GEM have been modified to implement the K d value according to the cesium concentration in solution (sorption isotherm). The derivative of the sorption isotherm has to be calculated from the measured sorption isotherm and read from a file as direct input to the differential equation to be solved within COMSOL.
Herewith, the K d values are updated at each time step according to the cesium concentration in solution. The equation for the sorption isotherm, as it had to be implemented in the different codes in form of a look-up table was different (see Table 2). It includes values for pairs of cesium data: Cs in solution versus K d , Cs sorbed (S), or dS/dC values, but all based on the measured sorption isotherm (see Appendix 1). It is important to note that the K d values are associated to the elements of the finite elements mesh, whereas concentrations are computed in a node wise manner, e.g. in CORE 2D V5, which is different e.g. for MCOTAC, where concentrations are volume element based. It should also be mentioned that sensitivity analysis with different mesh discretization (between Δx = L/25 and Δx = L/1000) has been performed with all the codes (see Appendix 3) to understand the differences between the results of the different codes, and to estimate the individual sensitivity of the codes results on spatial discretization and mesh size. We also estimated the maximum error calculated for an early Cs breakthrough approximated by the mean diffusion path during time T into an infinite model domain: , where D and R are the diffusion coefficient and retardation factor, respectively (see for example [38] and references therein), which might be induced due to node or volume based grids. For both see Appendix 3.  Ten primary species and 22 aqueous species were considered for describing the artificial initial pore water composition of the Opalinus clay. Table 3 shows the primary species and their concentrations in mol/L and the calculated aqueous composition. The equilibrium constant for the aqueous complexes are listed in Table 4.

Multi-species transport model concept
The sorption model of Bradbury and Baeyens [37] has been used to simulate the non-linear sorption of cesium in Opalinus clay. Cesium uptake is dominated by cation exchange and three different sorption sites are considered: s plan ("planar sites"), s Type-II ("higher affinity sites") and s FES ("frayed edge sites"). The Gaines-Thomas convention [43] has been used for all codes, whereas the activity coefficient for sorbed species has been set to 1. Table 5 lists cation exchange capacity (CEC) of the cation exchange sites. Cation exchange reactions and their selectivity coefficients constants for the three cation exchange sites are listed in Table 6. These three sites cation exchange sorption model was used in the multi-species transport models PHREEQC, CORE 2D V5, Flotran, MCOTAC to simulate the non-linear of cesium on Opalinus clay.

Results
As indicated in Table 1, COMSOL and OpenGeoSys-GEM were used for the single-species benchmark, whereas Flotran and PHREEQC v.3 were used for the multi-species benchmark. CORE 2D V5 and MCOTAC were used for both benchmarks (single-and multi-species approaches). In the following, first the individual benchmarkssingle-and multispecies conceptsare compared among each other. Then both concepts are compared with each other, and finally codes specific behavior as well as calculation time are mentioned. All comparisons performed in the following subsections are based on calculated Cs concentrations in solution at four different locations within the model area (x = 1, 5, 7 and 9 mm). Representation of the results are shown in a logarithmic scale (log scale) for displaying all simulated data over the wide range of time (1-3650 days) and concentration (10 −10 -10 −3 mol/L). This logarithmic scale representation also provides a greater resolution of the differences between codes at the non-steady states of Cs diffusion.

Single-species concept
The results obtained with the single-species transport model concepts with CORE 2D V5, MCOTAC, OpenGeoSys-GEM and COMSOL Multiphysics are compared by using calculated cesium concentrations in solution at several locations within the model area (x = 1, 5, 7 and 9 mm). Figure 2 shows the cesium breakthrough curves calculated with the four transport codes at different locations for the single-species transport model concept assuming the cesium concentration at the "high" concentration boundary equal to 10 −3 , 10 −5 and 10 −7 mol/L. The calculations cover a cesium concentration range in solution of four orders of magnitude to test the whole range of the non-linear sorption isotherm. Cesium breakthrough computed when cesium concentration is equal to 10 −3 mol/L at the "high" boundary are earlier than those computed when cesium concentration is equal to 10 −5 and 10 −7 mol/L at the "high" cesium concentration boundary. As expected, "high" Cs concentrations defined at the left boundary are never reached at x = 5, 7 or 9 mm due to the effect of the boundary condition at x = 10 mm where Cs concentration was fixed to 10 −10 mol/L. The cesium retardation is stronger for lower cesium concentration levels at the "high" cesium concentration boundary due to the non-linear cesium sorption behavior through Opalinus clay (strong cesium sorption for lower cesium concentration levels in solution and lower cesium concentration gradient). In fact, there are no cesium breakthrough curves calculated at 7 and 9 mm in the model domain  The results computed with the four models/codes with different "high" cesium concentrations (10 −3 , 10 −5 and 10 −7 mol/L) at the boundary show a good agreement taking into account the large cesium concentration range and its nonlinear sorption behavior with differences below 1% in plateaus of the calculated cesium breakthrough curves. There are small discrepancies depending on the position for which the cesium breakthrough curves are calculated. The largest differences are found near the model inlet and the relevance of spatial discretization decreases when moving away from the "high" concentration boundary. The cesium concentration plateaus computed with all the codes at later times (at quasi steady state conditions) are also in good agreement with differences below 1%, and the early steep breakthrough times are within a range from 2% to 6% for the different locations for all codes compared to node-and volume based grid errors of maximum 1% of mean arrival time. One specific point should be mentioned already here, when looking at the COMSOL calculations for the Cs boundary condition of 10 −5 mol/L, calculated concentration fluctuations for Cs concentrations below 10 −9 mol/L at x = 0.001 m disappeared only when using a very fine discretization (see Appendix 3). Such numerical oscillations of finite element codes are typical when introducing steep concentration gradients [69]. Other finite element based reactive transport codes, like OpenGeoSys for example, use so called "mass lumping" to avoid such concentration fluctuations (see e.g. [62,70], especially at sharp concentration fronts, where these fluctuations may reach negative concentrations, which do not fit to sequential geochemical calculations in the coupled reactive transport modeling. However, procedures like this also may induce some disagreement of the codes results (see also Appendix 3). To identify the reasons for the small discrepancies between the computed results, computed sensitivity analysis to the spatial discretization used by all the codes were carried out to select suitable spatial discretization for each code, so that the comparison between the codes was reliable. The analysis of the spatial discretization is presented in Appendix 3. Finer discretization for an individual code does not always yield a better convergence for the calculated breakthrough curves. Physically adapted grids, non-uniform grids with higher resolutions at the boundaries yield a better convergence than just a finer equidistant grid. CORE 2D V5, MCOTAC and OpenGeoSys-GEM calculations show different results depending on the spatial discretization used. For all three codes, the largest differences between the Cs breakthrough curves calculated with  Na + + X-K ⇔ K + + X-Na + 1.99·10 −1 Na + + 0.5 X 2 -Ca ⇔ 0.5 Ca 2+ + X-Na 4.62·10 −1 Na + + 0.5 X 2 -Mg ⇔ 0.5 Mg 2+ + X-Na 5.07·10 −1 Na + + X-Cs ⇔ Cs + + X-Na 2.51·10 −2 Second site (Type II) Na + + X-K ⇔ K + + X-Na 7.94·10 −3 Na + + X-Cs ⇔ Cs + + X-Na 6.31·10 −4 Third site (FES) Na + + X-K ⇔ K + + X-Na 3.98·10 −3 Na + + X-Cs ⇔ Cs + + X-Na 1.00·10 −7 several grid sizes are found near the model inlet and the relevance of spatial discretization decreases for Cs breakthrough curves calculated for larger distances to the "high" concentration boundary. According to the sensitivity analysis, finer discretization to more than 100 nodes in x-direction does not yield significant improvements. For the COMSOL single- species calculation, there was no influence of the grid size noticeable for the higher Cs concentration levels of the calculated breakthrough curves. Another reason for the small disagreement might be related to the implementation of the sorption isotherm that has been done differently in the reactive transport codes and different interpolation procedures used by the codes. The cesium sorption isotherm measured in the laboratory experiments was directly implemented in CORE 2D V5, MCOTAC and OpenGeoSys-GEM. However, the derivative of sorption isotherm had to be interpolated for COMSOL for specific cesium concentration levels as direct input to the differential equation to be solved within COMSOL, while K d values were assumed to be constant in predefined intervals of the cesium concentration in solution in CORE 2 D V5, MCOTAC and OpenGeoSys-GEM.

Multi-species concept
The multi-species transport approach has been implemented in CORE 2D V5, Flotran, MCOTAC and PHREEQC v.3. The code comparison has been performed with respect to the main "perturbation" of the system, which is the Cs concentration gradient across the Opalinus clay sample. This had to be chosen because multi-species data outputs of the different codes were very different. For example, for some codes one may define the output data specifically (e.g. specific ion(s) concentration(s) as a function of time and space) others produce an output of all ion(s) concentration(s) for each (automatically generated) time step calculated, which are not necessarily identical for all codes for comparison of the results for a given time. For instance, PHREEQC generates huge output files (~GB), which are sometimes difficult to handle with some text editors and graphics programs. Therefore we focus on comparison of calculated Cs breakthrough curves. A more detailed modeling of all the processes involved in the multi-species diffusion of Cs through Opalinus clay is shown in Appendix 2 exemplified for CORE 2D V5 calculations or in [11] for MCOTAC.
Cs breakthrough curves calculated at different distances from the inlet (x = 1, 5, 7 and 9 mm) are shown in Fig. 3 considering three different initial cesium concentrations (10 −3 , 10 −5 and 10 −7 mol/L) at the "high" concentration boundary. The cesium concentration as a function of time computed with the different codes maintain, in general, a good agreement in the trend but still some differences can be appreciated. As for the single-species concept comparison (see section 3.1), higher differences between cesium breakthrough curves occur near the "high" cesium concentration boundary (at x = 1 mm) for all three Cs boundary conditions. Sensitivity analysis of the spatial discretization (see Appendix 3) showed that the effect of spatial discretization was especially relevant near the left boundary for all the codes. After the sensitivity analysis, the optimal spatial discretization was selected to compare between codes. The largest differences were identified for calculated Cs breakthroughs at concentrations levels below 10 −8 mol/L (for all three Cs boundary concentration levels) at 1 mm. The calculated breakthroughs of all codes are within a range of less than 10% (concentration levels or breakthrough time) for Cs concentration levels larger than 10 −8 mol/L, and the Cs plateaus are calculated within a range less than 1%. These differences are probably due to their different numerical methods and coupling schemes used.

Comparison between single-species and multispecies transport models
Single-species and multi-species model concept where calculated by CORE 2D V5 and MCOTAC only, therefore the following comparison is limited to these codes. The comparison of calculated cesium breakthrough curves is related to two different "high" cesium concentration boundary conditions 10 −3 and 10 −5 mol/L and at x = 5 mm and 9 mm in the Opalinus clay. The same finite element mesh with 100 nodes in x-direction has been used for calculations with CORE 2D V5, whereas MCOTAC calculations have been done with 100 volume elements. Ignoring this small node/volume meshing difference, the results computed with the two numerical models yield a good agreement for both cesium concentrations levels at the boundary (see Fig. 4), despite the fact that the simulation of the non-linear cesium sorption in Opalinus clay is completely different in both models (single-species and multi-species models). Also, the stronger retardation for lower "high" cesium concentration at the boundary is reproduced by both concepts. But there is the communality that for both codes the single-species model yield a faster breakthrough than for the multi-species approach which is about the same for both codes: 10 days earlier breakthrough at 5 mm and 20 days earlier breakthrough at 9 mm for the multi-species model concept. An explanation for this was mentioned already in Jakob et al. [11], who argued that the single-species model concept with the look-up table does not include cation exchange Cs + with K + and Na + and the related K + and Na + diffusion in the pore water, which can only be taken into account by multi-species model concepts. This is confirmed by our benchmark too.

Codes computation times
A short glance on calculation times of the different codes used in the benchmark is given in Table 7. Although these numbers could not be directly compared one to one, because different processors, machine usage etc. might have been used, it is obvious that the single-species concepts are calculated at least 25 times faster than the multi-species concepts, and the multispecies concepts are calculated about an order of magnitude, in between 1 h and 40 h roughly, also depending on discretization.

CORE 2D V5 sensitivity analysis to the convergence tolerance
The system of chemical equations is solved with an iterative Newton-Raphson method. The unknown concentrations of the chemical components at the (s + 1) iteration, x s + 1 j , are computed for those of the previous iteration x s j for j = 1, 2, …. N, N being the number of unknowns. The iterative process stops when the maximum number of allowed iterations is reached or when where ω is a prescribed convergence tolerance. The time evolution of computed cesium concentrations are very sensitive to the convergence tolerance (% relative error), especially in the case with a cesium concentration in the inlet equal to 10 −7 mol/L. In CORE 2D V5 it is suggested to use a value for the relative convergence tolerance between 10 −3 and 10 −6 to solve the chemical system. However, when a relative tolerance value between 10 −3 and 10 −6 is used in the multispecies transport model, the cesium breakthrough curves computed with CORE 2D V5 change significantly. For this reason, several sensitivity runs were performed with CORE 2D V5 to analyze the influence of the relative convergence tolerance (ω) to solve the chemical reactions. Figure 5 shows the cesium breakthrough curves calculated with CORE 2D V5 at x = 1 and 5 mm by using different values for the relative convergence tolerance (10 −3 , 10 −5 , 10 −6 , 10 −7 and 10 −11 ) considering a cesium concentration at the "high" concentration boundary equal to 10 −7 mol/L.
With this analysis, it was decided to use a value of the relative convergence tolerance equal to 10 −7 to achieve reasonable agreement with the other codes results (see Figs. 3 and 4), although longer calculations times were required (results computed with ω = 10 −7 and ω = 10 −11 are similar). To investigate such code behavior is a direct consequence of this benchmark exercise, as it was found an initial disagreement when comparing codes with a high tolerance ω = 10 −3 . Without such a possibility of comparison other code results, specific code behavior might not have been observed.

Summary and conclusions
Six different numerical codes -CORE 2D V5, Flotran, COMSOL Multiphysics, OpenGeoSys-GEM, MCOTAC and PHREEQC v.3 -were used to simulate non-linear sorption behavior of cesium in Opalinus clay as a benchmark. Two different sorption model concepts, single-species with sorption isotherm and multi-species approach with a complex sorption model were employed. The benchmark setup is related to laboratory cesium diffusion experiments carried out by Jakob et al. [11] and considers cesium migration at a large range of cesium concentrations levels (10 −3 , 10 −5 and 10 −7 mol /L) to test the correct handling of the non-linear cesium sorption.
The four numerical codes CORE 2D V5, OpenGeoSys-GEM, MCOTAC and COMSOL have been modified to implement a tabulated cesium sorption isotherm and contributed to the single-species model concept benchmark. Code comparison is related to cesium concentrations calculated in solution at several locations within the model area. Good agreement has been obtained in all cases after an appropriate spatial discretization of the systems have been used accordantly to the different numerical methods, interpolation procedures, time discretization and implementation of the sorption isotherm used by each code. Small differences could be attributed to different discretization, different implementation of the sorption isotherm and different interpolation procedures used by the codes. Sensitivity analysis showed the importance of spatial discretization on the computed results, which confirmed that applying a coarser mesh leads to greater cesium dispersion. For the COMSOL calculation, there was no influence of the grid size noticeable for the higher concentration parts of the calculated breakthrough curves. However, at the lower  For multi-species concept benchmark using the more complex cesium sorption model with cation exchange on three sorption sites [37], four codes, CORE 2D V5, Flotran, MCOTAC and PHREEQC v.3, were used. In general, good agreement was again obtained in this cesium diffusion setup, however, specific discrepancies have been identified. When comparing CORE 2D V5 and MCOTAC, agreement was shown for both but the multi-species concept was better than for the single-species concept. Some differences between results were more obvious for other codes, especially between cesium breakdown curves near the "high" cesium concentration boundary. Sensitivity analysis to spatial discretization performed showed that the effect of spatial discretization is especially relevant at 1 mm of the "high" cesium concentration boundary. Reasons for these small discrepancies might be similar to those for the single-species model concept and include different spatial discretization (node and volume centered), time stepping, or different implementation of the boundary conditions, In addition, numerical methods used to solve the differential algebraic and partial differential equations or coupling procedures between transport and geochemical calculations were treated differently between the codes. These last points correlate with time stepping or tolerance parameter values, among others, and were used to achieve better agreement, which would have been impossible without this benchmark.
When comparing results from both benchmark conceptual models (single-and multi-species transport models), only done by CORE 2D V5 and MCOTAC calculations, good agreement has been achieved for both concepts. However, the singlespecies model yielded a faster breakthrough than for the multi-species approach for both codes. This was due to the K d look-up table in the single-species model concept, which does not include the cation exchange of Cs + with K + and Na + and the related K + and Na + diffusion in the pore water. We note that this can only be taken into account by multi-species model concepts. This benchmark confirms what has already been proposed by Jakob et al. [11]. Main differences between both approaches also affect to the computational time needed, with simulation times at least 25 times faster with the single species concept than with the one obtained with the multi-species one. Although this can be seen as a clear advantage for the application of the single species concept in long term prediction of radionuclides migration, the multi-species approach can give better information/process understanding of the geochemistry and a deeper understanding of the system to be analyzed. Comparing computation times of the codes for the multispecies benchmark case, it should be noted that they vary surprisingly within a factor of more than ten. From the experience of this benchmark study, especially according to the multispecies concept benchmark, we point out where code agreement has been achieved. After improvement of discretization, time stepping, choice of tolerance value or equivalent implementation of boundary condition, it can be concluded that results for complex geochemical system setups published in the literature may include some code specific properties, yielding different results when comparing codes results. Therefore, more benchmark studies are necessary even if they are disputed between different modelers.

Measured cesium sorption isotherm
For the modeling of this benchmark, the pore water and the mineral composition are in line with that from Opalinus clay samples from the Mont Terri test site (see Fig. 6). The cesium sorption isotherm measured by Van Loon et al. [16] in both, crushed and compacted Opalinus clay samples. The sorption isotherm is used in form of a look-up table with given values for pairs of cesium data: cesium in solution versus logarithm of Cs_sorbed divided by cesium in solution ( 10 log Rd_Cs = 10 log(S/C), see Table 2). These table entries had to be recalculated with respect to units or derivatives dS/dC to be used in the different code calculations. Figure 6 indicates about the different implementations of the look-up tables for the different codes.

Appendix 2
Detailed modeling of cesium diffusion through Opalinus clay by using the multi-species reactive transport concept A more detailed analysis of the processes happening during the cesium diffusion into Opalinus clay is exemplified for Fig. 6 Sorption isotherm for cesium measured on crushes and compacted Opalinus clay from the Mont Terri site [16] and used in the single-species transport model concept CORE 2D V5 modeling of the specific case by using a cesium concentration level of 10 −3 mol/L at the "high" concentration boundary for a distance of x = 5 mm within the Opalinus clay (middle of the sample). Figure 7 shows the time evolution of the total aqueous concentration in mol/L of Na + , Ca 2+ , Mg 2+ , K + and Cs + . The increase of cesium concentration after approximately 80 days represents the arrival of the diffusing cesium front from the "high" concentration side at 5 mm in the Opalinus clay. The concentrations of all the other cations remain about constant, at least no significant differences can be appreciated on logarithmic scale. This can clearly be explained as major cations concentration in solution are much higher than cesium and therefore, the influence of the migrating cesium front is almost negligible for the major cations over time and only slightly changes visible for potassium present iñ 10 −3 mol/L concentration.
Additionally, the time evolution of the exchange composition of the planar sites is shown in Fig. 8. One can see that the composition of the exchange sites including K, Na, Mg and Ca keep constant with time and only after approximately 100 days the composition is slightly modified due to cation exchange reaction with cesium. Figure 9 shows the time evolution of the concentrations of exchanged cations in the second cation exchange site computed. Na + is exchanged with K + during the first 100 days of the simulation at x = 5 mm. The concentration of exchanged K + increases while the concentration of exchanged Na + decreases. After approximately 100 days, Cs + exchanges with K + and Na + . Figure 10 shows the time evolution of the concentrations of exchanged cations in the third cation exchange site computed with the multispecies transport model at in the middle of the sample (x = 5 mm) assuming a "high" cesium concentration equal to 10 −3 mol/L. As in the second sorption site type, Na + is exchanged with K + during the first 100 days of the simulation at x = 5 mm in the third sorption site type. When the cesium diffusion front reaches x = 5 mm from the "high" concentration boundary, K + and Na + are exchanged with Cs + . The cesium breakthrough curves is mainly determined by the temporal and spatial potassium and sodium distribution. Other cations such as calcium or magnesium play only a minor role. Figure 11 shows the time evolution of the computed concentration of dissolved Na + and Cs + at x = 1, 5 and 9 mm assuming a cesium concentration at the "high" concentration boundary equal to 10 −3 mol/L. The time evolution of the sodium concentration is not as closely correlated with that of cesium as in the case of potassium. The concentration of dissolved Na + increases with time throughout the Opalinus clay sample and subsequently it decreases to its initial value in solution (0.24 mol/L). However, the decrease of Na + begins before the concentration of dissolved Cs + begins to increase due to the cesium diffusion from the "high" concentration boundary. Despite this, the decrease of Na + is faster when the concentration of dissolved Cs + increases. Figure 12 shows the time evolution of the computed concentration of dissolved K + and Cs + at x = 1, 5 and 9 mm assuming a cesium concentration at the "high" concentration boundary equal to 10 −3 mol/L. The time evolution of the potassium concentration correlates to that of cesium. The concentration of dissolved K + increases with time throughout the Opalinus clay sample. The increase of K + is greater at points located near the "high" concentration boundary (x = 1 mm). However, when the concentration of dissolved Cs + begins to increase due to the cesium diffusion from the "high" concentration boundary, the concentration of dissolved K + begins to decrease. Finally, the computed concentration of dissolved K + decreases to its initial value in solution (1.60·10 −3 mol/L).

Appendix 3
Sensitivity analysis to the spatial discretization used by the codes for the single-and multi-species approaches Sensitivity analysis of the spatial discretization used by all the codes were carried out for both the single-and multi-species approaches. This analysis was used to select suitable spatial discretization according to the mathematical treatment used by each code to solve the partial differential equations, so that the comparison of the results does not depend on the discretization and the error associated to discretization is minimum. Obtaining reasonable computational times was another factor that was taken into account in selecting the spatial discretizations used as a reference in the benchmark exercise. The sensitivity analyses carried out with each code are shown below.  Fig. 12 Time evolution of the concentration of dissolved K + and Cs + computed with CORE 2D V5 at x = 1, 5 and 9 mm. The cesium concentration at the "high" concentration boundary is equal to 10 −3 mol/L Fig. 13 Cs breakthrough curves calculated with CORE 2D V5 at different locations of the singlespecies model by using three different spatial discretizations. The cesium concentration at the "high" concentration boundary was assumed to be 10 −3 mol/L (top), 10 −5 mol/L (middle) and 10 −7 mol/L (bottom) Fig. 14 Cs breakthrough curves calculated with CORE 2D V5 at different locations of the multispecies model by using three different spatial discretizations. The cesium concentration at the "high" concentration boundary was assumed to be 10 −3 mol/L (top), 10 −5 mol/L (middle) and 10 −7 mol/L (bottom)

CORE 2D V5
Single-species approach Sensitivity runs were carried out by using three different mesh size: a) Standard, Δx = L/50 (50 nodes in x-direction), b) fine, Δx = L/100 (100 nodes in x-direction) and c) very fine: Δx = L/ 1000 (1000 nodes in x-direction) for the three cesium concentrations levels. Figure 13 shows the cesium breakthrough curves at different locations at the single-species transport models (x = 1, 5, 7 and 9 mm) by using three different mesh sizes assuming cesium concentration levels of 10 −3 , 10 −5 and 10 −7 mol/L at the "high" concentration boundary. In all three cases, the largest differences between model results calculated with several mesh sizes are found near the model inlet (1 mm) during the nonsteady state. For example, in the models with the "high" concentration boundary of 10 −3 mol/L, Cs concentration at x = 1 mm begins to increase from the background Cs in the Opalinus clay (10 −10 mol/L) after 0.8 days if the standard discretization is used, while Cs concentration does not increase until the third day in the model with a mesh size of Δx = L/1000. Similar errors of 70% occurs when low cesium concentration levels are calculated and 10 −5 and 10 −7 mol/L of Cs at the "high" concentration boundary is used. At 5, 7 and 9 mm from the inlet of the model, differences between the model results computed with several mesh sizes are much less significant and below 10% of error. In addition, these sensitivity runs revealed some errors committed by the fact that for CORE 2D V5 K d values are associated to the elements of the mesh, whereas concentrations are associated to nodes of the mesh. This mismatch was reduced for finer finite element meshes highlighting the relevance of fine discretization for this code and the single species concept. As observed in Fig. 13, the importance of spatial discretization decreases when moving away from the "high" concentration boundary and when reaching quasisteady state conditions. The single-species models using the very fine discretization took 2.5 times longer than those that used the standard discretization (see Table 7). However, calculation times is not critical in the single-species approach since less than 26 min with the most refined mesh) is totally acceptable for a simulation case of 10 years (see Table 7).

Multiple-species approach
Three different mesh sizes: a) standard Δx = L/50 (50 nodes in x-direction), b) fine, Δx = L/100 (100 nodes in x-direction), and c) finer: Δx = L/200 (200 nodes in x direction) have been used. Figure 14 shows the cesium breakthrough curves calculated at different locations (x = 1, 5, 7 and 9 mm) by using three different mesh sizes assuming cesium concentration levels of 10 −3 , 10 −5 and 10 −7 mol/L at the "high" concentration boundary. As in the single-species model, the effect of spatial discretization is especially relevant near the left boundary of the model. In the models with the "high" concentration boundary of 10 −3 mol/L, Cs concentration at x = 1 mm begins to increase from the background Cs in the Opalinus clay (10 −10 mol/L) after 1 day in the model using the standard discretization and after 2.5 days in the model with a mesh size equal to Δx = L/200 (finer). For the cesium boundary concentration level of 10 −7 mol/L, Cs concentration at x = 1 mm begins to increase significantly after 40 and 70 days using mesh sizes of Δx = L/50 and Δx = L/200, respectively. Cesium breakthrough curves calculated at x = 5, 7 and 9 mm with multi-species models with several mesh sizes hardly differ. Calculation times are greatly increased by reducing the mesh size in the multi-species approach (Table 7). If a spatial discretization of Δx = L/200 is used instead of Δx = L/50, calculation times could be increased more than seven times. Since calculation times with CORE 2D V5 can exceed 60 h, choosing a suitable mesh size for the multi-species models is essential to avoid excessive calculation times.

Single-species approach
Four different grids in x-direction were used for the COMSOL single-species calculation. Surprisingly, there was no influence of the grid size noticeable for the higher concentration parts of the calculated breakthrough for all 4 discretisations (25 to 500 nodes) tested (see Fig. 15). However, at the lower concentration level part of the breakthrough curves, i.e. at initial arrival, larger fluctuations -partly more than an order of magnitude-were calculated for the Cs breakthroughs for the 10 −3 and 10 −5 mol/L Cs boundary conditions, but not for the 10 −7 mol/L Cs boundary condition. The fluctuations disappeared only for the fine discretization. Very likely, these behavior might have been improved by explicitly induced time stepping within the COMSOL software package, but the convergence of all Cs breakthrough curves for Cs concentrations larger than 10 −7 mol/L (in case of 10 −3 mol/L Cs boundary concentration); 10 −8 mol/L (in case of 10 −5 mol/L Cs boundary concentration) and 10 −10 mol/L (in case of 10 −7 mol/L Cs boundary concentration) seemed to be reasonable for the single-species benchmark. However, such a behavior might be disadvantageous when using COMSOL coupled to a geochemical equilibrium reaction code.

Flotran
Multiple-species approach FLOTRAN grid investigations show a similar behavior as observed for CORE 2D V5 (see Fig. 16). Increasing the number nodes to 100 for the domain yield convergence of the breakthrough curves, also for the breakthrough curve calculated for 1 mm, which shows a similar convergence behavior as for MCOTAC.

MCOTAC
MCOTAC grid investigations show a similar behavior as observed for CORE 2D V5. Increasing the number nodes to 100 yield a convergence of the breakthrough curves, also for the breakthrough curve calculated for 1 mm, which shows the largest difference when comparing CORE 2D V5 and MCOTAC 50 nodes calculations. There was no difference observed for the single-and multi-species MCOTAC convergence behavior (Figs. 17 and 18).  17 Cs breakthrough curves calculated with MCOTAC at different locations of the single-species model by using two different spatial discretizations in x direction. The cesium concentration at the "high" concentration boundary was assumed to be 10 −3 mol/L (top), 10 −5 mol/L (middle) and 10 −7 mol/L (bottom). The straight parts of some breakthrough curves are due to output data intervals Fig. 18 Cs breakthrough curves calculated with MCOTAC at different locations of the multispecies model by using two different spatial discretizations in x direction. The cesium concentration at the "high" concentration boundary was assumed to be 10 −3 mol/L (top), 10 −5 mol/L (middle) and 10 −7 mol/L (bottom). The straight parts of some breakthrough curves are due to output data intervals Error estimation due volume and node based grids Finite element based reactive transport codes like CORE 2D V5, Flotran and OpenGeoSys-GEM have node based chemical equilibrium calculations, whereas MCOTAC has volume based geochemical reactions. When comparing calculated breakthrough curves by these codes one has to keep in mind that there might be a mismatch of the calculated breakthrough curves of half of the grid size, i.e. the offset node and volume based grids, if no interpolation to in one or the other direction is performed. In order to estimate such a node/volume grid error, we assumed a mean Cs diffusion into the clay approximated by the mean diffusion path during time D is the Cs diffusion coefficient and R is the retardation factor,  (Table 8). It is obvious that the node and volume grid based errors decrease with smaller grid size. The error increases with increasing breakthrough location and decreases with increasing Cs boundary concentration and related decreasing sorption coefficient. The maximum calculated error varies, depending on grid size, assumed Cs boundary concentration and location of calculated Cs breakthrough, between 0.18 days and 166 days (dx = 0.0002 m, high Cs concentration (10 −3 mol/L) at the boundary, low Cs retardation, at x = 0.001 m and low Cs concentration (10 −7 mol/L) at the boundary, high Cs retardation, at x = 0.009 m) and between 0.09 days and 83 days (dx = 0.0001 m, high Cs concentration (10 −3 mol/L) at the boundary, low Cs retardation, at x = 0.001 m and low Cs concentration (10 −7 mol/L) at the boundary, high Cs retardation, at x = 0.009 m).

Single-species approach
The investigation of grid size and type and using mass lumping option in OpenGeoSys is shown in Fig. 19. The calculated breakthrough curves show some dependency for the 1 mm breakthrough for the high Cs concentration boundary condition in case of uniform coarse (59 nodes) and the fine (100 nodes) grid and the coarse (59 nodes) non-equidistant grids used. For getting the correct (converging) solution, it is important to properly represent the change in concentration gradients in space and time. Changes of concentration are faster at early times near the high concentration boundary. Therefore, results for OGS can be improved by using a nonequidistant grid with a better spatial resolution near the high concentration boundary instead of an equidistant finer grid.  Figure 20 shows the cesium breakthrough curves calculated at different locations (x = 1, 5, 7 and 9 mm) by using different mesh sizes and  . 19 Cs breakthrough curves calculated with OpenGeoSys at different locations of the singlespecies model by using three different spatial discretizations in x direction. The cesium concentration at the "high" concentration boundary was assumed to be 10 −3 mol/L (top), 10 −5 mol/L (middle) and 10 −7 mol/L (bottom). The straight parts of some breakthrough curves are due to output data intervals assuming cesium concentration levels of 10 −3 , 10 −5 and 10 −7 mol/L at the "high" concentration boundary. In all the cases, the effect of spatial discretization is especially relevant near the left boundary of the model as observed with the other codes. In fact, cesium breakthrough curves calculated at x = 5, 7 and 9 mm with different mesh sizes in the range of 0.008-0.04 cm hardly differ. One of the particularities of PHREEQC is that cell-centered concentrations are calculated when transport calculations are performed. Concentrations to other distances of the cell can not be interpolated by the code, which means that adjusted spatial discretization were initially selected in order to have calculated concentrations at the exact locations defined in this benchmark (i.e. 1, 5, 7 and 9 mm). For example, the standard discretization with 25 cells (Δx = 0.04 cm) has its 3rd cell calculations at exactly the distance of 1 mm, however this is not the case with 50 cells (mess size Fig. 20 Cs breakthrough curves calculated with PHREEQC at different locations by using different spatial discretizations in x direction. The cesium concentration at the "high" concentration boundary was assumed to be 10 −3 mol/L (top), 10 −5 mol/L (middle) and 10 −7 mol/L (bottom) fine_2), when the center of the 5th cell is at a distance of 0.9 mm and the center of the 6th cell is at a distance of 1.1 mm. Therefor there is a clear mismatch of 0.1 mm with the exact position considered for the comparison with other codes in this benchmark. This mismatch is also present for the discretization with 40 and 100 cells (see Table 9 for details). However, as can be seen in Fig. 20, this mismatch of the exact position when 40, 50 and 100 cells are used have no influence at 5, 7 and 9 mm and it is minimum at 1 mm, with maximum differences of 4%, only when concentrations <10 −6 mol/L are calculated.
On the other hand, PHREEQC results show a difference of 10% (slightly higher cesium concentration) for the breakthrough calculated for 9 mm for cesium boundary concentration levels of 10 −3 and 10 −5 mol/L for later times when a quasi-steady state is reached when standard discretization (Δx = 0.04 cm) is used. Differences are attributed to the required new implementation of the Dirichlet boundary conditions (Cs =10 −10 mol/L) in the high concentration side (see input file in the supplementary information and Appendix 4 for more details). The effect of the selected boundaries in the modeling results of small thorough diffusion systems, like the one described in this benchmark, has been clearly shown in the literature [71]. Large diffusive fluxes across the boundary, linearly related to a given concentration gradient, makes that even very low concentration differences not taken adequately into account in the downstream boundary give considerable errors [71] which can be reduced when refining the spatial discretization close to the boundaries Due to all the explanations above and considering that simulation with finer discretization (Δx = 0.008 cm), increase the computational time four times and create very large output data files making difficult the visualization of the data in post-processors, it was decided to use the finer discretization of PHREEQC (Δx = 0.025 cm) to compare the results with the other codes Table 9 Mismatch of the exact position defined in this benchmark (i.e. 1,5,7,9) vs discretization used, due to cell centered concentration calculations of the code PHREEQC Discretization (number of cells) 25 40 50 100 125 Mismatch with the exact position defined in the benchmark (i.e. 1, 5, 7 and 9 mm) 0 mm 0.125 mm 0.1 mm 0.05 mm 0 mm "dummy boundary" has no influence in the numerical results. The chemical composition of the dummy boundary cells were defined as Opalinus clay pore water only (see Table 3) with a fixed concentration of Cs of 10 −10 mol/L. In order to fix the Cs concentration and keep it constant during the simulation a "Fix_Cs", fictitious solid phase was created as.

Fix_Cs Cs+ = Cs+ log_k 0
and equilibrated with the Opalinus Clay solution with the Keyword EQUILIBRIUM_PHASES in the two "dummy boundary cells "until the cesium activity of 10 −10 mol/L (log activity = −10) was reached:

EQUILIBRIUM_PHASES 26-27
Fix_Cs -10 In this approximation, it is assumed that the activity of Cs in Opalinus clay is equal to its concentration due to the low ionic strength of the system (see Table 3). However, it is true that at the ionic strength of Opalinus Clay (I = 0.36) and considering the Davies equation [64] for ionic strength corrections, the calculated concentration of Cs is fixed to 1.37 · 10 −10 mol/L instead of 1 · 10 −10 mol/L, though the convergence tolerance (see section 3.5) of the numerical solver in the case of PHREEQC was fixed to 10 −8 which means that this small difference in the fixed concentration is not causing additional numerical errors in the results as can be seen in Fig. 22. SO-093) within the iCross project for partial funding. Partial financial support for the PSI contribution was provided by the National Cooperative for the Disposal of Radioactive Waste (Nagra), Wettingen, Switzerland. We thank the comments and corrections of the reviewers who contributed to greatly improve the paper.
Funding Open Access funding provided by Lib4RI -Library for the Research Institutes within the ETH Domain: Eawag, Empa, PSI & WSL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.