On the concentration-dependent diﬀusion of sorbed cesium in Opalinus Clay

Cation diﬀusion coeﬃcients in clayey materials partly appear to be greater than diﬀusion coeﬃcients of water tracers. The measured values vary between experiments performed at diﬀerent salinities or diﬀerent tracer concentrations. This eﬀect is especially pronounced for cations that sorb strongly on the clay surfaces, such as Cs. The observations illustrate the diﬃculties in applying Fick’s law to cation diﬀusion in clays and demonstrate the need to ﬁnd a consistent description of cation diﬀusion in clays that can be used to predict experiments performed at diﬀerent conditions. In order to consistently describe Cs diﬀusion in Opalinus Clay, a multi-site surface diﬀusion model was implemented in the continuum-scale reactive transport code Flotran. The model combines pore and surface diﬀusion in one single diﬀusion coeﬃcient, which accounts for the diﬀusion of sorbed cations along the clay surfaces.. The contribution from surface diﬀusion to the diﬀusion coeﬃcient is directly coupled to the sorption behavior via the derivative of the sorption isotherm. The model parameters include the surface mobilities, which are speciﬁc for each cation and sorption site. To derive surface mobilities for Cs, in-diﬀusion experiments were conducted at eight diﬀerent stable Cs background concentrations. A set of surface mobilities for Cs on three sorption sites in Opalinus Clay was estimated by ﬁtting the surface diﬀusion model simultaneously to these experimental data. Moreover, the sensitivity of the model to sorption parameters and surface mobilities was evaluated. The surface diﬀusion model with the estimated surface mobilities was then successfully tested against independent experimental data for Cs in Opalinus Clay, illustrating the model’s predictive capabilities.


