## 2 Methods
### Modelling Multi-Component Diffusion

Table 2: Porosity (e), dry and wet bulk density (\(\rho_{ab}\) and \(\rho_{wb}\)) as well as wet water content (\(wc_{wet}\)) used for the scenarios with the minimum and maximum amount of rock per kilogram porewater for the facies of the Opalinus Clay [10]. the remaining net surface charge and anions are excluded. Consequently, the pore space accessible for anions is reduced and so their diffusive flux. This is known as the anion exclusion effect. In the simulations, anion exclusion is controlled by the composition and amount of water in the DDL. The composition of the DDL is calculated as average using the fast and robust method of the Donnan approximation implemented in PHREEQC [13], which assumes a single potential for the DDL as an entity [5]:  
\[c_{i,DDL}=c_{i}\cdot exp\left(\frac{-z_{i}F\psi_{DDL}}{RT}\right) \tag{4}\]  
with \(c_{i,DDL}\) the concentration of species \(i\) in the Donnan pore space (mol/m\({}^{3}\)) and \(\psi_{DDL}\) the potential in the DDL (V) to ensure zero charge of the paired cells. The potential \(\psi_{DDL}\) of the Donnan volume is adapted to counterbalance the surface charge [13]. Hence, the concentrations in the DDL are linked to the concentrations in the free porewater \(c_{i}\) (mol/m\({}^{3}\)) and the surface charge that must be neutralized [39]. The thickness of the DDL is determined via the ionic strengths (Table 1) [36]. MC diffusive transport in a model element is calculated by explicit finite differences for each interface among the cells and summed up for the diffusive transport between two model elements. A more detailed description of the approach can be found in [PERSON] [5], [PERSON] [38] and [PERSON] et al. [6]. In porous media, diffusive transport is impeded by tortuosity of the pores, the reduced cross-sectional area available for diffusion and the pore sizes [40], which is reflected by the effective diffusion coefficient \(D_{\epsilon}\) (m\({}^{2}\)/s). PHREEQC assumes that a model element exclusively contains water [13] and therefore the pore water diffusion coefficient \(D_{p}\) (m\({}^{2}\)/s) is applied instead to calculate the diffusive flux. Equation (5) gives the relation between \(D_{p}\) and \(D_{\epsilon}\). In PHREEQC, \(D_{p}\) is derived analogous to [PERSON]'s law via the self-diffusion coefficient in water \(D_{w}\) (m\({}^{2}\)/s), the accessible porosity \(\epsilon\) (-) and an empirical exponent \(n\) (-):  
\[\frac{D_{\epsilon}}{\epsilon}=D_{p}=D_{w}\cdot\epsilon^{n} \tag{5}\]  
The empirical exponent \(n\) is a constant, medium-specific parameter [41] and needs to be determined for the Opalinus Clay (Section 3) to use the MC approach in PHREEQC. In  
Figure 1: Concept of the one-dimensional finite-difference diffusion model in PHREEQC. order to individually calculate the diffusive flux for each species present in the system, \(D_{w}\) must be assigned in the database. The \(D_{w}\) of all relevant uranyl species were calculated by molecular-dynamic simulations and are taken from [PERSON] and [PERSON] (2014) as well as from [PERSON] et al. (2015). In the case of a missing \(D_{w}\) for an uranium species, the diffusion coefficient of a chemically similar one was chosen to be used instead. All other necessary self-diffusion coefficients stem from the database phreeqc.dat provided with PHREEQC ([PERSON], 2013). The selected \(D_{w}\) are given in the Supplementary Materials S1. On the metre-scale, MC simulations were performed for each facies (Table 1, Section 2.1) as well as for varying mineral quantities and compositions (Table 2, Section 2.2) for models of 1 m with a spatial resolution of 0.01 m and a simulation time of 100 years. Simulations on the host rock scale were conducted for each facies with a model of 50 m and a spatial resolution of 0.5 m and a simulation time of one million years. Grid independence of the diffusion models on both scales was confirmed by simulations with finer resolutions.