Modelling the kinetic reactive transport of pollutants at the sediment-water interface. Applications with atmospheric fallout radionuclides

Modelling the kinetic reactive transport of pollutants at the sedimentwater interface. Applications with atmospheric fallout radionuclides J.M. Abril a , * , H. Barros b a Departamento de Física Aplicada I, ETSIA, Universidad de Sevilla, Carretera de Utrera Km 1, 41013, Sevilla, Spain b Universidad Sim´on Bolívar, Departamento de Física, Apartado, Postal 89000, Caracas, 1080, Venezuela Keywords: Surficial sediments Trace elements Uptake kinetics Reactive transport Early compaction Eddy diffusivity ABSTRACT Understanding the behaviour of particulate matter and chemicals at the sediment-water interface (SWI) is of interest in environmental studies and risk assessments. These processes are still poorly understood, and this work aims to gain relevant insights by using a kinetic reactive transport model. It merges early diagenetic processes and box models for the uptake kinetics. Numerical solutions have been found for synthetic scenarios and for studying real cases from the literature ( 210 Pb and Chernobyl fallout radionuclides in Lake Sniardwy, Poland, and 7 Be in sediments from Tema Harbour, Ghana). The study identifies a series of factors that dynamically interact to govern the final fate of tracers in the SWI region, leading to a wide diversity of behaviours. When a term of eddy diffusivity is included in the upper regions of the pore fluid, which seems feasible for some energetic scenarios, it is possible to explain the observed large penetration depths for Cs and Be, while high particle-reactive elements are retained in thinner sediment layers. Desorption from the sediment occurs through the pore fluid as diffusive fluxes. Transient depth profiles of tracer concentrations can last from months up to a year, and they can show subsurface maxima at positions unrelated with the accretion rate. In the application cases, the model explained a wide set of observational data that was beyond the capabilities of other approaches involving physical mixing of solids and equilibrium k d . This modelling study could provide useful guidance for future research works. 1. Introduction Anthropogenic activities are increasingly altering the landscape and affecting the fluxes of pollutants and sediments that, in turn, result in the degradation of aquatic environments. At the sediment-water interface (SWI) meet the fluxes of matter, pollutants, and trace elements from diverse provenances, including atmospheric deposition. Understanding the behaviour of particulate matter and chemicals in the sediment region below the SWI is of great scientific interest. Among other topics, it is essential for assessing the fate of pollutants, including their potential remobilization and bioavailability, and for providing suitable boundary conditions for radiometric dating models for recent sediments (Mabit et al., 2014). These latter can then support a large variety of interdisciplinary studies. Examples include, among others, decoding past anthropogenic impacts and pollution history (Kennish, 2001; Bellucci et al., 2007; Cundy and Croudace, 2017), quantifying the carbon sequestration capacity of vegetated coastal sediments (Marchand, 2017), or assessing sea level changes at centennial scale (Cundy and Croudace, 1996; Barnett et al., 2017). Particulate matter in natural aqueous suspensions and surficial sediments show in their surfaces uncompensated electrical charges (mostly negative) so they can uptake ions of opposite charge from the dissolved phase (Eisma, 1993). Along with electrostatic adsorption, other physicochemical processes contribute to further fixation of trace elements at different time scales (Duursma and Carroll, 1996). As a result, for some trace elements the distribution or partition coefficient, k d (defined as the dimensionless ratio at equilibrium between the mass concentrations in solids and in solution), is of the order of 10 3 to 10 7 (IAEA, 2004). Consequently, although sediments are saturated porous media, for most trace elements, almost 100% of them are expected to be bound to the solid phase. This justifies the common approach of considering the sediment as a continuous medium, allowing for the following continuity equation (Abril, 2003): ∂ ( ρ A) ∂ t= ∂ ∂ z(D ρ∂ A ∂ z)− ∂ ∂ z(wA)(1) Abbreviations: SWI, sediment-water interface; NID, non-ideal deposition. * Corresponding author. E-mail address: [email protected] (J.M. Abril). 2 In Eq. (1) A is the mass concentration of the studied element (with physical dimensions of MM −1 ), ρ is the bulk density [ML −3 ], w [ML −2 T −1 ] is the advective mass flow due to sediment accretion (see more details below), D [L 2 T −1 ] is a diffusion coefficient, z [L] is the depth coordinate, measured positive downwards the SWI, and t [T] is time. In the case of unsupported radionuclides, A is the activity concentration (in units of Bq/kg or similar), and the term −λ ρ A must be added to the right side of Eq. (1), where λ [T −1 ] is the radioactive decay constant. Solving Eq. (1) requires providing initial and boundary conditions. At the SWI (z =0) the continuity of fluxes is commonly imposed: F(t)=− D ρ∂ A ∂ zz=0+wA(0,t)(2) where F(t)[ML −2 T −1 ] is the flux of the studied element at the SWI. The continuity equation for solids is written in terms of a compaction potential ψ , which takes into account that the potential gravitational energy of the system decreases when the water pores are progressively replaced by the overlaying solids (Abril, 2003, 2011): ∂ρ ∂ t= − ∂ w ∂ z;w= ρ (v+q); ρ q= − K(z) ∂ψ ∂ z(3) In Eq. (3) v is the velocity at which the SWI displaces up due to sediment accretion, and q is the mean velocity of the solids due to early compaction. ψ [L] has units of energy per unit weight, and K(z) [ML −2 T −1 ] can be interpreted as a conductivity function. The natural boundary condition for Eq. (3) in the SWI is the continuity of the mass flow, w0. The fundamental differences of this formulation with respect to the classical diagenetic equations of Berner (1980) are discussed in Abril (2003). Equation (1) is often simplified by assuming steady state bulk density ( ∂ρ ∂ t=0; and thus, w(z) = w0), which allows rewriting it in terms of a mass depth scale (Abril, 2003). All the radiometric dating models for recent sediments are particular solutions of this simplified Eq. (1), with the consistently simplified boundary condition of Eq. (2). The most widely used 210 Pb-based models assume D=0 (see an updated summary of these models in Iurian et al., 2021). The above approach has been applied in a huge number of scientific studies to date. However, concerns about the potential mobilization of pollutants through the pore water in marine and lacustrine sediments also have a long track. Scholkovitz and Mann (1984) first reported empirical studies on the pore water chemistry of 239,240 Pu and 137 Cs in marine sediments. Comans et al. (1989) reported in situ k d values (the ratio of tracer concentrations in solids and pore water) for 137 Cs in lacustrine sediments affected by the Chernobyl accident, providing convincing evidence for post-depositional mobilization of particle-bound 137 Cs. Benoit and Hemond (1991) carried out similar studies reporting evidence for diffusive redistribution of 210 Pb in lacustrine sediments. Although being elegant and convenient from a mathematical viewpoint, in some cases Eq. (2) may not be a reliable description of the microcosms in the SWI region. Thus, in some lacustrine and marine sediment cores sampled shortly after the Chernobyl accident in areas affected by its radioactive fallout, it was observed that 134 Cs had penetrated into the sediment up to 10–20 cm depth (e.g., Petersen et al., 1990; Kramer et al., 1991; Robbins and Jasinski, 1995; Holby and Evans, 1996). Similar observations have been documented in some sedimentary systems affected by the Fukushima accident (e.g., Yoshimura et al., 2014; Otosaka et al., 2020). Diffusion models provided apparently satisfactory descriptions of the data, but in some cases, the required diffusion coefficients had unreliable values (Holby and Evans, 1996). Furthermore, when applying these diffusion models to a sequence of cores taken on different dates, it was necessary to use values for the mixing coefficient that decreased with time (Robbins and Jasinski, 1995). Abril and Gharbi (2012) noticed that these observations represented a failure of the boundary condition of Eq. (2). Reasons argued were a fast initial redistribution in depth of the radioactive fluxes when an important fraction of them was associated to the dissolved phase, which was referred to as non-ideal deposition (NID). They used a simple numerical model for the uptake kinetics in the porous media for concept demonstration. They then suggested a consistent adaptation of Eqs. (1) and (2) by introducing a depth-distributed source term and a practical formulae for its parameterization. This allowed a successful description of some empirical profiles of 137 Cs and unsupported 210 Pb. Although being a step in the right direction, a series of concerns remain regarding the conditions for its applicability to a wide variety of tracers and sedimentary conditions, and on the involved time scales and depth-distribution patterns. Laboratory experiments studying the uptake of radiotracers by sediments from overlying water columns have shown quite limited penetration depths after several months of incubation (e.g., ~4 cm after 1 year of incubation for 134 Cs with lacustrine sediments, after Smith et al., 2000; and ~2 cm after 221 days for 133 Ba in experiments with estuarine sediments, after Barros et al., 2003, 2004). These results have been explained in terms of molecular diffusion through interstitial water and fast equilibrium in the partitioning, although inconsistencies due to the required time-dependent k d have been highlighted, as in the cited references. Furthermore, these penetration depths are far from the 10–20 cm range above commented. A good number of open questions remain unsolved, particularly identifying those factors linked to both the sedimentary system and the studied tracer elements that govern NID. How do different tracers behave under the same sediment environment? Can the behaviour of a given tracer change with depth or with time? What is the effect of varying porosity? How do percolation flows enhance the penetration depth of the tracers? What is the role of eddy diffusion? What is the effect of dynamic changes in tracer concentration in overlaying waters? How are the transient distributions and the involved time scales? When can be expected to occur pollutant desorption through the pore water? How generalized can NID be to represent a serious correction to existing radiometric dating models? The above questions can be faced by using reactive transport models, which have been a fruitful research approach for the Earth sciences during the last four decades. They combine the transport equations for one or multiple chemical species moving through a geophysical fluid with the physic-chemical reactions that they undergo; these last have been described with equilibrium and kinetics approaches (Steefel et al., 2005; Berkowitz et al., 2016). Examples of applications are, among others, studies on radioactive waste disposal (e.g., Glassley et al., 2003), reactive transport in fractured rocks (e.g., MacQuarrie and Mayer, 2005), dispersion of reactive pollutants in the marine environment (e.g., Peri´ a˜ nez, 2005), and geochemical reactions in petroleum subsurface reservoirs (e.g., Abd and Abushaikha, 2021). This paper aims to obtain relevant insights into the processes that governs the dynamic of pollutants at the sediment-water interface by using a simple kinetic reactive transport model specifically adapted to this environment. 2. Methods This section comprises a physical description of the SWI region, the mathematical approach used for the uptake kinetics, and the governing equations for the kinetic reactive transport in the porous media, along with the implementation of suitable boundary conditions for the mathematical problem. 2.1. The physical environment at the SWI region Fig. 1 shows a sketch to support discussion. Stacked solids conform the structure of the porous medium, with the pore spaces occupied by water. The porosity, φ (the fraction of the bulk volume occupied by 3 pores), can be higher than 90% in some cases. Solids are of varying composition, shapes and sizes, and organic matter may be present (Boudreau, 2000). The bulk density relates with porosity, so that ρ = ρ s(1−φ), where ρ s is the density of solids. Because of the stacking, a fraction of the free surface of solids remains inaccessible to pore water. Solid particles have rugose and porous surfaces where uncompensated electric charges and more specific reaction sites are distributed (Eisma, 1993). The uptake of trace elements from the dissolved phase is a kinetic process that can be described by a set of parallel and consecutive reversible reactions, along with some irreversible components (Barros and Abril, 2008). The SWI is not a mathematical planar surface, but an irregular, rugose, and porous interface. The trace elements reaching the SWI can be associated with a large diversity of carriers. A fraction of them can be free ions in solution, while other fractions can be bound to solid particles of different compositions and sizes. Tracers in solids may have reversible and irreversible bound components. Due to the previous history, these composite fluxes are tracer-specific and may show short- and long-term temporal variability (Abril and Brunskill, 2014). The sediment structure is dynamically rearranged in a gravity field, with solids replacing underlying water pores. This process is known as early compaction and makes porosity asymptotically decrease over a certain depth interval in the sediment to reach a virtually uniform value (Abril, 2003). Early compaction is driven by overbounding pressure and mechanical perturbations such as pressure waves, micro-earthquakes, etc. It can progress as a uniform decrease in porosity without perturbing the isochronous horizons or with a differential translocation of particles of different sizes. In dynamical sedimentary scenarios the friction stress of water currents at the SWI can lead to sediment resuspension, with episodic or periodic (e.g. in tidal systems) destruction and re-composition of the upper SWI region. The threshold friction stress for resuspension is a function of particle size and the roughness length. Typical values of 0.06 Pa are required for the clay fraction, which roughly corresponds to depth-averaged water currents of 5 cm s −1 (Eisma, 1993). Higher values are required for larger particle sizes. The diffusion term in Eq. (1) accounts for sub-grid scale advection. This is, displacements of solids up and down a horizontal planar surface, in such a way that there is no net mass transport, but the concentration of the tracer changes when there are spatial gradients. These null net mass flows do not affect the bulk density, so diffusion is absent in Eq. (3). The physical reworking of solids against the gravity field requires a forcing agent, such as bioturbation. Differential translocation of solids by early compaction can represent a tiny contribution (of the order of 10 −9 cm 2 h −1 after a Fermi estimate, based upon a Reynolds approach). The resuspension and deposition cycles could support in some cases a mathematical description as equivalent diffusion. All the above processes are clearly depth dependent. The chemical environment progressively changes with depth towards anoxic conditions and with varying concentrations of competing ions, such as NH 4+ (Comans et al., 1989). 2.2. Uptake kinetics in a porous media The scientific literature reports numerous results from experimental and modelling studies on the uptake of radionuclides and other trace elements by natural aqueous suspensions. Some examples are shown in Fig. 2, and they are briefly explained in the figure caption. The time series of tracer concentrations in the dissolved phase can be mathematically fitted to an analytical function involving a set of exponentials Fig. 1. Sketch illustrating the physical environment of the SWI region and the composite nature of the mass flows of solids and trace elements (see text). Dots represent ions, in a dynamical exchange among dissolved and solid-bound phases, and different geometric figures entering the SWI represent different carriers conforming the composite flux of the tracer. The artwork uses a template from the Pore Network Comparison Forum (ANU) – the structure of the porous media (left), and the subfigure of the granular medium (Boudreau, 2000) –(right up). 4 that are associated with different characteristic times in the uptake, ranging from minutes up to weeks (Barros and Abril, 2005). This function arises as the mathematical solution of box models involving the dissolved phase and up to three different reaction sites within the particle. They are connected by parallel or consecutive reversible reactions, or by a mixture of them, which are described by constant kinetic coefficients. Some examples of box models are shown in Fig. S1, in the electronic supplementary material (ESM). In some cases, models include some irreversible reactions and saturation terms (e.g., El-Mrabet et al., 2001). It has been shown that the above empirical data do not allow unambiguously identifying the involved processes in terms of parallel or consecutive reactions, since it is always possible to find the same analytical solution using any possible variation of the reactions schemes (Barros and Abril, 2004, 2005). The factors affecting the kinetic coefficients and the partitioning at equilibrium have been profusely studied from both experimental and theoretical approaches. Abril (1998a; 1998b) presented a basic microscopic theory for the uptake kinetics of trace elements by suspended particulate matter. The direct transfer coefficient is proportional to the free surface of solids per unit volume of solution and to a geometric probability of reaction. Thus, its value depends on the experimental setup (e.g., on the volumetric concentration and the granulometric and mineralogical characteristics of suspended solids). The more physically significant proportionality factor could be isolated, but there is no unified approach to estimate the specific surface area (SSA) of solids. Some works use the so-called external SSA, estimated from the size spectrum and the assumption of spherical shapes. Other works report SSA measured by gas adsorption using the BET isotherm, which accounts for non-spherical shapes, and the rugose and porous surfaces, leading to much larger values than the above method (see examples in the references cited in the caption of Fig. 2). The geometric probability of reaction depends on the surficial concentration of reaction sites but also on the concentration of competing ions into solution, which explains the dependence of k d with pH, salinity, and NH 4+ (Abril and Fraga, 1996; Comans et al., 1997). The coefficients for reverse reactions are related to the probability of rupture of the unions, which can be mediated by the composition of the aqueous solution and the temperature, but they do not depend on SSA. As the uptake is a genuine surface phenomenon, the inner portion of the solids are prevented for any sorption or desorption, which leads to the well-known particle size effects in k d (Abril and Fraga, 1996). Because of the mathematical equivalence of models, in this work we have adopted for simplicity a scheme of consecutive reactions (case A in Fig. S1, in ESM) that for two components in the solids can be written as follows: daw dt = − k11aw+k21as,1 das,1 dt =k11aw− (k21 +k12)as,1+k22as,2 das,2 dt =k12as,1−k22as,2 (4) In Eq. (4), aw is the volumetric concentration of the tracer in dissolved phase, while as,1 and as,2 are their volumetric concentrations in reaction sites 1 and 2 in the solids; k1i and k2i (i =1,2) are the kinetic coefficients for direct and reverse reactions, with units of inverse time. It is worth noting that the radioactive decay constant for radionuclides does not need to be explicitly included in Eq. (4) when half-lives are much longer than the duration of the experiment. More complex Fig. 2. Experimental data and fitted curves with box models (see Fig. S1, in ESM) for a set of adsorption experiments taken from the scientific literature. The tracer is added to the solution and the remaining fraction in this phase is plotted versus the contact time. Am (I) and Am (II) correspond to the uptake of 241 Am by aqueous suspensions from the Guadalquivir River (SW Spain) of 43 and 172 ppm with external SSA of 1.33 m 2 g -1 (experimental data from El-Mrabet et al., 2004, and box models from Barros and Abril, 2005). Ba corresponds to the uptake of 133 Ba by suspended estuarine sediments of 20 g l −1 with an external SSA of 0.010 m 2 g -1 (experimental data from Barros and Abril, 2004). Co corresponds to the uptake of 60 Co by marine sediments (2 g underlying 20 mL of seawater) with a SSA by the BET method of 13 m 2 g -1 (experimental data from Børretzen and Salbu, 2000, and box model from Barros and Abril, 2008). Pu corresponds to the uptake of 239 Pu by natural aqueous suspension from the Guadalquivir River of 12.6 ppm with external SSA of 1.63 m 2 g -1 (experimental data and box model from El-Mrabet et al., 2001). Cs corresponds to the uptake of 134 Cs by marine sediments (2 g underlying 20 mL of seawater) with a SSA by the BET method of 13 m 2 g -1 (experimental data from Børretzen and Salbu, 2002). Pb corresponds to the uptake of 210 Pb by suspended marine sediments of 52.6 ppm and size <63 μ m (experimental data from Wang et al., 2019). 5 kd=k11(k22 +k12) k21k22 ms (5) where ms is the concentration of suspended solids (kg/L) and the density of the water (1.0 kg/l) is implicitly included, without corrections in volume because ms, which are assumed to be small. Equation (4) can be locally adapted to the porous sediment medium, assuming homogeneity of concentrations in a region characterized by a porosity φ, a mean density of solids, ρ s, and in absence of advectivediffusive transport. The volumetric concentrations in Eq. (4) can be expressed in terms of their corresponding mass concentrations (denoted with upper ˆ): aw= ρ wφaw;as,i= ρ s(1−φ)as,i ( ρ w is the density of the pore fluid). Solving for k d : kd=k11(k22 +k12) k21k22 ρ wφ ρ s(1−φ)(6) From the above discussion, the state-of-the-art does not allow a priori reliable estimates of the kinetic coefficients for the trace elements of interest. They will take site-specific values that capture the physicchemical properties of the solid and aqueous phases, and they will vary with space and time as such properties do. Therefore, transient stages are expected to be ubiquitous. The previous presentation focused on sorption experiments. There is wide experimental evidence that after long contact times, an important fraction of the tracer becomes irreversibly bound to the solids. As an example of these studies, Børretzen and Salbu (2000) conducted desorption experiments with Cd and Co in marine sediments. They used a sequence of reagents having increasing displacement and dissolution power. Seawater only could desorb about 10% of Cd and a negligible amount of Co. The selected approach for this study uses a set of synthetic kinetic coefficients that define several classes of trace elements with increasing affinity for solids (Table 1). These values are inspired in experimental data (e.g., classes A and B are adaptations from experiments with 133 Ba reported in Fig. 2, while class D uses results from 241 Am). Table 2 summarises the recommended k d values for a set of elements in various aquatic environments. In each case of interest, the behaviour of the target element can be considered to be close to one of the classes in Table 1. Increasing k d values are often associated with a faster uptake kinetics, as shown in Fig. S2 (in ESM) and Table 1. Different equilibrium times are possible for a given k d , as illustrated in Fig. S2 and Table 1 with tracer A*. This last uses the set of kinetic coefficients of tracer A, but multiplying k22 and k12 by a factor (1 4), what leaves invariant the k d (Eq. (6)). The effects of irreversibility and saturation will be incorporated some further as a modification of data in Table 1. It is worth noting that for all the tracers in Table 1, after equilibrium more than 99.8% of the tracer was associated to solids. Within solids, sites 2 stored approximately 98% of tracer A and 99.99% of tracer D. 2.3. Constituent equations for the kinetic reactive transport of trace elements in the SWI region The conceptual model tries to capture the essential processes described above while keeping complexity at a minimum, but with flexibility for necessary extensions. The first point is to solve the velocity field for the physical fluid carrying the tracers studied. On the time scales of interest, both the aqueous and solid phases of the sediment can undergo transport processes. A steady state will be assumed with a mass flow wo, and percolation velocities, vp,0, impossed as boundary conditions at the SWI. Then, Eq. (3) allows estimating the accretion and compaction velocities from the empirical bulk density profile. The reactive transport model will be formulated for a single component that is presented in dissolved form and in two reaction sites in solids, connected by consecutive reversible reactions. Some flexibility will be necessary to include irreversibility in the uptake. The transport equations for the tracer use the velocity field of the physical medium and incorporate a diffusion term in the fluid phase, along with the terms describing the kinetics of the uptake. The kinetic coefficients in Table 1 are assumed to account for the effective surface area of solids accessible to the dissolved phase. A first goal is studying the conditions that can lead to a rapid initial redistribution over large depths of trace elements that reach the SWI into the dissolved phase. When the involved time scales are of the order of few days to few weeks, sediment accretion can be neglected, and radioactive decay is negligible in most cases. The equations for the reactive transport can then be written as follows: d( ρ wφaw) dt =− k11 ρ wφaw+k21 as,1 ρ s(1−φ)+ ∂ ∂ z(Deff ρ wφ ∂ aw ∂ z)− ∂ ∂ z(vp ρ wφaw) d dt(as,1 ρ s(1−φ))=k11 ρ wφaw−[(k21 +k12)as,1+k22 as,2] ρ s(1−φ) d dt(as,2 ρ s(1−φ))=[k12 as,1−k22 as,2] ρ s(1−φ) (7) In Eq. (7), vp is the percolation velocity, and Deff is an effective diffusion coefficient. This last can be written in terms of the molecular diffusion, Dm, corrected by the tortuosity of the porous media (the factor φn, with n in the range 1–2 is a practical empirical correction – see Krezoski et al., 1984; Barros et al., 2004), and it may include a term for eddy diffusivity De, as discussed further below: Deff =De+φnDm(8) At the SWI (z =0), a flux of tracers, Fw, is entering into the dissolved phase, fulfilling the boundary condition Fw(t)=− Deff ρ wφ ∂ aw ∂ zz=0+vp,0 ρ wφaw(0,t)(9) Any analytical function or numerical time series can be imposed as Fw(t). Alternatively, it can be estimated as the advective-diffusive transport with an overlaying water column where the concentration of dissolved tracers is known as a function of time, a0,w(t). This work will show that, in some cases, transient stages persist from several months up to more than one year. Thus, sediment accretion and radioactive decay must be explicitly included in the equations. The last is a trivial task. For the former, a simplified approach will assume a constant mass sedimentation rate, w, and a steady state bulk density profile with a known early compaction limit with bulk density ρ ∞. The SWI displaces up, due to sedimentation, with a linear velocity v=w ρ ∞. From a reference frame anchored to this SWI, the velocity of solids, vs, and pore water, vw, at depth z are, respectively (Abril 2003, 2011): vs(z)=v+q=w ρ (z)(10a) Table 1 Kinetic transfer coefficients defining various categories of trace elements used for the numerical modelling of reactive transport in the SWI region. Class k11 (h −1 ) k21 (h −1 ) k12 (h −1 ) k22 (h −1 ) kdteq (h) A 12.5 1.0 0.25 0.005 1.0⋅10 3 50.5 A* 12.5 1.0 0.0625 0.00125 1.0⋅10 3 120 B 25.0 2.0 0.5 0.001 1.0⋅10 4 31.0 C 200 6.0 1.0 0.0005 1.1⋅10 5 15.0 D 500 10.0 1.3 0.0001 1.0⋅10 6 11.0 kd values estimated from Eq. (6) with φ=0.8, ρ s=2.5 g cm −3 and ρ w=1.0 g cm −3 . Equilibrium times from Fig. S2 (ESM). reaction schemes are possible, but Eq. (4) accounts for the major features of the uptake kinetics and is versatile enough for the present goals. For the simplest case of constant coefficients, the partition coefficient at equilibrium (the time derivatives in Eq. (4) become null) is: 6 vw(z)=vp(z)+v−q′;q′=q1−φ(z) φ(z)(10b) In Eqs. (7) and (9) vp must be replaced by vw(z), and an advective term must be added to Eq. (7) for as,i (i =1,2): − ∂ ∂ z[vs ρ s(1−φ)as,i](11) At the SWI the flux of tracer associated to the reaction sites in solids must be specified: Fs,i(t)=vs(0) ρ s[1−φ(0)]as,i(0,t)(12) The system of Eqs. (7)–(12) conforms a model being general enough as to account for composite fluxes with dissolved and particle-bound fractions. Eq. (12) allows one to impose any degree of transient partitioning among phases. The occurrence of an irreversibly bound fraction in the mass flow reaching the SWI can be included as a third reaction site in solids, as,3, that only undergoes accretion (Eq. (11)) and radioactive decay. This third reaction site can be connected to sites 2 as a possible formulation to model the progression of irreversibility. For numerical solutions, we used an upwind finite-difference scheme, corrected by numerical dispersion when relevant. The spatial discretization, Δz, and the time step, Δt, are selected to meet the criteria for numerical stability (Peri´ a˜ nez, 2005). The lower boundary is positioned deep enough so that the used boundary conditions (continuity of gradients) have negligible effects on the solution. The quality of the solution is checked by monitoring the mass balance of the tracer over time. Nyffeler et al. (1986) used a similar modelling approach including basic processes in the overlaying water column and a diffusion term for solids, but in most of the applications the contribution of as,2 was neglected. Robbins (1986) presented a two-phase model for radionuclide transport within sediments assuming instantaneous equilibrium between the dissolved and particle bound phases. The same approach appears in the work by Smith et al. (2000). The previous work by Smith and Comans (1996) included sediment accretion and physical diffusion affecting both aqueous and solid phases, and a k d value scaled by the depth-dependent concentration of NH 4+ . More recently, a multi-fractional sediment transport module similar to Eq. (7) has been incorporated into a more general 3D radioactivity transport model coupled with multiscale circulation (Maderich et al., 2017). 3. Results and discussion 3.1. Contact time The following setup has been selected for studying the effect of increasing contact times in the absence of sediment accretion: Uniform value of porosity, φ=0.8, kinetic coefficients for classes A and C (Table 1), vp=0; ρ w =1 g cm −3 ; ρ s =2.5 g cm −3 ; molecular diffusion (Eq. (8)) with Dm=0.05 cm 2 h -1 (Smith et al., 2000) and n =1.75. As the boundary condition (Eq. (9)), Fw(t)was analytically formulated as a Gaussian pulse that peaks at t =15 h and with σ =4 h. The total amount of tracer injected was 100, in arbitrary units (a.u.). The simulation times were set as 5 d, 10 d, and 15 d. The core length was 10 cm for tracer A and 5 cm for tracer C. The results are shown in Fig. S3 (in ESM). The aggregated amount of tracers per unit dry mass (a.u. g −1 ), A(z), has been selected as the meaningful magnitude, since it would correspond to the empirical value measured in a real sediment core after slicing and drying. The kinetic reactive transport driven by molecular diffusion allows a depth distribution of the Gaussian pulse of tracers. After 5 days, the remaining amount of tracer in water is unable to promote any further redistribution for tracer C, with a higher nominal k d value (Table 1), while some tinny evolution can be observed for tracer A (Fig. S3, ESM). The penetration depth is defined here as the depth above which 99% of the inventory is located. It is 0.18 cm for tracer C and 0.70 cm for tracer A after 10 d. The in situ k d , locally estimated from the ratio of mass concentrations of the tracer in solids and pore water, changed with depth and time (Fig. S3, ESM). Unlike in the absence of transport processes (Eq. (6)) the tracer concentration in dissolved phase changes with time because uptake by solids, diffusion, and the time-dependent fluxes at the SWI. As a result, k d reflects transient stages. A similar decrease in in situ k d with depth has been described for lacustrine sediments by Benoit and Hemond (1991). 3.2. Porosity and percolation The numerical setup is as the one described in the previous section, but using different values for the porosity. The simulation time was set to 10 d. In all cases, the depth distribution of A(z)was similar to that shown in Fig. S3 (ESM). For the sake of brevity, only the penetration depth is reported in Table 3. Porosity, through the parameterization of tortuosity, affects the magnitude of the effective diffusion coefficient (Eq. (8)). A second effect is that higher porosities enhance the role of the dissolved phase. In the absence of percolation, increasing porosities promote slightly larger penetration depths for tracer A, but their effects become negligible for tracers with higher nominal k d values. Even with very high porosities, the solid phase acts as an effective filter trapping the high particlereactive tracers in a few millimetres. Percolation velocities through sediment to the groundwater table usually take low values. In their study, Benoit and Hemond (1991) reported measured values of 0.6 mm/day. Here we explored the effect of a higher value of 3.0 mm/day, and the extreme case of 3 cm/day. The effect in the penetration depth is negligible, but for tracer A, with increments of 1.3% and 27%, respectively (Table 3). This happens again because the solids act as an efficient filter, retaining the particle-reactive tracers. The performance of the numerical scheme has been tested for a Environment Co Ce Cs Ba Pb U Pu Am Open ocean [1] 5⋅10 7 7⋅10 7 2⋅10 3 9⋅10 3 1⋅10 7 5⋅10 2 1⋅10 5 2⋅10 6 Ocean margin [1] 3⋅10 5 3⋅10 6 4⋅10 3 2⋅10 3 1⋅10 5 1⋅10 3 1⋅10 5 2⋅10 6 Freshwater [2] 4.3⋅10 4 2.2⋅10 5 9.5⋅10 3 2⋅10 3 1.1⋅10 5 - - 7.9⋅10 4 2.1⋅10 5 [1] Data from IAEA (2004); [2] data from IAEA (2010) but for Pb, which is reported from Mudbidre et al. (2014). Table 3 Penetration depth (given in cm) of different tracers for varying uniform porosities and percolation velocities (vp, in mm/day). Results from numerical simulations with the model setup of section 3.1 Tracer vp=0 φ=0.90 φ=0.65 φ=0.80 φ=0.90 vp=3 vp=30 A 0.59 0.70 0.78 0.79 0.99 B – 0.43 0.47 0.48 – C – 0.18 0.19 0.19 0.20 Gaussian pulse for Fw(t);Dm=0.05 cm 2 h -1 and n =1.75. Simulation time 10 days. Table 2 Recommended kd values for some trace elements in aquatic environments. 7 De=De,0e−βz (13) The following simulation setup is as in Section 3.1, with vp=0, φ= 0.90, and a simulation time of 10 days. In addition to molecular diffusion, two formulations of eddy diffusivity have been studied: E1 uses De,0=3 cm 2 h −1 (8.3⋅10 −8 m 2 s −1 ) with β=2.0 cm −1 and E2 uses De,0= 6 cm 2 h −1 with β=1.0 cm −1 . The results for tracers A and C are shown in Fig. S4 (ESM). Tiny eddy diffusivity in the SWI region has a major effect on the penetration depth of particle-reactive tracers, including those with high k d values. 3.4. Boundary conditions with a constant tracer concentration in the overlying water column In the numerical scheme, Fw(t)in Eq. (9) is now estimated as advective-diffusive transport with an overlaying water column where the concentration of dissolved tracers, a0,w, remains constant over time. The simulation setup uses vp=0, φ=0.80, a0,w=100 a.u.g −1 , De=0, and varying contact times of 1, 5, 10, and 20 days. The results are summarized in Fig. S5 (ESM). Total inventory and the penetration depth increased with time (Fig. S5, ESM). Analytical solutions assuming instantaneous partitioning can be obtained as follows: In Eq. (7), the second reaction is neglected (as,2=0, and as,1 is denoted as as), vp=0, porosity and Deff are assumed to be constant in time and uniform with depth. Summing the two equations and using that kd=as aw lead to: ∂ as ∂ t=Da ∂ 2as ∂ z2;Da=Deff 1+kd ρ s(1−φ) ρ wφ (14) For initially clean sediments with a boundary condition at the SWI of constant a0,w, the analytical solution is (Krezoski et al., 1984): as(z,t)=a0,wkderfc(z 2 Dat √)(15) In Eq. (15) erfc is the complementary error function. Taking into account the typical large values of kd, A(z) ≅ as(z,t). Fig. S5 (ESM) plots the analytical solution for tracer A at t s =20 d with a fitted value of k d = 750. The agreement seems quite reasonable, but doing the same for other contact times requires different fitted values for k d , namely 75, 300 and 500 for 1, 5 and 10 days, respectively. Similar results are found for tracer C, with fitted k d values of 400, 1800 and 3300 for contact times of 1, 5 and 10 days, respectively. The failure of Eq. (15) because the needed time-dependent k d has been previously reported (e.g., Barros et al., 2003). The uptake kinetics is not an instantaneous process, and sharp gradients appear in the liquid phase at the SWI because solids continuously adsorb the tracer. The in situ k d estimated by the model reflects transient stages, which follow a trend of decrease with depth similar to that shown in Fig. S3 (ESM). The calculated penetration depths for tracers A and C, at t s =10 d, were very close to the values reported in Table 3 (for φ=0.80). A variant of this model setup is found in the laboratory experiments reported by Barros et al. (2003, 2004), where the concentration in the overlaying water column varied with time but asymptotically tended towards a constant value. With water at rest and at constant room temperature, eddy diffusivity is absent, and after several months of contact time, the penetration of tracers of type A was of a few cm, as previously commented. The effects of percolation on inventories and penetration depths have been assessed for tracer A, as summarized in Table 4. These magnitudes are enhanced only by a ~3% for vp =3 mm d −1 , a representative value for environmental conditions. The results shown in Fig. S5 (ESM) are in arbitrary units for the tracer concentration, but the use of constant kinetic coefficients in Eq. (7) implicitly assumes that such concentrations are far from producing any significant saturation of the reaction sites in solids. In case of planned or accidental releases resulting in high concentrations of the tracer, the saturation effect can promote larger penetration depths. A saturation of the second reaction sites can be modelled by using a concentrationdependent kinetic coefficient for the second direct reaction: Table 4 Inventories and penetration depths computed for tracer A after t s =20 d using the model setup of Fig. S5 (ESM), but with varying percolation velocities. vp=0 vp=3 mm/d vp=30 mm/d Inventory (a.u. cm −2 ) 6222 6456 8866 Penetration depth (cm) 0.74 0.76 0.93 conservative tracer (what can be achieved setting k 11 = 0 in Eq. (7)). In this case, the centroid of the pulse of concentration in the pore fluid displaced down with the prescribed percolation velocity. The tiny role of vp when compared against molecular diffusion can be understood by handling the concepts of mixing length and mixing velocity. Thus, in the numerical scheme, with spatial discretization Δz, diffusion through a planar surface can be seen as the result of two opposite mass flows with a mixing velocity vD ∼ φnDm/Δz (Abril and García-Le ´on, 1991). For Δz = 0.1 cm, vD ∼ 96 mm d − 1 , which is much higher than vp. It is worth noting that with this boundary condition, where Fw(t) is prescribed, vp only affects the transport of the tracer. For boundary conditions providing the tracer concentration in the overlaying water column, vp will also affect to the flux entering through the SWI (see Eq. (9)). This second configuration is studied in Section 3.4. 3.3. Eddy diffusivity The above results suggest that molecular diffusion cannot explain large penetration depths even when porosity and vp take very high values (Table 3). It is worth noting that experimental studies on depth distribution of tracers in sediments use typical slicing resolutions of the order of 1 cm. Eddy diffusivity is generally used for representing the nonlinear turbulent advection of scalar properties of the fluid. The vertical eddy diffusivity coefficient in water bodies can be estimated from time series of temperature vs depth profiles (e.g., Li, 1973; Salas de Le´on et al., 2016). It largely depends on the depth of the water and on the season. For lacustrine environments, empirical values have been reported within the range 10 − 8 to 10 − 3 m 2 s − 1 (e.g., see Salas de Le´on et al., 2016, and references therein). In marine environments the vertical eddy diffusivity can take larger values (of the order of 10 − 2 m 2 s − 1 ), which can be estimated as a function of the water current, the water depth, and the bed friction coefficient (Lane, 2005). The vertical eddy diffusivity in the bottom boundary layer is still several orders of magnitude higher than the molecular diffusion (Wakata, 2018). Thus, assuming a sharp drop from eddy diffusivity above the SWI to only molecular diffusion just below it may not be a realistic description of the natural world. A smooth transition can be more realistic in some cases. Displacements of a large collective of molecules can occur in the pore fluid in the SWI region induced, among other factors, by bottom friction, heat transfer between solids and water, early compaction, and biological activity. The first two factors are favoured in shallow waters. Extreme cases of energy dissipation in the SWI can involve sediment resuspension. In this work, we hypothesize that the energy dissipation in the overlying water column excites the pore water below the SWI, generating an eddy diffusivity that progressively declines with depth. For a first assessment of the effect of eddy diffusivity, the coefficient De in Eq. (8) has been parameterized as follows: 8 k12 =k0 12(1−as,2(t) aS s,2)(16) k0 12 is the initial value of the coefficient (Table 1), and aS s,2 is the saturation concentration at reaction sites 2. Note that Eq. (16) has a local effect, keeping high the uptake kinetics in the front of the cloud of tracers in the dissolved phase. The simulations of Fig. S5 (ESM) have been replicated for tracer C using Eq. (16) with aS s,2 =50,000 a. u. g −1 (a tentative value that is ~16% of the maximum concentration found without saturation for t s = 10 days). The following effects are observed when compared with the absence of saturation: i) The total uptake (i.e. the inventory in the sediment) is initially favoured (by a factor 1.5 in the first day), but then it is reduced by a factor 0.64 and 0.49 after 5 and 10 days, respectively. Ii) The shape of A(z)flattens in the upper layers and the penetration depth slightly increases. Similar results have been found for tracer A. Saturation must be understood in the context of the coexisting major ions in solution that compete for the reaction sites in the solids. El-M- rabet et al. (2001) reported a fitted value of 6.5⋅10 5 Bq/kg for the saturation concentration of Pu in suspended sediments from a reservoir. This value corresponds to about ~10 −4 meq/100 g, which is six orders of magnitude lower than the typical cation exchange capacity of soils and sediments. 3.5. Uptake kinetics at the SWI region in real scenarios (I): Lake Sniardwy, Poland The above results show how low k d tracers can undergo ideal deposition in some sedimentary scenarios, while in others they can show large penetration depths. As an example close to the first situation, Tylmann et al. (2016) reported the 137 Cs profile in a varved sediment core from Lake Zabinskie, Poland, sampled at 40 m water depth. The water column is stratified during most of the year with a well-isolated and anoxic hypolimnion, allowing the formation and preservation of annually laminated sediments. Although porosity was high (~90%), eddy diffusivity is presumed to be absent. The 137 Cs profile was perfectly consistent with the history of atmospheric deposition, without signs of depth penetration, and with its first traces consistently dated by varve counting to AD 1952–1957 (Tylmann et al., 2016). More properly, below this age interval, there was a last measurement in a sediment slice of 1.8 cm thickness, yielding 6.1 ±1.6 Bq kg −1 . The varve date interval was 1949–1952 AD (see their Fig. 4 and data in the supplementary material). This small penetration depth is what can be expected from molecular diffusion (Sections 3.2 and 3.4). This section studies an example for the second situation, with the relatively well-documented case of the Chernobyl fallout radionuclides in Lake Sniardwy (Robbins and Jasinski, 1995). Within four months of the fallout, 134 Cs had penetrated down to about 10 cm in a sediment core. In cores collected from three sites a year later, 134 Cs had penetrated from 14 to 24 cm. 3.5.1. Sedimentary environment and data summary All the information summarized below comes from the work by Robbins and Jasinski (1995). Sniardwy is a shallow lake (5.8 m and 23.1 m of mean and maximum depths, respectively), where frequent vertical mixing occurs favoured by wind-driven currents. The high concentration of suspended solids (annual mean of 230 mg 1 −l ) is maintained primarily by resuspension of sediments. The authors reported data for a set of six cores sampled at three sites: one core at site I, three at site II, and two at site III. They were sampled at elapsed times of 0.31–2.12 years after the Chernobyl accident. Depths at the sampling sites were not reported. They measured activity concentrations versus depth profiles for the Chernobyl fallout radionuclides 134 Cs, 137 Cs and 144 Ce. In all cases, they could disaggregate the contribution of 137 Cs from atmospheric nuclear weapon tests (NWT). For the core from site I (core I 8/87, sampled 1.34 years after the accident), they provided information on porosity and depth profiles of 210 Pb and 226 Ra. In this core, porosity was very high and fitted the analytical function φ(z)=(φ0−φf)e−βz +φf(17) with φ0 =0.978, φf =0.941 and β =0.038 cm −1 . Sediments were relatively homogeneous, with a mixture of clay-silt (70%), fine sand (15%) and coarse materials (15%). They showed a high content of organic matter (39% of dry weight in the surficial sediments and 27% at a depth of 20 cm). The last affects the mean density of solids, which has been set at 2.0 g cm −3 for the present work. The authors reported a total biomass at 9–10 m water depth of 3.5 g m −2 referring to studies carried out at the end of the 1970s. The authors conducted laboratory experiments to study the uptake of Cs in aqueous suspensions of fresh, sieved (<125 μ m), surface sediments, periodically mixed by mechanical agitation. They reported k d values of (4.3 ±0.3)10 3 for a contact time of 130 days for conditions reproducing the surficial sediments. This can represent an overestimation because of the sieving (which excluded larger particle sizes) and the mechanical agitation (which prevents the loss of free surface because of permanent contact among solids). With concentrations of 230 mg l −1 of suspended solids, the authors found that within about a week, short-term equilibration was nearly complete. Robbins and Jasinski (1995) reported few data on 137 Cs concentrations in the water column at six different dates, and used a simple batch reactor model for reconstructing time series of concentrations over the entire range of interest. They neglected the anomalous extreme concentration of 3.1 kBq m −3 found on May 4th , 1986. Here we adopted the same criteria and used a simple analytical fit to the data for defining a suitable a0,w(t)function. From t =0 (April 28th, 1986) until t =10.5 d it fits a Gaussian function, and after that it follows an exponential decline (in units of Bq l −1 and with time in days): a0,w(t)= 7  2 π √3.3e−(9−t)2 2⋅ 3.32;t<10.5d(18a) a0,w(t)=0.763e−0.071(t−10.5);t≥10.5d(18b) The maximum concentration of 0.846 Bq l −1 occurs at t =9 d, and is 0.044 Bq l −1 at t =50.5 d. The simple batch reactor model and the analytical function of Eq. (18) are both oversimplifications. Abril et al. (2004) presented a numerical simulation of the kinetic reactive dispersion of Pu in a lacustrine environment, following a hypothetical atmospheric deposition event. With typical water currents of the order of cm/s and a varying bathymetry, large spatial gradients in concentrations may persist over time scales of several days up to a few weeks. Robbins and Jasinski (1995) used a diagenetic model assuming instantaneous equilibrium between dissolved and particle-bound fractions, described by a uniform and constant k d value. The model included sediment accretion, molecular diffusion, a depth-dependent physical mixing (affecting both solid and dissolved phases), and radioactive decay. Accretion was described in terms of the theory by Berner (1980), which was latter revisited by Abril (2003, 2011). In a simplified version of the model, they neglected accretion and molecular diffusion and assumed a uniform value for the mixing coefficient. The full model was solved for Cs using a simple Gaussian pulse to describe the Chernobyl input, and for 210 Pb, assuming a steady-state profile. This modelling approach could provide an apparently acceptable description for most of the data set, although it left several inconsistencies and open questions: i) Support for the postulated physical mixing. The occurrence of sediment resuspension in shallow areas of the lake does not prove its occurrence at the sampling sites. It has not been quantitatively shown how a sequence of resuspension and deposition events could result in an almost uniform value of the mixing coefficient over a 15 cm depth interval. Bioturbation was also argued as an 9 involved mechanism, although the empirical support is weak. However, it remains unexplained how physical mixing affects Cs and Pb, but not Ce, at sampling sites II and III. ii) 134 Cs profiles in the sequence of cores from site II. A model with equilibrium k d and physical mixing could explain the shape of these profiles, but this required a value of the mixing coefficient which decreased with time. The reasons for this remained unexplained. iii) The treatment of 210 Pb data cannot be seen to be conclusive. Reasons are discrepancies at the top layers, as noticed by the authors, the large involved analytical uncertainties, and the use of Berner’s equations for sediment accretion. A different treatment of the data will be shown later. iv) In all the cases the model failed to simultaneously explain the NWT 137 Cs data. The application of the kinetic reactive transport model to this data set, although limited to a semi-quantitative analysis, will provide a quite different view of the involved processes, as shown below. 3.5.2. 137 Cs profile in core I 8/87 The kinetic reactive transport model has been applied for tracer A (Table 1), used as a proxy for Cs, with boundary conditions given by Eqs. (9) and (18). Porosity is described by Eq. (17), and the bulk density estimated with ρ s =2.0 g cm −3 , because the high percentage of OM in these sediments. Although a varying porosity φ(z)implies changes in the mass of solids and their exposed surfaces (likely linked to tortuosity), both affecting the kinetic coefficient k 11 , a constant value has been adopted for simplicity. The effect of varying concentrations of competing ions has also been neglected because of lack of data. The core was sampled after an elapsed time of 1.34 years. Consequently, and accordingly with the empirical evidence (e.g. Børretzen and Salbu, 2000), the role of irreversibility in the uptake must be explicitly considered. A possibility is to include a third reaction site in solids, as previously commented. An alternative and simpler modelling approach is modifying in Eq. (7) the coefficient for the reverse second reaction: k22 =k0 22e−γt′(19) k0 22 is the initial value (Table 1) and γ is a scaling factor that reduces the desorption rate to its half after a certain elapsed time of the order of 50–100 days. Ideally, the clock for measuring t′should be activated at a given depth, z, when the concentration at sites 2 within the solids reached a certain threshold value. Here, for the sake of simplicity, the clock is simultaneously activated in the entire core taking t′=t. The co-occurrence of high porosity in the sediment and the turbulent mixing by wind-driven currents in this shallow lake allows hypothesizing a certain contribution of eddy diffusivity in the interstitial water, which has been parameterized as follows: De=De,0 1+e α (z−za)(20) Note that this does not require any sediment resuspension. Eq. (20) has similar characteristics to the parameterization used by Robbins and Jasinski (1995): it allows for an almost constant value of eddy diffusivity in the range of depths from 0 to za, followed by a smooth decline to zero (scaled by the coefficient α ) below za. For the sake of simplicity, De is a function only of z, but not of the time, so it is intended to represent the average sedimentary and hydrodynamic conditions. The sediment accumulation rate estimated by the authors from excess 210 Pb returned a value of 0.12 g cm −2 y −1 when excluding flattening in the upper layers (mainly due to early compaction). The figure is 0.17 g cm −2 y −1 after a CF-CS model applied to the entire dataset using the mass depth scale (see details further below). The core length is too short to reach the early compaction limit, so only an upper limit can be estimated for the linear accretion velocity of the SWI using φf from Eq. (17) (see Abril, 2003), which yields v <1.4 cm y −1 . Note that this value is much lower than the >2.5 cm y −1 estimate by the authors above from Berner’s diagenetic equations. As the present interest is exploring the effects of the uptake kinetics, we used Eq. (7) without including sediment accretion. This avoids the need of additional assumptions for implementing the boundary conditions for reaction sites in solids (Eq. (12)). The percolation velocity has been set to zero. Radioactive decay was equally neglected for 137 Cs (T1/2=30.17 years). Several numerical experiments have been conducted with different parameter values. RUN-1 uses De,0 =5 cm −2 h −1 , za =4.5 cm, α =0.7 cm −1 , and γ =0.01 d −1 . The model used a spatial discretization of Δz= 0.25 cm and a time step of 0.003 h, which met the criteria for numerical stability. The depth profiles of activity concentrations, A(z), have been evaluated at different elapsed times after April 28th, 1986. The results are summarized in Figs. 3 and 4. The temporal sequence of A(z)profiles (Fig. 3) shows that total inventories (roughly proportional to the net area under the curve A(z)) first increased (t s =10 d and t s =20 d) and then decreased over time (t s =30 d and further). Although the tracer concentration in the overlaying waters remained higher than that of the pore fluid at the SWI, there was a net diffusive flux entering the sediment (positive values in Fig. 4), but the flux reversed when the concentration in the overlaying waters became lower (negative values in Fig. 4). Note that the concentrations in the overlaying water were prescribed from empirical observations. They were the likely result of a complex reactive transport of Chernobyl fallout radionuclides in the lake, including sorption and desorption from sediments. While the flux remained positive, A(z)showed an exponential pattern, with penetration depths of about 8 cm. After reversal, the concentration in the porefluid decreased in the upper sediment, resulting in desorption from solids, while in the deeper layers the diffusive transport progressed with adsorption. As a result, A(z)progressively flattened and developed a subsurface maximum. The prescribed irreversibility in the second reaction (Eq. (19)) reduced desorption and stabilized the A(z)profile. Thus, after 300 days, the changes became negligible. This numerical experiment highlights the importance of desorption through pore water. The maximum transient inventory in the sediment was 595 Bq cm −2 and was achieved at t =28.7 days. After desorption, the inventory stabilized around 225 Bq cm −2 (data for t =400 days). Sensitivity tests have been conducted to assess the effect of eddy diffusivity and the scaling time in irreversibility. RUN-2 uses De,0 =3 cm −2 h −1 , za =6.5 cm, α =0.5 cm −1 , and γ =0.007 d −1 . RUN-3 uses De,0 =2 cm −2 h −1 , za =8.0 cm, α =0.5 cm −1 , and γ =0.007 d −1 . RUN-2* is as RUN-2, but using the kinetic coefficients of tracer A* (Table 1). The results are summarized in Fig. 5, where the computed A(z)profiles after one year are compared against the empirical data. These last refer to the 137 Cs from Chernobyl. It is worth noting that increasing za does not have a direct effect on the depth of penetration. The magnitude of eddy diffusivity dynamically interplays with the uptake kinetics and the time scale of irreversibility. The goal here was not determining the best set of parameter values that describe the empirical profile, but just getting insights about the relevant basic processes at the SWI region. Fig. 5 suggests that the uptake kinetics coupled to the molecular and eddy diffusivity in the pore fluid can account for the main features of the empirical profile without invoking the physical mixing of solids. Major discrepancies are the scaling factor for concentrations with tracer A and the tail beyond 12 cm depth. Such discrepancies can be at least partially explained by the series of simplifying assumptions included in the model setup. Including accretion can increase the depth of penetration by ~1.5 cm, as commented on above. Analytical description of the tracer concentration in the water column (Eq. (18)) is a likely oversimplification. The set of kinetic coefficients for tracer A does not exactly correspond to that for Cs in this lake, which seems to be better described by tracer A*. The dependence of these coefficients with depth has not been considered. A decreasing 16 processes in sediments. A reviewed basis for radiometric dating models. J. Paleolimnol. 30, 363–370. Abril, J.M., 2011. Could bulk density profiles provide information about recent sedimentation rates? J. Paleolimnol. 46, 173–186. https://doi.org/10.1007/s10933- 011-9520-2. Abril, J.M., 2016. A 210 Pb-based chronological model for recent sediments with random entries of mass and activities: model development. J. Environ. Radioact. 151, 64–74. Abril, J.M., 2020. Multimodal-TERESA, a 210 Pb-based radiometric dating model for recent sediments under largely varying rates of supply. Quat. Geochronol. 55, 101032. Abril, J.M., Brunskill, G.J., 2014. Evidence that excess 210 Pb flux varies with sediment accumulation rate and implications for dating recent sediments. J. Paleolimnol. 52, 121–137. Abril, J.M., El-Mrabet, R., Barros, H., 2004. The importance of recording physical and chemical variables simultaneously with remote radiological surveillance of aquatic systems: a perspective for environmental modelling. J. Environ. Radioact. 72, 145–152. Abril, J.M., Fraga, E., 1996. Some physical and chemical features of the variability of kd distribution coefficients for radionuclides. J. Environ. Radioact. 30 (3), 253–270. Abril, J.M., Gharbi, F., 2012. Radiometric dating of recent sediments: beyond the boundary conditions. J. Paleolimnol. 48, 449–460. Abril, J.M., García-Le´ on, M., 1991. A mathematical approach for modelling radionuclide dispersion in the marine environment. J. Environ. Radioact. 13, 39–54. Appleby, P.G., 1998. Dating recent sediments by 210 Pb: problems and solutions. In: Illus, E. (Ed.), Dating of Sediments and Determination of Sedimentation Rate. STUK A-145, pp. 7–24. Finland. Barnett, R.L., Newton, T.L., Charman, D.J., Gehrels, W.R., 2017. Salt-marsh testate amoebae as precise and widespread indicators of sea-level change. Earth Sci. Rev. 164, 193–207. Barros, H., Abril, J.M., 2004. Experimental and modelling study on the uptake and desorption kinetics of 133 Ba by suspended estuarine sediments from southern Spain. Water Res. 38 (3), 749–755. Barros, H., Abril, J.M., 2005. Constraints in the construction and/or selection of kinetic box models for the uptake of radionuclides and heavy metals by suspended particulate matter. Ecol. Model. 185, 371–385. Barros, H., Abril, J.M., 2008. Kinetic box models for the uptake of radionuclides and heavy metals by suspended particulate matter: equivalence between models and its implications. J. Environ. Radioact. 99, 146–158. Barros, H., Abril, J.M., El-Mrabet, R., Laissaoui, A., 2003. Kinetically controlled radionuclide sorption by sediment cores from two different environments. Experimental studies using 133 Ba as a tracer. In: “Recent Advances in Multidisciplinary Applied Physics”. Proceedings of the APHYS-2003 International Meeting, Badajoz, Spain, pp. 531–542. Barros, H., Laissaoui, A., Abril, J.M., 2004. Trends of radionuclide sorption by estuarine sediments. Experimental studies using 133 Ba as a tracer. Sci. Total Environ. 319 (1–3), 253–267. Bellucci, L.G., Frignani, M., Cochran, J.K., Albertazzi, S., Zaggia, L., Cecconi, G., Hopkins, H., 2007. 210 Pb and 137 Cs as chronometers for salt marsh accretion in the Venice Lagoon and links to flooding frequency and climate change. J. Environ. Radioact. 97, 85–102. Benoit, G., Hemond, H.F., 1991. Evidence for diffusive redistribution of 210 Pb in lake sediments. Ceochim. Cosmochim. Ac. 55, 1963–1975. Berkowitz, B., Dror, I., Hansen, S.K., Scher, H., 2016. Measurements and models of reactive transport in geological media. Rev. Geophys. 54, 930–986. Berner, R.A., 1980. Early Diagenesis: A Theoretical Approach. Princeton University Press, Princeton, NJ. Børretzen, P., Salbu, B., 2000. Estimation of apparent rate coefficients for radionuclides interacting with marine sediments from Novaya Zemlya. Sci. Total Environ. 262, 91–102. Børretzen, P., Salbu, B., 2002. Fixation of Cs to marine sediments estimated by a stochastic modelling approach. J. Environ. Radioact. 61, 1–20. Botwe, B.O., Abril, J.M., Schirone, A., Barsanti, M., Delbono, I., Delfanti, R., Nyarko, E., Lens, P.N.L., 2017. Settling fluxes and sediment accumulation rates by the combined use of sediment traps and sediment cores in Tema Harbour (Ghana). Sci. Total Environ. 609, 1114–1125. Boudreau, B.P., 2000. The mathematics of early digenesis: from worms to waves. Rev. Geophys. 38 (3), 389–416. Comans, R.N.J., Middelburg, J.J., Zonderhuis, J., Woittiez, J.R.W., De Lange, G.J., Das, H.A., Van Der Weijden, C.H., 1989. Mobilization of radiocaesium in pore water of lake sediments. Nature 339, 367–369. Comans, R., Hilton, J., Cremers, A., Geelhoed-Bonouvrie, P.A., Smith, J., 1997. Interpreting and predicting in situ distribution coefficients of radiocaesium in aquatic systems. Stud. Environ. Sci. 68, 129–140. Cundy, A.B., Croudace, I.W., 1996. Sediment accretion and recent sea-level rise in the Solent, Southern England: inferences from radiometric and geochemical studies. Estuar. Coast Shelf Sci. 43, 449–467. Cundy, A.B., Croudace, I.W., 2017. The fate of contaminants and stable Pb isotopes in a changing estuarine environment: 20 years on. Environ. Sci. Technol. 51, 9488–9497. Duursma, E.K., Carroll, J., 1996. Environmental Compartments: Equilibria and Assessment of Processes between Air, Water, Sediments, and Biota. Springer, Berlin. modelling study of El Mrabet, R., Abril, J.M., G Manj´on, G., Garca-Tenorio, R., 2001. Experimental and modelling study of plutonium uptake by suspended matter in aquatic environments from southern Spain. Water Res. 35 (17), 4184–4190. El-Mrabet, R., Abril, J.M., 241 Manj´on, G., Garca-Tenorio, R., 2004. Experimental and Am uptake by suspended matter in freshwater environments from southern Spain. J. Radioanal. Nucl. Chem. 261 (1), 137–144. Eisma, D., 1993. Suspended Matter in the Aquatic Environemnt. Springer, Berlin, Heidelberg, ISBN 0-387-55825-X. Glassley, W.E., Nitao, J.J., Grant, C.W., 2003. Three-dimensional spatial variability of chemical properties around a monitored waste emplacement tunnel. J. Contam. Hydrol. 62–63, 495–507. Holby, O., Evans, S., 1996. The vertical distribution of Chernobyl-derived radionuclides in a Baltic Sea sediment. J. Environ. Radioact. 33, 129–145. IAEA, 2004. Sediment distribution coefficients and concentration factors for biota in the marine environment. IAEA Technical Reports Series No. 422, Vienna. ISSN 0074–1914. IAEA, 2010. Handbook of parameter values for the prediction of radionuclide transfer in terrestrial and freshwater environments. IAEA Technical Reports Series No. 472, Vienna. ISSN 0074–1914. Iurian, A.R., Millward, G.E., Blake, W., Abril Hern´ andez, J.M., 2021. Fine-tuning of 210 Pb-based methods for dating vegetated saltmarsh sediments. Quat. Geochronol. 62, 101153. https://doi.org/10.1016/j.quageo.2021.101153. Kennish, M.J., 2001. Coastal salt marsh systems in the U.S.: a review of anthropogenic impacts. J. Coast Res. 17, 731–748. Kramer, K.J.M., Misdorp, R., Berger, G., Duijts, R., 1991. Maximum pollutant concentrations at the wrong depth: a misleading pollution history in a sediment core. Mar. Chem. 36, 183–198. Krezoski, J.R., Robbins, J.A., White, D.S., 1984. Dual radiotracer measurement of zoobenthos mediated solute and particle transport in freshwater sediments. J. Geophys. Res. 89, 7939–7947. Lane, A., 2005. Development of a Lagrangian sediment model to reproduce the bathymetric evolution of the Mersey Estuary. Ocean Dynam. 55, 541–548. Li, Y.H., 1973. Vertical eddy diffusion coefficient in Lake Zürich. Schweiz. Z. Hydrol. 35, 1–7. https://doi.org/10.1007/BF02502061. Mabit, L., Benmansour, M., Abril, J.M., Walling, D.E., Meusburger, K., Iurian, A.R., Bernard, C., Tarjan, S., Owens, P.N., Blake, W.H., Alewell, C., 2014. Fallout 210 Pb as a soil and sediment tracer in catchment sediment budget investigations: a review. Earth Sci. Rev. 138, 335–351. Maderich, V., Jung, K.T., Brovchenko, I., Kim, K.O., 2017. Migration of radioactivity in multi-fraction sediments. Environ. Fluid Mech. 17, 1207–1231. Marchand, C., 2017. Soil carbon stocks and burial rates along a mangrove forest chronosequence (French Guiana). For. Ecol. Manag. 384, 92–99. MacQuarrie, K.T.B., Mayer, K.U., 2005. Reactive transport modeling in fractured rock: a state-of-the-science review. Earth Sci. Rev. 72 (3–4), 189–227. Mudbidre, R., Baskaran, M., Schweitzer, L., 2014. Investigations of the partitioning and residence times of Po-210 and Pb-210 in a riverine system in Southeast Michigan USA. J. Environ. Radioact. 138, 375–383. Novikov, A.P., Kalmykov, S.N., Utsunomiya, S., Ewing, R.C., Horreard, F., Merkulov, A., Clark, S.B., Tkachev, V.V., Myasoedov, B.F., 2006. Colloid transport of plutonium in the far-field of the Mayak production association, Russia. Science 314, 638–641. Nyffeler, U.P., Santschi, P.H., Li, Y.H., 1986. The relevance of scavenging kinetics to modeling of sediment-water interactions in natural waters. Limnol. Oceanogr. 31 (2), 277–292. Otosaka, S., Kambayashi, S., Fukuda, M., Tsuruta, T., Misonou, T., Suzuki, T., Aono, T., 2020. Behavior of radiocesium in sediments in Fukushima coastal waters: verification of desorption potential through the pore water. Environ. Sci. Technol. 54, 13778–13785. Peri´ a˜ nez, R., 2005. Modelling the Dispersion of Radionuclides in the Marine Environment: an Introduction. Springer, Berlin, Germany. ISBN 3-540-24875-7. Petersen, W., Knauth, H.D., Pepelnik, R., 1990. Vertical distribution of Chernobyl isotopes and their correlation with heavy metals and organic carbon in sediment cores of the Elbe Estuary. Sci. Total Environ. 97, 531–547. Robbins, J.A., 1986. A model for particle-selective transport of tracers in sediments with conveyor belt deposit feeders. J. Geophys. Res. 91 (C7), 8542–8558. Robbins, J.A., Jasinski, A.W., 1995. Chernobyl fallout radionuclides in lake Sniardwy, Poland. J. Environ. Radioact. 26, 157–184. Salas de Le´ on, D.a., Alcocer, J., Ardiles Gloria, V., Quiroz-Martínez, B., 2016. Estimation of the eddy diffusivity coefficient in a warm monomictic tropical Lake. J. Limnol. 75 (s1), 161–168. https://doi.org/10.4081/jlimnol.2016.1431. Scholkovitz, E.R., Mann, D.R., 1984. The pore water chemistry of 239,240 Pu and 137 Cs in sediments of Buzzards Bay, Massachusetts. Geochem. Cosmochim. Acta 48, 1107–1114. Smith, J.T., Comans, R.N.J., 1996. Modelling the diffusive transport and remobilisation of 137 Cs in sediments: the effects of sorption kinetics and reversibility. Geochem. Cosmochim. Acta 60 (6), 995–1004. Smith, J.T., Comans, R.N.J., Ireland, D.G., Nolan, L., Hilton, J., 2000. Experimental and in situ study of radiocaesium transfer across the sediment-water interface and mobility in lake sediments. Appl. Geochem. 15, 833–848. Steefel, C.I., DePaolo, D.J., Lichtner, P.C., 2005. Reactive transport modeling: an essential tool and a new research approach for the Earth sciences. Earth Planet Sci. Lett. 240 (3–4), 539–558. Taylor, A., Blake, W.H., Couldrick, L., Keith-Roach, M.J., 2017. Sorption behaviour of beryllium-7 and implications for its use as a sediment tracer. Geoderma 187–188, 16–23. Abril, J.M., 1998b. Basic microscopic theory of the distribution, transfer and uptake kinetics of dissolved radionuclides by suspended particulate matter. Part II: applications. J. Environ. Radioact. 41, 325–341. Abril, J.M., 2003. A new theoretical treatment of compaction and the advective-diffusive 17 dating results with varve chronology and independent chronostratigraphic markers: problems and implications. Quat. Geochronol. 32, 1–10. Wakata, Y., 2018. LES study of vertical eddy diffusivity estimation in bottom boundary layers. J. Phys. Oceanogr. 48 (8), 1903–1920. https://doi.org/10.1175/JPO-D-17- 0165.1. partitioning of Pb, Bi and Wang, J., Zhong, Q., 210 Baskaran, 207 M., 210 Du, J., 2019. Investigations on the time-series Po between marine particles and solution under different salinity and pH conditions. Chem. Geol. 528, 119275 https://doi.org/ 10.1016/j.chemgeo.2019.119275. Yoshimura, K., Onda, Y., Fukushima, T., 2014. Sediment particle size and initial radiocesium accumulation in ponds following the Fukushima DNPP accident. Sci. Rep. 4, 4514. https://doi.org/10.1038/srep04514. Tylmann, W., Alicja Bonk, A., Goslar, T., Wulf, S., Grosjean, M., 2016. Calibrating 210 Pb