INTRODUCTION
There is a worldwide interest to use clays and claystones as sealing components for waste disposal sites. In Switzerland, Opalinus Clay (OPA), an argillaceous rock, is considered as potential host rock for a radioactive waste repository (Nagra, 2002). The low hydraulic conductivity (10 -10 -10 -15 m s À1 ) restricts migration of solutes through clay formations, such as OPA, essentially to diffusive pro-cesses only (Gimmi et al., 2007;Mazurek et al., 2011;Yu et al., 2018). The sorption capacity of clay for various cations or contaminants of interest as well as the detailed mechanisms of the sorption process have been investigated in numerous studies (Elprince et al., 1980;Neal and Cooper, 1983;Shainberg et al., 1987;Fletcher and Sposito, 1989;Helios et al., 1995;Staunton and Roubaud, 1997;Tournassat et al., 2007;Robin et al., 2015;Montoya et al., 2018;Siroux et al., 2018;Fernandes and Baeyens, 2019). For Cs, sorption on clays has been found to be non-linear in a complex way, which was attributed to the existence of different sorption sites. Various multisite sorption models (Poinssot et al., 1999;Bradbury and Baeyens, 2000;Savoye et al., 2012;Benedicto et al., 2014; https://doi.org/10.1016/j.gca.2021.01.012 0016-7037/Ó 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Chen et al., 2014;Missana et al., 2014a;Missana et al., 2014b;Cherif et al., 2017) have been developed in order to characterize the sorption behavior of Cs onto different clays. Most of these models considered cation exchange as the main sorption process. It was generally found that soils and sediments have a comparably small number of high affinity sites where Cs sorbs preferentially, and a larger number of sites with lower affinity for Cs (Fuller et al., 2014). Often, the high affinity sites were attributed to frayed edges of illite (FES), whereas the sites with lower affinity were attributed to planar sites (PS), that is, to sites on external or internal basal planes of illite and smectites, with the charge arising from isomorphous substitution.
In Opalinus Clay, Cs sorption can be described by a generalized three-site ion exchange model Van Loon et al., 2009). This model is based on the sorption properties of illite (Poinssot et al., 1999) but also includes planar sites of other clay phases (e.g., illite/ smectite mixed layers) according to the measured cation exchange capacity (CEC). Following Poinssot et al. (1999), Type II sites (TTS) with intermediate capacity and strong to intermediate affinity for Cs are distinguished in addition to FES and PS. Zachara et al. (2002) used a two-site model, similar to Poinssot et al. (1999), to describe Cs sorption in Handford sediments. Steefel et al. (2003) modified this model by adding a second FES for a more sensitive fit. Chen et al. (2014) applied the generalized three-site model (FES, TTS, and PS) successfully to samples from the Callovo-Oxfordian (COx) clay-rich formation with differing mineral compositions. They concluded that only PS should be attributed to illite/smectite mixed layer minerals. More recently Cherif et al. (2017) compared sorption of Cs on many different clay minerals and developed a model with two sites (one ion exchange site, one surface complexation site) each on illite, montmorillonite, and kaolinite. As different sorption affinity parameters (selectivities, equilibrium constants) were derived for each of these minerals, this essentially represents a six-site model for a multi-mineral clay. Surprisingly, they considered a pHdependent SC model for frayed edge sites on illite, which is at odds with the widely accepted idea that Cs has a highly specific affinity to these sites due to its large size and low hydration energy (Sawhney, 1972;Brouwer et al., 1983;Poinssot et al., 1999;Bradbury and Baeyens, 2000). In any case, the distinction of different sites in the aforementioned multi-site sorption models is empirically based on the shape of the adsorption isotherms, mostly without any evidence from spectroscopic methods. This is especially true for the low capacity sites.
Model concepts for cation sorption have also been developed by considering interactions between a charged, planar surface and the surrounding electrolyte solution. In the classical Gouy-Chapman diffuse double-layer model (Huang and Stumm, 1973;Singh and Uehara, 1999), the negative surface charge (representing the first layer of the model) is compensated by a swarm of electrostatically attracted cations, with decreasing concentration with distance from the surface such that electrostatic and thermal forces are balanced. The modified Gouy-Chapman model (Carnie and Torrie, 1984;Sposito, 1992;Tournassat et al., 2009) accounts for the fact that ions have a finite size; the origin of the diffuse layer is thus shifted by one ion radius from the charged surface. The Gouy-Chapman-Stern model (Bowden et al., 1977;Westall and Hohl, 1980) includes an innermost compact layer (Stern layer) with more specifically adsorbed cations, followed by the diffuse layer; resulting in a triple-layer model (cf. Fig. 1). The Grahame model or other triple-layer models (Davis et al., 1978;Leroy and Revil, 2004;Leroy et al., 2007) also assume a specific layer of cations next to the surface. The layer of sorbed cations next to the surface in these models can be considered as Stern layer, without referring explicitly to the strength or type of interaction. This surface layer may consist of cations forming inner-sphere complexes (i.e., without an intervening water molecule) or outer-sphere complexes. Cations in the surface layer and in the diffuse layer may be considered as sorbed in the sense that they compensate the surface charges. One should keep in mind, however, that the diffuse layer typically contains an excess of cations, depending on the equilibrium electrolyte concentration. Far away from the surface, the ion concentrations reach the values of this electrolyte concentration. Thus, cations can generally occur in the two counter layers (Stern or surface layer, diffuse layer) as well as in so-called 'free' pore water, which is the charge neutral electrolyte solution (Fig. 1).
The traditional approach to model cation diffusion in clays is based on the simple Fick's law. In this approach, only cations in 'free' pore water contribute to diffusion, while any sorbed cations are considered as immobile, with Fig. 1. Schematic representation of the concept of a triple-layer model: Cations are sorbed (light blue domain) on the negatively charged clay surfaces (grey) in the Stern layer (SL) and the diffuse layer (DL) close to the clay surface. The latter contains an excess of cations. In the bulk water (or 'free' pore water; dark blue domain) ions are not affected by electrostatic forces resulting from the negatively charged clay surfaces. Each domain may contribute to the overall mass flux j in a surface diffusion model. Here, we only distinguish between sorbed cations (red frame) and cations in the pore water, with an (average) surface mobility of sorbed cations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) no distinction between different sites or types of sorption. This means that sorbed cations, including those in the DL, do not add to the overall mass flux. However, this approach often appears to be problematic when evaluating cation diffusion experiments in clays. Especially for strongly sorbing ions, such as Cs, diffusion coefficients are often higher than diffusion coefficients of neutral tracers (e.g. tritiated water, HTO). For example, Jakob et al. (2009) found the effective diffusion coefficient of Cs in Opalinus Clay to be one order of magnitude greater than that of HTO, Melkior et al. (2005) found similar results for Callovo-Oxfordian clay, whereas Wersin et al. (2008) found a 3 times higher and Van Loon et al. (2004b) a four times higher effective diffusion coefficient for Cs than for HTO in Opalinus Clay. These observations are inconsistent with simple Fick's law.
One explanation for this phenomenon is the excess of cations in the diffuse layer which leads to an enhanced mass flux. However, models including DL diffusion on the basal planes of the clay platelets also cannot describe the diffusion behavior of Cs in OPA consistently, and bulk diffusion coefficients still have to be adapted to the experimental data (Appelo et al., 2010;Soler et al., 2019). Tachi and Yotsuji (2014) found the modeled Cs diffusion coefficient at an ionic strength of 0.01 M about matching the experimental value when applying a DL diffusion model to Cs diffusion in compacted sodium montmorillonite, while it underestiamted the values at higher ionic strengths.
A more general explanation for the discrepancy between experimental results and classical Fickian diffusion theory is the movement of any sorbed cations along the clay surfaces, without referring to a specific type of sorption (Jenny and Overstreet, 1939;van Schaik et al., 1966;Cheung, 1990;Berry and Bond, 1992;Eriksen et al., 1999;Gimmi and Kosakowski, 2011). This process is called surface diffusion, with the term referring to the diffusion of any sorbed cations, including those in the DL. Although DL cations are mobile, Stern layer cations are generally considered to be immobile. However, some mobility has been attributed to these cations in order to be able to match conductivity or electrophoretic data (Lyklema et al., 1998;Revil et al., 1998;Lyklema, 2001;Leroy and Revil, 2004;Weber and Stanjek, 2017). From a mechanistic perspective, surface diffusion can be considered as ions moving from one adsorption site to another (Jenny and Overstreet, 1939;Revil et al., 1998;Lyklema, 2001). Gimmi and Kosakowski (2011) compiled a large set of cation diffusion data in clays and soils in order to evaluate the mobility of sorbed cations. They found a distinct surface mobility for each cation. However, Cs data were concentration-dependent, and different surface mobilities were given for low and high Cs background concentrations. The concentration-dependent mobility was attributed to the different sorption sites for Cs, which dominate at different Cs concentrations. Therefore, they proposed a generalized multi-site surface diffusion model. This model considers a (partial) mobility of the cation on each of the sorption sites or domains.
In this study, we present the application of a multi-site surface diffusion model to Cs diffusion data in Opalinus Clay. The model was implemented in the reactive transport code Flotran (Lichtner, 2007). Specific in-diffusion experiments at different Cs background concentrations were conducted in order to estimate the sorption-site-specific surface mobilities of Cs. The sensitivity of simulations to model parameters was evaluated and the model was tested against other Cs diffusion data in Opalinus Clay.

Formalism of the multi-site surface diffusion model
At the pore scale, surface diffusion occurs in parallel to diffusion in the bulk pore water (Gimmi and Kosakowski, 2011). Therefore, the total cation mass flux can be written as the sum of the contributions from pore diffusion, described by Fick's first law, and surface diffusion: where j p and j s are the fluxes for the pore and surface regions ( Fig. 1), C is the concentration in solution, S is the sorbed concentration per mass of dry solid, D 0 is the cation diffusion coefficient in bulk water, D s0 is the intrinsic surface diffusion coefficient for a flat surface, e is the total (water filled) porosity, q bd is the bulk dry density, s p and s s account for the tortuous pathway in the pore water and on the surface, respectively, and @C=@x and @S=@x are the concentration gradients in the pore water and on the surface. The balance equation is: Neglecting the second term on the right hand side in Eqs. (1) and (2) results in the classical Fickian diffusion model with the effective diffusion coefficient D e being D e ¼ eD p ¼ eD 0 =s p , and the pore diffusion D p ¼ D 0 =s p . Eq. (1) can also be used to describe the flux at the sample scale as long as the tortuosities are representative for this scale. If we assume local equilibrium between pore and surface concentrations, pore and surface diffusion fluxes can be combined in a single term that depends on the concentration gradient @C=@x only, using the relation @S=@x ¼ @S=@C ð Þ@C=@x ð Þ . This relationship is, however, only valid if S depends uniquely on C, as is typically the case for trace concentrations where the background concentrations of other cations participating in sorption remain about constant. More generally, the surface concentration gradient of a cation k, @S k =@x, can be given as: where @S k =@C j are the derivatives of the sorbed concentration of cation k with respect to the concentrations of the cations j in solution and @C j =@x are the concentration gradients of the cations j. Here, we focus on the case of constant background cations where @C j-k =@x % 0, but the surface diffusion model was implemented according to the general form of Eq. (3).
Substituting Eq. (3) into Eq. (1), the total flux of the surface diffusion model depends then on a single, combined diffusion coefficient. If multiple sorption sites i are present, we have In the above relation, l sk;i is the relative surface mobility of the cation k on site i. When species concentrations or solution composition are not constant over the domain or time span of interest, the local derivative terms, and therefore D e;comb , become space-and time-dependent. As different D e;comb have to be used for different ions in the surface diffusion model, a multi-species diffusion approach is required. Thus, diffusion has to be described by the Nernst-Planck equation, which is implemented in Flotran. The relative surface mobility l sk;i is equal to the ratio D s0;i =D 0 . In other words, it is a measure for the average diffusive mobility of a cation sorbed on site i on a flat surface in comparison to its diffusive mobility in water. For constant background concentrations (@C j-k =@x % 0 or @C j-k =@C k % 0), this relation simplifies to (Gimmi and Kosakowski, 2011): where the subscript k for the sorbing cation of interest has been dropped. Gimmi and Kosakowski (2011) assumed that s s ¼ s p , which was justified in view of the closeness of surface and pore diffusion pathways. Also, the effective cation porosity is set equal to the total (water-filled) porosity e ¼ e w , because cations and water occupy the same pore space (external and interlayer pores). The first term in Eqs.
(4) and (5) stands for diffusion in the pore water, while the second term accounts for diffusion along the clay particle surfaces. Surface diffusion may occur on each of the available sorption sites i for a certain cation, according to the sorption model and the value of the surface mobility l s;i . The multi-site surface diffusion coefficient sums up the contributions of all sorption sites i to the overall mass flux.

Cs sorption model for Opalinus Clay
The sorption behavior of Cs in OPA can be described by the generalized three-site sorption model (GCS) of Bradbury and Baeyens (2000), where Cs sorption is linked to the illite mineral content of the clay. The site-specific derivatives of the sorption isotherm @S i =@C needed in Eq. (5) (or the more complicated sum term over all cations j in Eq. (4)) are calculated from the GCS. The calculation of these terms based on a general cation exchange reaction is explained in the Appendix A. The GCS proposes that Cs sorbs via cation exchange reactions onto three different sorption sites: Planar sites (PS), Type II-sites (TTS) and Frayed-edge sites (FES). Frayed-edge (0.25%) and Type II (20%) site capacities were fixed according to the illite mineral content in the OPA specimen and a reference illite with a CEC of 0.2 eq kg À1 . The planar site capacity was taken to be 0.095 eq kg À1 for Mont Terri Opalinus Clay (Bradbury and Baeyens, 1998), which accounts for the sum of all clay minerals (e.g. illite-smectite mixed layers) participating in Cs sorption. TTS and FES site capacities were calculated to 9.2 Á 10 -3 eq kg À1 and 1.15 Á 10 -4 eq kg À1 , respectively, taking an illite content of 23%. The mineralogical and pore water composition of OPA are summarized in Table 1. Selectivity coefficients (uncertainty of ±0.2 log K units) for the cation exchange reactions are given by Bradbury and Baeyens (1998) Table 2).

Experimental setup
Diffusion experiments were performed in order to determine the surface mobilities of Cs in Opalinus Clay. The idea was to determine Cs tracer diffusion at several different, about constant stable Cs (Cs-133) background concentrations. In each experiment, Cs tracer diffusion was expected to be dominated by sorption and surface diffusion of certain sites only. For this purpose, small rock columns were presaturated with a given Cs background concentration. Then, Table 1 Mineralogical composition (range according to Bradbury and Baeyens (1998) and Lauber et al. (2000)) and pore water composition of OPA.

Mineralogical composition
Pore water composition in-diffusion of a Cs tracer (Cs-134) into the same column was monitored. The Opalinus Clay specimen used in these experiments were from the Mont Terri Rock Laboratory (St. Ursanne, Switzerland). The specimen originate from the borehole drilled for the DI-A field experiment (Wersin et al., 2004). The setup of the experiment is similar to that of Van Loon and Mü ller (2014) (Fig. 2). Small cylindrical specimen were subcored parallel to the bedding plane from a single larger core. The small bore cores were embedded in an epoxy resin (Epofix, Streurs, Germany). After hardening, the resin at one front side of the core was removed. This end of the specimen was contacted for 2 months at 25 ± 1°C with 100 mL synthetic pore water suited for OPA at the sampling site (Table 1) containing different concentrations of stable cesium between 10 -2 M and 3.5 Á 10 -8 M (Table 3). Simplified calculations showed that after 2 months clearly different levels of Cs-133 concentrations are reached in the clay specimen close to the reservoir. The geometrical parameters of the specimens are given in Table 4.
After 2 months of resaturation and equilibration, the specimens were contacted with 100 mL solution of similar composition but spiked with Cs-134 tracer (Table 3). After 7 days in-diffusion, the profiles of Cs-134 in the OPA specimens were analysed using the abrasive technique described in Van Loon and Mü ller (2014). The abrasive paper used was P400 (Siawatt fc, SIA abrasives, Switzerland). The thickness of the sampled layers varied between $100 mm and 200 mm. The uncertainty on the diffusion distance was calculated as described in Van Loon and Mü ller (2014). The total concentration of Cs-134 per dry mass of the removed rock was determined by c-counting.

General model parameters
Diffusion was modeled by applying species-dependent pore diffusion coefficients for all ions, while the contribution of surface diffusion was only added to the Cs diffusion coefficient. Aquatic complexes and competitive sorption were taken into account by considering geochemical equilibrium in solution and between solution and clay surfaces using the activity model of Davies (Eq. (A2)) (Stumm and Morgan, 1996).
A Cs diffusion coefficient in bulk water D 0;Cs of 2 Á 10 -9 m 2 s À1 (Flury and Gimmi, 2018) was used for all simulations. The porosity e and the bulk dry density q bd were calculated from the geometric in-diffusion data (Table 4) using a grain density of OPA of 2.7 g mL À1 . The tortuosities s p and s s were approximated with s w for a water tracer (Gimmi and Kosakowski, 2011) taking the same value of 6.14 for all experiments from Appelo et al. (2010). At the interface between reservoir and OPA specimen a constant concentration boundary condition was taken, whereas a no flux boundary was used to describe the boundary at the other end of the specimen.

Sorption model
All experiments were modeled using the given pore water composition and a single sorption isotherm according to the generalized Cs sorption model of Bradbury and Baeyens (2000). The default parameters for sorption (site capacities, selectivities) were those from the literature given in Section 2.2. The Cs-134 concentration in the reservoir was approximately constant during the experiment (less Table 2 Cation exchange reactions and selectivity coefficients used in the GCS.

Specimen
Cs-133 (M) Cs-134 Bq dm À3 Cs-134 (M) 1 1 0 -2 10 7 1.56 Á 10 -9 2 1 0 -3 5 Á 10 6 7.70 Á 10 -10 3 1 0 -4 3 Á 10 6 4.68 Á 10 -10 4 1 0 -5 2 Á 10 6 3.11 Á 10 -10 5 1 0 -6 3 Á 10 5 4.68 Á 10 -11 6 3 Á 10 -7 3 Á 10 5 4.68 Á 10 -11 7 1 0 -7 10 5 1.56 Á 10 -11 8 3.5 Á 10 -8 10 5 1.56 Á 10 -11 than 1% of the initial mass of Cs-134 in the reservoir diffused into the specimen). The comparison of the reservoir Cs-134 concentration with that of the first data point of the in-diffusion profiles, which was adjacent to the reservoir solution, allows an independent estimation of the Cs-134 sorption at different Cs background concentrations. These points were compared with the default sorption isotherm based on the literature data, and site capacities were adapted when considered necessary. To do so, the stable Cs sorbed concentration was calculated from the corresponding stable Cs aqueous background concentration by using the distribution coefficient derived from the Cs-134 aqueous and sorbed concentrations. The natural heterogeneity of OPA at the centimeter scale can lead to local variations of the illite or smectite contents, which then directly affect the local sorption capacities of the different specimen drilled from the same core sample. Furthermore, selectivity coefficients were also modified -within their given uncertainty -in order to better match the isotherm data. As stated, a single sorption isotherm was finally used for all eight experiments, but it may also be justified to use slightly different site capacities for each separate specimen in response to local heterogeneity.

Surface mobilities
The surface mobilities were estimated by fitting the surface diffusion model visually to the experimental indiffusion profiles. A sequential procedure was used. At the highest C Cs;bg of 10 -2 M, FES and TTS are mostly occupied by stable Cs and @S FES =@C and @S TTS =@C tend to be negligibly small. Therefore, only the mobility on the planar sites significantly contributes to the overall mass flux. Accordingly, the data at this highest background concentration are most sensitive to l PS , and the value of this parameter was determined from these profile data and kept constant for the modeling of the other experiments. Second, the value for l TTS was estimated mainly based on the data at intermediate to high concentrations, where @S FES =@C still is negligible. Finally, the l FES value was derived mainly with a focus on the data at the lowest background concentration of 3.5 Á 10 -8 M, leading in the end to a single set of surface mobilities.

Sensitivity to model parameters
Some correlation exists between the three surface mobilities, as well as between these parameters and the isotherm parameters. In order to get a feeling for these correlations, three sets of sensitivity simulations were performed, where the effect of a variation of the sorption capacities, the selectivities, and the surface mobilities were investigated. These sensitivity simulations served also to check effects of possible specimen heterogeneities with respect to these parameters.

Global parameter estimation
In order to better judge the validity of the sorption and diffusion parameters, the site capacities, selectivity coefficients and surface mobilities for Cs (in total 9 parameters) were also estimated with a global optimization procedure. To do so, the open-source library NLopt (Johnson, 2019) for non-linear optimization was used. The sum of the squared error between model and experimental data (diffusion profiles) was minimized applying the ESCH algorithm (da Silva et al., 2010). This is an evolutionary algorithm for global optimization implemented in the NLopt library. The squared errors between model and data values were weighted by the corresponding squared uncertainty (error) of the diffusion data points. The optimization procedure was limited to 20,000 iterations to find the global optimum of all parameters simultaneously for seven out of the eight experiments (excluding that at C Cs;bg 10 -7 M, which has some peculiarities, see below).

Modification of the sorption isotherm
The usage of the original (default) isotherm parameters (Section 3.2.2) was not satisfactory, especially for low stable Cs background concentration C Cs;bg ; sorption as derived from the first three profile points was overestimated for most experiments (Fig. 3). Therefore, the site capacities as well as the Cs Na K c values were slightly modified to better match the experimental data, the latter within the range of their uncertainty (Table 5). The CEC of the PS and TTS had to be decreased by $25%, compared to the original isotherm , and the CEC of the FES by $50%, respectively, in order to approximately match all eight experiments simultaneously (Fig. 3). The modification significantly lowered the sorbed amount of Cs at low C Cs;bg . However, sorption remained slightly overestimated for C Cs;bg of 3.5 Á 10 -8 M, while for the experiment with C Cs;bg of 10 -7 M the large difference between the first three data points, which may hint to heterogeneity within the specimen, precludes a unique match anyway. The modified isotherm parameters were then used for the simulations carried out to estimate the surface mobilities. The used specimen were small (Table 4). Thus, the scatter of the data with respect to a single isotherm may be partly due to heterogeneity between the specimen (or heterogeneity within the specimen for C Cs;bg of 10 -7 M). Heterogeneity is expected to affect mainly site capacities through local variations of mineralogical composition. Selectivities, even though not representing real thermodynamic parameters, are considered as more fundamental than capacities, but it is clear that they also can depend on mineralogical variations.
One has to keep in mind that the sorption isotherm affects also the effective diffusion coefficient in a surface diffusion model, because of the coupling to its derivative. Accordingly, a variation of the sorption parameters will also influence the respective diffusion profiles.

Surface mobilities determined from in-diffusion profiles
From the data at the highest C Cs;bg of 10 -2 M, a value of 3.3 Á 10 -2 was estimated for l PS ; best agreement with the experiments at intermediate to high concentrations was found when using a mobility of zero for the Type II sites, whereas a mobility on the FES of 1.5 Á 10 -3 was necessary to fit the data for low C Cs;bg . Fig. 4 shows the experimental and modeled Cs-134 diffusion profiles for the eight experiments. Reasonable matches between experimental data and model were obtained for the high C Cs;bg of 10 -2 M, 10 -3 M, 10 -4 M and 10 -5 M (with the first point appearing as an outlier at 10 -5 M). For 3.5 Á 10 -8 M Cs, sorption near the inlet (first three points) was slightly overestimated, whereas for C Cs;bg of 10 -6 M it was slightly underestimated near the inlet and the diffusion profile was too flat. For 3 Á 10 -7 M C Cs;bg the model matches the experimental data only near the inlet but then overestimates the sorbed Cs-134 concentration (except for the last point). Finally, a reasonable match for C Cs;bg 10 -7 M could not be achieved. At this C Cs;bg , the first measurement point exhibits a very high concentration compared to the rest of the profile. This behavior cannot be described appropriately by diffusion into a homogeneous rock specimen.

Global parameter estimation
The globally estimated model parameter values are in good agreement with the values obtained by the visual fit (Table 5). Site capacities of PS and FES from the global fit were found to be about 20% lower, those for the TTS about 10% higher than those from the visual fit. Both fitting procedures led to the same selectivity coefficients. The different site capacities result in a slightly different Cs sorption isotherm (Fig. 3). At low C Cs;bg of 3.5 Á 10 -8 M to 3 Á 10 -7 M the global fit tends to better match the isotherm data, whereas for 10 -6 M sorption is underestimated. At midrange concentrations both modified isotherms are about identical and for highest C Cs;bg sorption is slightly lower for the global fit. The surface mobilities were found to be nearly identical on the PS and about 20% lower on the FES compared to the visual fit. A very small mobility of 5.3 Á 10 -4 was estimated for the TTS. The diffusion profiles for the global fit exhibit only small, hardly significant differences compared to those of the visual fit (Fig. 4). For all C Cs;bg the global fit parameters lead to slightly lower total Cs-134 concentrations over the diffusion distance. The differences of both fits increase with decreasing C Cs;bg . These results show that a global fitting routine can be successfully applied to optimize the parameters for a complex diffusionsorption behavior. In addition, the results corroborate the validity of the step-wise visual fit procedure, on which the following discussion will focus.

Sensitivity of the surface diffusion model to sorption parameters
Because of the direct coupling of the effective diffusion coefficient of the surface diffusion model to the sorption isotherm through the derivative term @S=@C, the diffusion behavior is sensitive to the choice of sorption parameters. In order to evaluate the influence of the sorption parameters on the simulated diffusion profiles, the CEC of the sorption sites was varied by ±40% and the selectivity coefficients by ±0.2 log Cs Na K c units (Figs. 5-7). A positive variation of the sorption parameters leads to higher effective diffusion coefficients and vice versa (Eqs. (5), (A5) and (A6)).  Table 5); Experimental data: black and grey dots with error bars for the first three data points of the measured indiffusion profile. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) A variation of the CEC of each site appears as a shift of the isotherm in a log-log plot along the y-axis (Fig. 5a). This leads to identical trends in the profiles at all C Cs;bg (Fig. 6). At the specimen-reservoir interface, where the OPA specimen is in equilibrium with the constant concentration in the reservoir, the difference between the base case and the upper and lower variation is equal to ±40%. With increasing distance from the interface, the variations approximate the base case faster or slower depending on the C Cs;bg .
The variations of the sorption isotherm parameters affect the diffusion profiles differently in the surface diffusion model compared to a classical Fickian diffusion model having a single, constant diffusion coefficient. For a model with a constant diffusion coefficient a positive variation of CEC and/or Cs Na K c leads to a steeper diffusion profile, a negative variation to a flatter profile. For the surface diffusion model, the same trend appears when varying the CEC. However, the trend is less pronounced, because it is partly compensated by the larger effective diffusion coefficient calculated for a higher CEC (and vice versa) (see Eqs. (5), (A5) and (A6)).
The upper and lower variation of Cs Na K c lead to a different behavior. An increase or decrease of all Cs Na K c tends to shift the isotherm in a log-log plot along the x-axis to the left or to the right, respectively (Fig. 5b). Thus, little change in sorption (and thus in the effective diffusion coefficient) occurs in regions where the isotherm is flat (e.g., near saturation of a site) as compared to steep regions. As a consequence, the effect of a variation of all Cs Na K c on the simulated profiles clearly varies depending on C Cs;bg (Fig. 7). It does so in a complex manner because of the superposition of the three sorption sites and because C Cs;bg is not fully constant but decreases slightly along each diffusion profile. Therefore, depending on C Cs;bg , the variations depart from the base case (10 -3 M, 10 -5 M, 10 -6 M, 3 Á 10 -7 M, 10 -7 M) or approximate to the base case (10 -2 M, 10 -4 M, 3.5 Á 10 -8 M) with increasing distance from the interface.
The variable influence of the sorption parameters on the shape of the diffusion profiles over the entire penetration distance clearly illustrates the difficulty to find an appropriate set of CEC and Cs Na K c . We limited the range of selectivity coefficients to the given uncertainty of ± 0.2 log K units, and required the simulation to match the first data point of the profile within its uncertainty. It is possible that other sets of sorption parameters (and correspondingly other values for the surface mobilities) exist that describe the data similarly. However, it appears justified to restrict the values of the sorption parameters to the range of well-established data from the literature, which then restricts the values of surface mobilities to a relatively small range. In any case, the set of sorption parameters that was finally used led to an appropriate fit of the diffusion profiles with the estimated surface mobilities.
Because of the small size of the OPA rock specimens (diameter: 2.6 mm; length: 12 mm), the mineralogical composition, and therefore the illite/smectite content as well as the selectivities, can vary from specimen to specimen or within a specimen due to the natural heterogeneity of Opalinus Clay at the centimeter to (sub)millimeter scale. From Figs. 6 and 7 we can see that local variations of CEC and/or Cs Na K c could indeed be responsible for the scatter between the specimens or within a single specimen.

Sensitivity of surface diffusion model to surface mobilities
The site-specific surface mobilities determined from the experimental data were varied by ±20% in order to evaluate their influence on the diffusion profiles. At the specimenreservoir interface the variations do not differ from the base case (Fig. 8), because the reservoir volume (100 mL, ratio to pore volume $13 0 000) was large enough such that the reservoir concentration of Cs remained approximately constant, independent of the used l si . Larger and smaller l si lead to flatter and steeper profiles, respectively. Accordingly, the variations diverge with increasing penetration distance for all experiments. The surface mobilities are constant model parameters, such that their variation by ±20% leads to a 20% higher or lower contribution of the surface diffusion term to the effective diffusion coefficient. The effect is thus similar to a corresponding variation of the constant diffusion coefficient in the classical Fickian diffusion model, but not fully equivalent, because of dependence of D e;comb in our model on the slightly variable background Cs concentrations.

Average mobilities and comparison to existing data
One set of site-specific surface mobilities was found to describe the Cs diffusion behavior in Opalinus Clay. There are no reference data against which the present set of surface mobilities can be directly compared. However, it is possible to calculate the average mobility of Cs in OPA by weighting the site-specific mobilities with their average sorption capacity ratio. The average mobility can then be directly compared to the surface mobilities determined for  [-] 3.3 Á 10 -2 3.4 Á 10 -2 $0 5.3 Á 10 -4 1.5 Á 10 -3 1.2 Á 10 -3 a large set of diffusion data in clays in Gimmi and Kosakowski (2011). The average mobility can be written as (Gimmi and Kosakowski, 2011): with j ¼ P j i and j i being the sorption capacity ratio for site i: Table 6 shows the average mobilities for the different indiffusion experiments. For high C Cs;bg , the contribution of surface diffusion is dominated by the mobility on the planar sites (l PS = 3.3 Á 10 -2 ), as expected. The calculated average mobilities for high Cs are in good agreement with the range of 1-4 Á 10 -2 determined in Gimmi and Kosakowski (2011). This finding confirms the validity of the estimated mobility on the planar sites. Kosakowski et al. (2008) performed molecular dynamics simulations determining Na and Cs diffusion in montmorillonite. They found values of l PS for monohydrated interlayers similar to ours ($0.002-0.1). This is plausible, because considering the high bulk dry density of OPA ($2.3 g mL À1 ), mostly monohydrated interlayers are expected in OPA. With the estimated l PS and the derived (close to) zero mobility of Cs on the TTS, the experiments can be described down to a C Cs;bg of 10 -6 M. A larger mobility on the TTS would result in a worse match of the data for the mid-range concentrations. The average mobilities derived here for Cs at trace conditions ($0.002 at very low Cs concentrations) in OPA are significantly larger than average values (2 Á 10 -4 -4 Á 10 -4 ) found by Gimmi and Kosakowski (2011), but similar to their maximum values. The surface mobility of 1.5 Á 10 -3 attributed here to the FES, which is dominant at low Cs, is comparably large. Such a high surface mobility on FES may not be expected considering the strong sorption of Cs on these sites (high Cs Na K c ), but is clearly required to describe the experimental data at low C Cs;bg (<10 -6 M).
The site-specific surface mobilities may incorporate geometrical effects of the surface topology. By assuming s s ¼ s w we neglect the differences of surface and pore diffusion pathways for Cs. As a result, any deviation of the tortuosities between these pathways is accounted for in the fitted values of the surface mobilities.

Comparison with a classical diffusion model
The experimental data were also modeled with a classical Fickian diffusion model using the modified Csisotherm, in order to demonstrate the differences to a surface diffusion model. For every experiment, an effective diffusion coefficient of Cs D e;Cs was fitted separately to match the corresponding profile data. All fitted D e;Cs were larger than typical effective water tracer diffusion coefficients parallel to bedding, by a factor of 3 to 21 (Table 7; with an effective diffusion coefficient of a water tracer D e;HTO being 6 Á 10 -11 m 2 s À1 ). The fitted D e;Cs show a concentration dependency with increasing values from high to low C Cs;bg . The different values of the diffusion coefficients demonstrate the difficulties using a classical Fickian diffusion model for describing Cs diffusion in OPA. The experiments cannot be described consistently with a single diffusion coefficient.
In contrast, the observations can be easily explained with the surface diffusion model. From C Cs;bg of 10 -2 M to 10 -5 M or 10 -6 M, the relative diffusion coefficient D e;Cs =D e;HTO increases from 3 to 5 (Table 7). At 10 -2 M C Cs;bg planar sites approach saturation. This means that the derivative @S PS =@C is smaller at 10 -2 M C Cs;bg compared to the lower C Cs;bg of 10 -3 M, 10 -4 M and 10 -5 M (Table 6). At the latter concentrations, Cs sorption is about linear for the PS and the value of @S PS =@C is about constant, which is the reason for the constant effective diffusion coefficient in this concentration region. The ratio of @S PS =@C at 10 -2 M and at 10 -3 M to 10 -5 M is $0.6, which equals about the ratio of the relative diffusion coefficients in Table 7. The increasing diffusion coefficients at low C Cs;bg corroborate the findings of a non-negligible surface mobility of Cs sorbed on the FES in the surface diffusion model. Without a mobility of Cs on FESs, the surface diffusion model is not able to reproduce the experimental data or the large relative diffusion coefficient at the lowest Cs background concentrations.

Indications for diffuse and Stern layer diffusion
According to Appelo et al. (2010) 45% of the charge of the total CEC in OPA is compensated in the diffuse-layer. Wigger and Van Loon (2018) found that about 40% of the total porosity of OPA has to be attributed to DL (or Donnan) pores in order to describe anion exclusion effects. Therefore, DL effects certainly play a relevant role in diffusion of charged species. The surface diffusion model here does not directly account for diffusion in the DL; any DL effects are included indirectly in the average surface mobilities. A DL model can predict a dependency of diffusion on the ionic strength of the solution. However, a model including just DL sorption and diffusion cannot account for a dependence of diffusion on the Cs background concentration as observed here, and it would also predict equal effects for equally charged ions such as Cs and Na. In the experiments here, the solution composition with respect to the major ions was identical, except for the different Cs back-ground, which is small compared to the total ion concentration. Then, neither a significant variation of the thickness of the diffuse layer for the different experiments, nor a variation of the diffusion coefficients in the DL, nor a variation of the distribution ratio of Cs between the DL and external pore water is expected. Accordingly, a DL-only model cannot reproduce the observed Cs concentration dependency of diffusion.
Furthermore, in Cs diffusion studies with Cs reservoir concentrations of 10 -3 M (Appelo et al., 2010) and 2 Á 10 -4 M (Soler et al., 2019), where a model with DL diffusion was applied besides specific Cs sorption, bulk diffusion coefficients of Cs still had to be increased. This observation suggests that an additional diffusion process on PSs, which are dominant at these concentrations, contributes to the overall diffusive flux. Appelo et al. (2010) attributed this additional process to interlayer diffusion. The process may generally be related to a mobility of cations sorbed more specifically. Experimental data for Cs sorption onto muscovite (Lee et al., 2012) and molecular dynamics studies of Cs sorption on illite (Lammers et al., 2017) showed that Cs sorbs as inner-sphere complex on the respective basal planes, that is, in the Stern layer. This means that l PS incorporates not only DL diffusion, but also (to some degree) Stern layer diffusion. Such (partial) Stern layer diffusion becomes even more distinct at low C Cs;bg , where FESs dominate sorption. There, a mobility of sorbed Cs on the FES (and therefore a higher effective diffusion coefficient) is required to match the experimental data. At these low C Cs;bg , the proportion of specifically sorbed Cs (with higher selectivity) becomes larger, and Stern layer diffusion (or diffusion of more specifically sorbed Cs) can contribute proportionally more to the overall diffusive flux. However, to evaluate in more detail the contributions from diffuse and Stern layer diffusion, separate modeling of the experiments with a model considering also a DL is necessary.

Test of the model: Application to different Cs diffusion data in Opalinus Clay
The surface diffusion model with the estimated set of surface mobilities was finally tested against data of an independent, radial through-diffusion experiment for Cs in Opalinus Clay (Appelo et al., 2010). The experimental setup is described in detail in Van Loon et al. (2004a). The initial Cs concentration in the high concentration reservoir was 10 -3 M. During the course of 1500 days, the reservoir 10 -2 M 1.8 Á 10 -10 3 10 -3 M 3 Á 10 -10 5 10 -4 M 3 Á 10 -10 5 10 -5 M 3 Á 10 -10 5 10 -6 M 3 Á 10 -10 5 3 Á 10 -7 M 4.8 Á 10 -10 8 10 -7 M 6 Á 10 -10 10 3.5 Á 10 -8 M 1.3 Á 10 -9 21 concentration decreased to approximately 10 -5 M. The Cs flux at the outer boundary of the OPA specimen was measured. For our modeling the original Cs-isotherm parameters  were used, applying a CEC of 0.117 eq kg À1 for the PS as derived in Appelo et al. (2010). The sorption selectivity coefficients  were slightly adapted (±0.05 log K units) within the range of their uncertainty. Other model parameters (porosity, bulk dry density, tortuosity of water tracer) were taken from Appelo et al. (2010), and the surface mobilities as determined here. The surface diffusion model shows good agreement with the trend of the experimentally determined Cs flux (Fig. 9). However, the simulated arrival of the Cs front is somewhat delayed and the flux from the peak until day 1000 is slightly overestimated compared to the experimental data. To account for the delay either the pore diffusion coefficient can be increased or sorption can be decreased, but both modifications lead to a too high peak of the breakthrough curve. Appelo et al. (2010), when fitting a model including diffusion in the DL to the experimental data, observed exactly the same discrepancies. They explained the early arrival front compared to their simulation by the existence of dead-end pores, leading to physical kinetics of sorption. Also, they had to distribute the surface charges unequally between the open and the postulated dead-end pores. We did not include any dead-end pores in our model, because the early breakthrough could also result from experimental details not considered appropriately in the modeling. Irrespective of the dead-end pores, their fitted Cs diffusion coefficient is rather large, leading to a tortuosity much smaller than that for HTO. They interpreted this high diffusion coefficient (low tortuosity) as being the result of additional interlayer or surface diffusion of Cs (besides some ion pairing in the DL), a process not considered in their model. Having largely different tortuosities for Cs and HTO is unexpected, as (pore) diffusion of both occurs in the whole pore water. The overall good match between the experimental data and our surface diffusion model, using a tortuosity for Cs identical to that of HTO, is a proof that indeed this process is relevant for Cs. Furthermore, the test of the surface diffusion model with independently estimated surface mobilities was successful. There was no need to change surface mobilities compared to our calibration experiments, nor to adapt the Cs pore diffusion coefficient to any value that is difficult to defend (value of D p;Cs used here: 3.36 Á 10 -10 m 2 s À1 ). Therefore, the surface diffusion model appears to be capable to predict Cs diffusion in Opalinus Clay. The experiment was additionally modeled using the classical Fickian diffusion model. Here, the effective Cs diffusion coefficient had to be fitted to 3.75 Á 10 -10 m 2 s À1 , a value about 7 times larger than the effective diffusion coefficient of HTO. This again shows the inconsistency that arises when using the classical Fickian diffusion model for modeling Cs diffusion in Opalinus Clay.

SUMMARY AND CONCLUSIONS
A multi-site surface diffusion model was implemented in the reactive transport code Flotran. In-diffusion experiments for Cs-134 in Opalinus Clay were conducted at different stable Cs background concentrations. Diffusion properties of Cs-134 clearly depended on the stable Cs background concentration. A single set of three sitespecific surface mobilities could be estimated by fitting the surface diffusion model simultaneously to the data of the eight in-diffusion experiments. Parameters for the adsorption isotherm, only slightly adapted compared to literature values (within the given uncertainty range), led to a good match with the isotherm data points derived from our experiments. The results show that a multi-site surface diffusion model can describe the concentration-dependent Cs diffusion in Opalinus Clay consistently based on a single set of surface mobilities, while the Cs pore diffusion coefficient could be calculated from that in water and the assumed tortuosity of a water tracer. The observed concentration dependency of Cs diffusion cannot be explained by a classical diffusion model or by models including diffusion in the DL only. With a classical Fickian diffusion model, diffusion coefficients that depend on the Cs background concentration would be obtained. Sensitivities of the surface diffusion model to sorption parameters and surface mobilities were evaluated. Finally, the surface diffusion model with the estimated surface mobilities could successfully predict an other, independent Cs diffusion experiment in Opalinus Clay.
Even though the surface diffusion model was successfully applied to two data sets (the eight in-diffusion and the radial-diffusion experiments), it is clear that the presented surface diffusion model has limitations and has to be further improved. Site-specific surface mobilities of cations that compete with Cs for sorption in Opalinus Clay should probably be included for a more accurate description of Cs diffusion and sorption in situations with variable Fig. 9. Cs flux at the outer boundary of an Opalinus Clay specimen used in a radial diffusion experiment described in Appelo et al. (2010); Red line: surface diffusion model; Cyan dots with error bars: experimental data. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) pore water chemistry. Additionally, the application of molecular modeling of cation diffusion at clay surfaces could support the determination and/or validation of cation surface mobilities (e.g., Bourg and Sposito, 2011). Finally, at present the model focuses only on cations and does not consider effects of anion exclusion. Accordingly, a combination with a model that takes into account anion exclusion, such as DL models, is required. The combination of diffusion of cations that are specifically sorbed (e.g., in the Stern layer), diffusion in the DL and diffusion in free pore water could eventually lead to a consistent and more detailed model including all known diffusive processes in clays.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

ACKNOWLEDGMENTS
Financial support by the Swiss National Science Foundation (grant no. 200021_166287) is kindly acknowledged. We thank Dr. Bart Baeyens, Dr. Dan Miron and Dr. Wilfried Pfingsten for the fruitful discussions.

APPENDIX A
A.1. Calculation of the derivative term ›S=›C of the sorption isotherm A general cation exchange reaction between cations A and B onto a sorption site can be written as: where C A and C B are the concentrations in solution, z A and z B are their charges and S A and S B are the sorbed concentrations on a cation exchange site. The corresponding mass action law after the Gaines-Thomas convention (Gaines and Thomas, 1953) is: where k 1 K c is the selectivity coefficient of cation k with respect to the reference cation 1, z is the cation charge, N is the fractional occupancy on the exchanger solid, i.e., the moles charge of cation k on the exchanger per cation exchange capacity, and a is the cation activity. The activity is related to the concentration C via: where c is the activity coefficient. The sorbed concentration S k of cation k (in moles per mass of dry solid) is related to the fractional occupancy according to Rearranging Eq. (A1) for N k leads to an expression for the fractional occupancy N k of cation k: Thus, N k depends on the activity a k and charge z k of cation k in solution, on the selectivity coefficient k 1 K c , and on the fractional occupancy, the activity and the charge of the reference cation 1. The partial derivative @N k =@a k is then (applying the chain rule): The dependency of a 1 on a k or the partial derivative @a 1 =@a k , respectively, is most often very small and can be neglected. The combination of Eq. (A3) and Eq. (A5) results in the derivative of the sorption isotherm: The sum of the fractions of all competing cations for an exchange site is equal to 1 (full surface site coverage): Combining Eq. (A4) and Eq. (A7) we obtain an equation for N 1 as a function of all cation activities and selectivity coefficients: For homoionic exchange (z i =z 1 ¼ 1), Eq. (A8) can be solved analytically, but in general, N 1 is obtained numerically by using the secant method (or any other suitable method) to find the root. The partial derivative @N 1 =@a k is then obtained also numerically by varying a k in Eq. (A8). With N 1 and @N 1 =@a k , we can then derive the required derivatives according to Eqs. (A5) and (A6). The derivatives @S k =@C j appearing in the more sophisticated Eqs. (3) and (4) are evaluated numerically accordingly by varying a j instead of a k .