Stepping beyond perfectly mixed conditions in soil hydrological modelling using a Lagrangian approach

A recent experiment of Bowers et al. (2020) revealed that diffusive mixing of water isotopes (δ²H and δ¹⁸O) over a fully saturated soil sample of a few centimetres in length required several days to equilibrate completely. In this study, we present an approach to simulate such time-delayed diffusive mixing processes, on the pore scale, beyond instantaneously and perfectly mixed conditions. The diffusive pore mixing (DIPMI) approach is based on a Lagrangian perspective on water particles moving by diffusion over the pore space of a soil volume and carrying concentrations of solutes or isotopes. The idea of DIPMI is to account for the self-diffusion of water particles across a characteristic length scale of the pore space using pore-size-dependent diffusion coefficients. The model parameters can be derived from the soil-specific water retention curve, and no further calibration is needed. We test our DIPMI approach by simulating diffusive mixing of water isotopes over the pore space of a saturated soil volume using the experimental data of Bowers et al. (2020). Simulation results show the feasibility of the DIPMI approach for reproducing the measured mixing times and concentrations of isotopes at different tensions over the pore space. This result corroborates the finding that diffusive mixing in soils depends on the pore size distribution and the specific soil water retention properties. Additionally, we perform a virtual experiment with the DIPMI approach by simulating mixing and leaching processes of a solute in a vertical, saturated soil column and compare the results against simulations with the common perfect mixing assumption. The results of this virtual experiment reveal that the frequently observed steep rise and long tailing of breakthrough curves, which are typically associated with non-uniform transport in heterogeneous soils, may also occur in homogeneous media as a result of imperfect subscale mixing in a macroscopically homogeneous soil matrix.

1 Introduction

Water isotopes are widely used as tracers to investigate a variety of hydrological processes (Sprenger et al., 2016). While they were originally used to separate pre-event and event water contributions to storm runoff (Bonell et al., 1990; Sklash et al., 1996), they are now more frequently considered as a continuous source of information to infer the travel time distributions of water through hydrological systems (e.g. McGlynn et al., 2002; McGlynn and Seibert, 2003; Weiler et al., 2003; Klaus and McDonnell, 2013). Early analyses often relied on time-invariant transfer functions, whereas some of the more recent approaches are time-dependent and, for example, use an age-ranked storage as a “state” variable in combination with StorAge Selection (SAS) functions for streamflow and evapotranspiration to infer their respective travel time distributions (Harman, 2015; Rodriguez and Klaus, 2019; Rodriguez et al., 2021). This inference of transit times from water isotopes commonly implies a distinct relation between water age and its isotopic composition.

However, recent laboratory and field experiments suggest that this relation and the fate of water isotopes in the soil–plant–atmosphere system may in fact be more complex than frequently assumed. Mennekes et al. (2021), for example, used in situ probes to measure isotopic signatures of water in soil and tree xylem, during tracer irrigation experiments on the plot scale, and discussed that the travel times of water fractions in soils and plants may be distinctly different. This is in line with the findings of Benettin et al. (2021), who performed lysimeter experiments with isotopic-labelled waters to not only close, but also trace, all fluxes in the water balance. They found that the isotopic composition of transpiration fluxes was significantly different compared to breakthrough fluxes in soil drainage. On the pore scale, Orlowski and Breuer (2020) investigated how the isotopic composition of water depends on the retention characteristic of a soil. Their experimental results highlight “a need to better characterize processes that govern isotope fractionation with respect to soil water retention characteristics” because they found fractionation of δ²H and δ¹⁸O isotopes during their diffusive movement over different pore sizes, especially under high tensions in small pores. In this context, Bowers et al. (2020) performed an experiment using a combination of extraction methods to sample isotopically defined water fractions from a saturated soil sample over the complete water retention curve to explore how fast water isotopes (δ²H and δ¹⁸O) mix diffusively over the entire pore size distribution. They showed that mixing and fractionation processes of water isotopes depend on different tensions at which water is held in pores of different size. The most interesting insight of the experiment was that the isotope tracer required up to 3–4 d until it was distributed uniformly over the entire pore space, even though the studied soil sample was only a few centimetres in length.

In particular, the experimental findings of Bowers et al. (2020) suggest that ignoring self-diffusion processes of water isotopes within the soil pore space can result in incorrect estimates of water ages and travel times, which emphasizes the requirement for including these pore-scale processes into soil hydrological models. Common soil hydrological models average over pore-size-dependent differences in the flow field and concentration gradients in control volumes (Berkowitz et al., 2016) to describe diffusive mixing of water and solutes. This implies that incoming “new” event water and “old” pockets of pre-event water in soil mix perfectly and instantaneously over the subscale pore size distribution in a single time step. This common perfect mixing assumption is, thus, not in line with the recent experimental findings (e.g. Bowers et al., 2020; Orlowski and Breuer, 2020). Further studies also have shown that these different pockets of water may indeed co-exist, even in close spatial distances, without perfect mixing and that water and solutes repeatedly travel through the same pathways (memory effect), even after several infiltration cycles of new precipitation (Gouet-Kaplan and Berkowitz, 2011; Kapetas et al., 2014). In this way, the establishment of stable water pockets in soils is possible, which may comprise significantly different isotopic and chemical compositions, depending on the properties of infiltrating water. This imperfect mixing of water and solutes in the pore space is frequently discussed but instead in the context of rapid preferential flow in macroporous structures (Beven and Germann, 1982, 2013), which is also commonly assumed to be the main reason for the characteristic steep rise and long tailing of corresponding breakthrough curves (e.g. Berkowitz et al., 2006; Edery et al., 2014). Based on the findings of Bowers et al. (2020), we hypothesize that imperfect mixing of water and solutes over the pore sizes of a macroscopically homogeneous and saturated soil matrix will also yield such typical shapes of breakthrough curves, even without the presence of macroporous soil structures.

To account for subscale diffusive mixing of solutes or water isotopes over pore sizes, in line with the findings of Bowers et al. (2020), we propose that the recent particle-based Lagrangian approaches (e.g. Berkowitz et al., 2006; Zehe and Jackisch, 2016; Jackisch and Zehe, 2018; Engdahl et al., 2017, 2019; Schmidt et al., 2019) offer a series of new possibilities in this regard. Zehe and Jackisch (2016) showed that the conceptualization of fluid flow in partially saturated soils as a Lagrangian advective–diffusive random walk of water particles is feasible for successfully reproducing observed soil water dynamics and distinguishing explicitly pre-event and event waters. The key was to account for a variable, pore-size-dependent mobility of water particles, which was achieved by discretizing the pore space into pores of different sizes with specific hydraulic conductivities and water diffusivities (cf. Sect. 2.1). In follow-up studies (Sternagel et al., 2019, 2021), we extended this model approach and developed the Lagrangian Soil Water and Solute Transport (LAST)-Model for simulations of (reactive) solute transport combined with water motion in heterogeneous, partially saturated 1-D soil domains. These former versions of the LAST-Model, however, assumed the instantaneous, perfect mixing of solutes among water particles in a control volume, which implied that the model may have smoothed out concentration gradients too quickly (Sternagel et al., 2021).

In this study, we eliminate this perfect mixing assumption and introduce the diffusive pore mixing (DIPMI) approach to provide a Lagrangian method to improve our ability to describe diffusive mixing processes on the pore scale. The idea of DIPMI is to account for the self-diffusion of water particles across a characteristic length scale of the pore space, using pore-size-dependent diffusion coefficients. Its model parameters can be derived from the soil-specific water retention curve, and no further calibration is needed. We initially test the DIPMI approach by simulating the experiment of Bowers et al. (2020), using the respective dataset. Furthermore, we implement the DIPMI approach into our LAST-Model framework and perform a virtual experiment to test our hypothesis. To this end, we simulate diffusive mixing and the breakthrough of a representative solute in a vertical 1-D soil column during a steady-state saturated flow and compare the results to simulations using the perfect mixing assumption.

2 Lagrangian approach for soil hydrological and subscale diffusion processes

2.1 Underlying concept of the LAST-Model

The Lagrangian perspective describes a mobile observer travelling along the trajectory of a fluid particle through a system (Currie, 2016). As mentioned above, we have applied the Lagrangian perspective before, in our LAST-Model (Sternagel et al., 2019, 2021), to describe the vertical displacement of water particles with related (reactive) solute transport in interacting domains of soil matrix and macropores. Water particles are defined discretely by constant water mass and volume. They additionally carry time-dependent information about, for example, their vertical position in both domains and solute concentrations. The two flow domains of soil matrix and macropores are vertically subdivided into layers. This vertical discretization is required to quantify and translate the number of water particles, in combination with the water particle mass and density, into a soil water content per vertical soil layer. The soil water content in turn corresponds to the sum of the volume fractions of soil water, which are stored in soil pores of different sizes. Water particles travel at different velocities in these pore fractions that are characterized by the shape of the water diffusivity and hydraulic conductivity curve. These curves are partitioned into a certain number of pore size classes or bins (“binning”) between the residual and saturated water content. Depending on the pore size class/bin in which a water particle is located, it experiences different displacements in the vertical direction by means of pore-size-specific advection and diffusivity, i.e. the water particles in smaller pores experience a smaller vertical displacement step than those in coarser pores. Hence, this approach accounts for the combined effects of gravity and capillarity on water flow in partially saturated soils, as well as the subscale variability in flow velocities across different pore sizes (Zehe and Jackisch, 2016). However, in former versions of LAST (Sternagel et al., 2019, 2021), we assumed that the timescale for diffusive mixing is smaller than the simulation time step, and hence, solutes perfectly and instantaneously mix over all pore size classes/bins in a soil layer. Thus, after the non-uniform, vertical movement of particles, the solute concentrations were averaged over all the present water particles within a vertical soil layer per time step (perfect mixing assumption).

Furthermore, the LAST-Model allows for the simulation of sorption and degradation processes during the transport of reactive substances. Non-linear adsorption and desorption processes are realized by an explicit transfer of dissolved solute mass between water particles and surrounding sorption sites of the soil phase in a certain depth. Adsorbed solute masses are then degraded by means of first-order kinetics.

Previous test simulations revealed that the LAST-Model effectively describes solute concentration profiles and leaching behaviours of both conservative tracers and reactive substances, both on the plot and field scale, under various flow conditions. In particular, the structural macropore domain of LAST is an asset in the capturing of the typical pattern of the preferential bypassing of solutes in macroporous soils. Despite these promising results, we also showed that our former assumption of perfect mixing of solutes within a vertical soil layer was a strong simplification that could lead to smoothing of pore-size-dependent differences in the flow field and concentration gradients, called over-mixing (Sternagel et al., 2021).

2.2 The DIPMI approach: concept to represent subscale diffusion in a Lagrangian model

In this study, we step beyond the use of the perfect mixing assumption by developing a Lagrangian approach to simulate the self-diffusive mixing of water and solutes over the pore space. We explain this diffusive pore mixing (DIPMI) approach based on the schematic sketch in Fig. 1.

Figure 1 Schematic sketch of the diffusive pore mixing (DIPMI) approach. See the descriptions of Eqs. (1)–(4) in Sect. 2.2 for further information on the parameters.

The rectangle in the left box at t₀ schematically illustrates a control volume with height dz (e.g. a soil layer) of a 1-D vertical soil profile with total depth z. The width of this rectangle illustrates the entire subscale extent of pore space L_D in which fluid particles can move by self-diffusion. L_D represents a characteristic flow length in the pore space, which is related to the tortuosity of flow paths, the subscale distribution of pore sizes, and, thus, to the soil-specific water retention curve (see the right box at t₀). The extent of the pore space L_D and the soil water retention curve is subdivided equally into a certain number N of bins i, which represents water storage in different pore size classes with corresponding matric potentials ψ. N generally depends on soil-specific properties (Talbot and Ogden, 2008; Ogden et al., 2017), and we assign N=200. This value of N is in line with Talbot and Ogden (2008), who used a comparable method and suggested that the soil moisture domain of most soil types can be discretized sufficiently by 200 bins. Furthermore, Zehe and Jackisch (2016), who used a similar Lagrangian approach to simulate soil water dynamics, performed an analysis of the sensitivity of N and found that N>50 is favourable for producing good simulation results, compared to a Richards solver.

Based on the Young–Laplace equation and the subdivided (“binned”) soil water retention curve, we can determine the total subscale extent of the pore space L_D (L) by the integral (Eq. 2) over the corresponding distribution of pore radii r_i (Eq. 1), as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{r}_i={\displaystyle \frac{-\mathrm{2}\cdot \mathit{\sigma }}{g\cdot \mathit{\rho }\cdot \mathit{\psi }\left(i\right)}},\text{(2)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{D}}=\underset{i=\mathrm{1}}{\overset{N}{\int }}{r}_i,\end{array}$$

where r_i (L) is the radius of a pore size class, σ (F L⁻¹) is the surface tension of fluid, g (L T⁻²) is the gravitational acceleration of the Earth, ρ (M L⁻³) is the fluid density, and ψ(i) (L) is the matric potential of a pore size class/bin derived from the soil water retention curve. Hence, the Young–Laplace equation represents a connection between measurable matric potentials and corresponding pore radii of each pore size class.

Each bin is defined by a constant width $\mathit{\delta }{L}_{\mathrm{D}}=L_{\mathrm{D}}\cdot N^{-\mathrm{1}}$ (L) and a corresponding location within L_D. In our example, each bin is saturated by fluid particles carrying two different isotopic signatures (illustrated by the green and dark yellow particles). These fluid particles have a position within L_D, and accordingly, they are located within a certain bin (i=1 … N). This means that, at t₀, the pore space is filled by fluid particles with different isotopic signatures. Coarse and medium pore size classes/bins (on the left) are filled by particles with one isotopic signature (green particles), while small pore size classes/bins (on the right) are filled by particles with another isotopic signature (dark yellow particles). Hence, fluid particles with different isotopic signatures are distinctly unmixed in the pore space before the self-diffusive mixing starts at t₁.

At t₁, the self-diffusive mixing of particles with different isotopic signatures begins. For each particle, a displacement step Δd_LD (L) along L_D is calculated using a random walk equation (Eq. 3), which is then subtracted from the current position of the particle.

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{d}_{L_{\mathrm{D}}}=Z\cdot \sqrt{\mathrm{2}\cdot D\left(i\right)\cdot \mathrm{d}t}-\left({\displaystyle \frac{\mathit{\delta }D\left(i\right)}{\mathit{\delta }{L}_{\mathrm{D}}}}\right),\end{array}$$

where Z [−1, 1] is a random number drawn from a standard normal distribution, D(i) (L T⁻¹) is the diffusion coefficient in a certain bin or pore size class, and dt (T) is the time step. The last term $\frac{\mathit{\delta }D\left(i\right)}{\mathit{\delta }{L}_{\mathrm{D}}}$ is a correction term to avoid artefacts in case of spatially variable diffusion coefficients (see the explanation below).

The random number between −1 and 1 enables the displacement of particles in the positive and in the negative direction along L_D, representing the undirected process of molecular self-diffusion due to Brownian motion. In this diffusion process, the particles are not displaced by the same diffusion coefficient D. Depending on the bin or pore size class in which a particle is located, it experiences a specific diffusivity. Each bin/pore size class has its own D(i) value, which is determined by its proportion on the total soil porosity (proportion factor), multiplied with the molecular diffusion coefficient of free water = $\mathrm{2.272}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹ (after Mills, 1973; Eq. 4) as follows:

$$\begin{array}{}\text{(4)}& D\left(i\right)=\mathrm{2.272}\times {\mathrm{10}}^{-\mathrm{9}}\cdot \left({\displaystyle \frac{\mathit{\theta }\left(i\right)-{\mathit{\theta }}_{\mathrm{r}}}{\mathit{\varphi }}}\right),\end{array}$$

where the proportion factor comprises the respective water content of a specific bin θ(i) (–) according to the binning of the water retention curve (cf. Fig. 1), the residual water content θ_r (–), and the total soil porosity ϕ (–).

In this way, larger pores/bins have higher D values, and thus, particles experience a larger diffusive displacement in these pores/bins, while particles in smaller pores/bins with lower D values experience a smaller diffusive displacement. This reflects the decline of the free path length for Brownian motions in smaller pores. With this approach of variable, i.e. pore-size-dependent diffusion coefficients, we account for the general controls in the diffusion rate in soil solution, e.g. the diffusion coefficient of a certain fluid or substance, pore size or water content, and tortuosity of flow paths (Chou et al., 2012). With these variable diffusivities in pore size classes, we add a correction term $\frac{\mathit{\delta }D\left(i\right)}{\mathit{\delta }{L}_{\mathrm{D}}}$ (Zehe and Jackisch, 2016) to the random walk equation (Eq. 3) to avoid artificial particle accumulation in the smallest pores/bins, as stated by Uffink (1990). The random number Z in Eq. (3) is either positive or negative and determines, in this way, in which direction a fluid particle is displaced along L_D, i.e. if it moves diffusively in the direction of smaller or coarser pore size classes/bins (minus sign shows direction of the coarse pores; plus sign shows the direction of the small pores). At the same time, the correction term has a constant negative sign, and thus, the diffusive displacement steps Δd_LD of particles in the direction of coarse pores are enhanced, while they are diminished in the direction of smaller pores. In this way, calculated displacement step lengths of particles are corrected with different strength. The overall greater displacement lengths of particles in coarser pores (due to higher D values) in the direction of smaller pores are balanced, and consequently, artificial particle accumulations in the smallest pores are prevented. According to its displacement step, a particle is assigned a new position within L_D, and if Δd_LD>δL_D, then the particle is also assigned a new bin number. At the left and right boundary of the entire pore space L_D, particles are reflected into the pore space to avoid particle accumulation at the boundaries.

Finally, after a certain mixing time t₂, the pore space L_D in our example has reached a final equilibrium state with an uniform isotopic signature in all bins. The subscale separation (binning) of the pore space allows for the calculation of mixing concentrations in single bins or tension areas (i.e. a certain number of adjacent bins/pore size classes within defined ranges of matric potentials).

3 Testing the DIPMI approach

We test our DIPMI approach by simulating the experiment of Bowers et al. (2020) with diffusive mixing of water isotopes over the pore space of a fully saturated soil volume. Furthermore, we perform a virtual experiment by simulating the mixing and leaching processes of a representative solute in a vertical, saturated soil column and compare results of the DIPMI approach against simulations that employ the common perfect mixing assumption.

Table 1 Experimental and simulation parameters and soil hydraulic parameters, after van Genuchten (1980) and Mualem (1976), of sandy loam, where θ_s is the saturated soil water content, θ_r is the residual soil water content, α is the inverse of an air entry value, and n is a quantity characterizing the pore size distribution.

3.1 Simulating the experiment of Bowers et al. (2020)

3.1.1 Original experiment

Bowers et al. (2020) used a combination of extraction methods to quantify the time-dependent mixing of different water isotopes (δ²H and δ¹⁸O) held at different tensions in fully saturated soil samples over the entire water retention curve. Their objective was to analyse how separate soil water fractions, stored in different pore sizes, interact by self-diffusion. They took oven-dried, homogenized soil samples (18–30 g) of a sandy loam (cf. Table 1) and initially wetted them with isotopically light water (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{130}$ ‰ and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{17.6}$ ‰) to a relative saturation of about 16 %–17 %. This water fraction represented an initial water content stored at high matric potentials in the smallest pores. The remaining free pore space was then completely saturated with isotopically heavy water (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{44}$ ‰ and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{7.8}$ ‰), representing new incoming water. Soil samples were then placed into horizontal cylinders and different equilibration time periods of 0 h, 8 h, 1 d, 3 d, and 7 d were applied to enable mixing of the two isotopically distinct waters over the pore space by pure self-diffusion (no advection). After each time period, the soil water samples were sequentially extracted from the soil samples at the following three subsequent tensions: (i) centrifugation at $\sim <\mathrm{0.016}$ MPa for waters in low-tension areas, (ii) centrifugation at ∼0.016–1.14 MPa for waters in mid-tension areas, and finally (iii) cryogenic vacuum distillation (CVD) at 100 MPa to capture residual waters in high-tension areas ($\sim >\mathrm{1.14}$ MPa). Isotopic compositions of extracted water samples were then analysed to assess the differences between the diffusive mixing behaviour in the three tension areas over 7 d. All experimental data with detailed information about soil hydraulic properties of the sandy loam are freely accessible via an Open Science Framework (Bowers and Mercer, 2020).

3.1.2 Simulation with the DIPMI approach

For our simulations with the DIPMI approach, we assume a representative, saturated soil layer volume with the same soil hydraulic properties and soil water retention characteristics as the sandy loam used by Bowers et al. (2020). From the given soil water retention curve in the study of Bowers et al. (2020), it is possible to infer the pore diameters corresponding to the respective matric potential ranges of the sandy loam (cf. Eq. 1; Sect. 2.2). The representative soil layer is defined only by the extent of the pore space L_D (cf. Sect. 2) and has no vertical extent, as we simulate the pure diffusion over the pore space without any vertical displacement of particles in this case. For saturation of the pore space, we generally use the same saturation procedure but we do not use the pure isotopically light and heavy waters as in the experiments (cf. Sect. 3.1.1). Instead, we use the upper and lower standard deviation (SD) values of isotopic concentrations, which Bowers et al. (2020) measured after the first extraction time at 0 h. This is necessary to enable an equal initial condition of isotopic concentrations in simulations compared to observation; we, thus, perform simulations with differing initial isotopic values for light and heavy water (cf. Table 1). First, we use the upper SD values for light water (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{79}$ ‰ and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{9.3}$ ‰) and heavy water (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{46}$ ‰ and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{7.2}$ ‰). Second, we use the lower SD values for light water (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{99}$ ‰ and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{12.3}$ ‰) and heavy water (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{48}$ ‰ and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{7.8}$ ‰). After both simulation runs, isotopic concentrations are averaged over the respective tension areas at each time point, resulting in the mean and SD values of our simulations (cf. Table 2). We use the given experimental data of the soil water retention curve and Eqs. (1) and (2) to quantify the total extent of pore space L_D of the soil volume by 21 000 µm and subdivide it into 200 bins or pore size classes (cf. Sect. 2.2). Additionally, we repeat the simulations with the (i) constant diffusivity D of $\mathrm{2.272}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹ (diffusion coefficient of free water) in all bins and the (ii) linear pore-size-distributed diffusivities D over all bins calculated by Eq. (4) (cf. Sect. 2.2). Pore-size-distributed D values thus range from $\sim \mathrm{1.9}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹ in bin 1 (i.e. the largest pore in a low-tension area) to $\sim \mathrm{9.0}\times {\mathrm{10}}^{-\mathrm{12}}$ m² s⁻¹ in bin 200 (i.e. the smallest pore in a high-tension area). We do not distinguish between specific diffusion coefficients for δ²H and δ¹⁸O, as Hasegawa et al. (2021) recently found generally equal diffusion properties of both isotopes in artificial and natural porous media.

Table 2 Mean and standard deviation of isotopic concentrations in each tension area and time point for simulation with constant (const.) and pore-size-distributed (distr.) D values, respectively. Bold values are within the standard deviation range of the measured values of Bowers et al. (2020).

We use a total of 10⁵ particles, which corresponds to 500 particles per bin at full saturation. A high number of particles is needed to enable a stochastically valid random walk process (cf. Eq. 3; Zehe and Jackisch, 2016). Initially, all particles are distributed randomly over all bins in L_D, and thus, each particle is assigned an exact position and bin number within the pore space prior to the start of the mixing process. To saturate the pore space with the two isotopically distinct waters, we further initially define the particles in each bin that contain the isotopically light or heavy water concentrations. According to the binning of the water retention curve (cf. Sect. 2.2), we identify that bins 168–200 correspond to a relative saturation of 16 %–17 % (cf. Sect. 3.1.1). Thus, these bins are filled with particles carrying the light water concentration, mimicking the initial water content stored at high matric potentials in the smallest pores. The residual bins 1–167 are filled accordingly with particles carrying the heavy water concentration representing new input water. Furthermore, we also link the three tension areas to bin numbers, with bins 1–143 as low-tension areas, bins 144–177 as mid-tension areas, and bins 178–200 as high-tension areas. Isotopic concentrations in these tension areas are calculated by averaging the concentrations of all particles present in the corresponding bin numbers.

3.2 Set-up of virtual experiment: simulating diffusive pore mixing and leaching of solute during steady-state saturated flow in a soil column

The virtual experiment serves as an additional evaluation of the capability of the DIPMI approach to simulate pore scale mixing in a more complex setting. We use the term virtual to emphasize that these numerical simulations do not rely on real, existing experiments, which is in contrast to the experiments of Bowers et al. (2020) in the first part of this study.

For the virtual experiment, we assume a vertical 1-D soil column of length z=1.0 m, which is subdivided into vertical layers of dz=0.1 m length (cf. the top of Fig. 1). The (fully water-saturated) soil column contains the same macroscopically homogeneous sandy loam with a saturated hydraulic conductivity K_s of 10⁻⁶ m s⁻¹ and has all other hydraulic properties and the definition of the three tension areas in each soil layer, as used in the experiment of Bowers et al. (2020). All other experimental and simulation parameters are also the same (cf. Table 1). Water particles initially located in the pore space of the surface soil layer (0–0.1 m) carry a concentration C=100 M L⁻³ of a representative conservative solute, while water particles in the other (lower) soil layers carry a zero solute concentration. The soil column is then irrigated by pure water without any solute. A steady-state flow through the soil domain is established for 7 d, driven by a free drainage condition at the bottom boundary, neglecting any evaporation effects at the soil–atmosphere interface. For the virtual experiment, the vertical displacement routine of LAST is used, assuming pure matrix flow without the influence of macropores or reactive transport processes (cf. Sect. 2.1, Zehe and Jackisch, 2016; Sternagel et al., 2021). It calculates a vertical displacement step, by means of advection and dispersion, for each fluid particle in all soil layers, starting from the bottom to the surface layer. Thus, a certain number of particles initially leaves the soil domain via the bottom boundary, and to maintain the saturation state, missing numbers of particles in soil layers are gradually refilled by particles from overlying layers until the soil domain is, at the top, finally re-saturated by adding new event particles to the surface soil layer (steady state). The length of a vertical displacement step is, therefore, also dependent on the bin/pore size (cf. Sect. 2.1). Particles in the coarse pores experience a larger vertical displacement than particles in smaller pores due to higher advective velocities and diminished capillary effects. Hence, particles in coarse pores/bins are more likely to travel into the next underlying soil layer.

We simulate this virtual experiment set-up with our LAST-Model using (i) the DIPMI approach, with constant and pore-size-distributed diffusivities D (cf. Sect. 2.2) over bins, respectively, and (ii) the common perfect mixing assumption used in the former versions of LAST (cf. Sect. 2.1; Sternagel et al., 2019, 2021).

4 Results and discussion

In Sect. 4.1 and 4.2, we present and discuss the results of the DIPMI simulations of the experiment of Bowers et al. (2020), followed by the presentation and discussion of the results of the virtual experiment in Sect. 4.3 and 4.4.

4.1 DIPMI simulations of the experiment of Bowers et al. (2020)

Table 2 contains the mean and standard deviation of isotopic concentrations in each tension area and the time point for each simulation with constant and pore-size-distributed D values, respectively. Values highlighted in bold are within the measured standard deviation range of the experimental values of Bowers et al. (2020), which are additionally given. It is obvious that most of the simulated isotopic values are in accordance with the observation. Deviations mainly occur in the high-tension area after 8 h, as simulations over-predict the degree of mixing (i.e. mixing is too fast) in this high-tension area. All three tension areas require different times to reach a mean equilibrium concentration δ_e of around −54 ‰ for δ²H and −8.0 ‰ for δ¹⁸O. These equilibrium concentrations are reached in the simulation with pore-size-distributed D values after (i) 8 h to 1 d in the low-tension area, (ii) ∼1 d in the mid-tension area, and (iii) 3 d in the high-tension area. Thus, our DIPMI approach simulates complete isotope mixing somewhat faster than Bowers et al. (2020), who found that complete mixing is achieved after around 4 d. However, mixing times and isotopic concentrations of our simulations, with pore-size-distributed D values in particular, are generally consistent with the measured values. Comparing results of simulations with constant and pore-size-distributed D values reveal differences in the mid- and, especially, high-tension areas, while isotopic concentrations in the low-tension area are quite similar. All tension areas reach the equilibrium state already between 8 h and 1 d when simulating with constant D values in all bins.

Additionally, Fig. 2 shows the isotopic concentrations, simulated with pore-size-distributed D values and the observations, for each tension area and time point in a dual-isotope space. The isotopic concentrations in the three tension areas gradually converge towards the mean equilibrium concentrations over time. Simulated concentrations in low- and mid-tension areas are in accordance with the observations, as indicated by the overlapping of simulated and measured standard deviation ranges. An exception is the high-tension area, where measured isotopic concentrations of especially δ²H have lower values compared to our simulations after 8 h; these findings show that the actual mixing process in the smallest pores is delayed. Furthermore, there are also differences between the two isotope species. For δ¹⁸O, simulated concentrations are consistent with the observations in all tension areas over the entire period, except the 7 d concentrations in the mid- and high-tension areas. However, these measured concentration values were unexpectedly high in the experiments of Bowers et al. (2020), as further discussed in Sect. 4.2. The value range of δ¹⁸O is also generally larger compared to δ²H.

Figure 2 Isotopic concentrations in a dual-isotope space for each tension area and time point, simulated (red) by the DIPMI approach with pore-size-distributed D values and measured (black) by Bowers et al. (2020). Note that, at t=0 h, the simulated and measured isotopic values in the low-tension and high-tension areas are identical.

4.2 Analysis of simulations of the experiment of Bowers et al. (2020)

The simulation results (cf. Sect. 4.1) indicate that our diffusive pore mixing (DIPMI) approach is capable of reproducing the experiment of Bowers et al. (2020). The results in Sect. 4.1 show that our approach is suitable for (i) simulating the measured mixing time (∼3 d) required to reach an equilibrium concentration δ_e in the entire pore space and (ii) resolving most of the isotopic concentrations in all tension areas with time (cf. Sect. 4.1, Table 2, and Fig. 2).

In general, the concept of DIPMI is consistent with other studies, which usually observe a lag time in the isotopic mixing of waters stored initially in different pore sizes (e.g. Adams et al., 2020; Bowers et al., 2020; Orlowski and Breuer, 2020), proving that diffusive mixing over an entire pore space is far from being a perfect, instantaneous process, even on the (small) scale of a few centimetre long soil sample. A realistic description of diffusive mixing processes is especially crucial for interpreting studies examining the origins of plant water (e.g. Sprenger et al., 2016; Penna et al., 2020), water ages (e.g. Hrachowitz et al., 2013; Sprenger et al., 2019), travel times (e.g. Klaus et al., 2015; van der Velde et al., 2015), or even groundwater (e.g. Berkowitz et al., 2016).

Our calculated extent of L_D=21 000 µm (cf. Eqs. 1 and 2) appears to feasibly resolve the internal pore space of the sandy loam soil volume used in the experiment of Bowers et al. (2020). The distributed bin-dependent diffusivities D (cf. Eq. 4) facilitate a realistic simulation of measured isotopic concentrations, which is superior to simulations with constant D values (cf. Table 2). The latter results in an over-prediction of the degree of mixing, which reaches the equilibrium state almost simultaneously in all tension areas after just 8 h to 1 d. A constant diffusivity of $\mathrm{2.272}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹ seems feasible only in coarser pores. The smaller a pore size class, the greater the capillary forces and friction caused by interactions between solid surfaces and the thin water layers directly adjacent to them. Both capillary forces and friction significantly decrease the water movement in the smallest pores. This effect corroborates our calculation of pore-size-dependent diffusivities. However, D values in the high-tension area are still too high because the measured concentrations are not matched to the same extent as for the case in the low- and mid-tension area. This implies that our approach still simulates overly strong mixing in the high-tension area.

The high-tension area of the water retention curve probably bears the highest uncertainty. In this area of the pore space, water is held at tensions ≫ 1.0 MPa and is usually extracted by cryogenic vacuum distillation (CVD). This method, however, may be prone to producing artefacts when analysing the isotopic concentrations of soil and plant water (Orlowski et al., 2013; Adams et al., 2020; Orlowski and Breuer, 2020; Allen and Kirchner, 2021). Beside possible uncertainties in the measurement procedure, specific subscale processes can be the reason for discrepancies between simulation and observation in the high-tension area. Bowers et al. (2020) suggested that the interactions/adsorption of water ions with clay minerals within smaller pores can have a significant effect on the mixing behaviour. This might also be the reason for the discrepancies between simulated and measured δ¹⁸O concentrations in mid- and high-tension areas after 3–7 d, as measured values are higher than expected, as stated by Bowers et al. (2020). Such δ¹⁸O enrichment was also reported in other studies (e.g. Oerter et al., 2014). Orlowski and Breuer (2020) further suggested that isotopic fractionation may occur during diffusive mixing, especially in high-tension areas. Reasons for such isotopic fractionations are difficult to define, and in addition to adsorption effects, further subscale processes may play a role, such as adhesion or evaporation effects. Evaporation is often regarded as a main driver for isotopic fractionation, especially at the interface of the soil–atmosphere system (Sprenger et al., 2016, 2018). However, it may also result in a fractionation effect on the pore scale during the water extraction process in experimental studies when a phase change, from liquid to gaseous, occurs at the interface of saturated and unsaturated pores, which in turn may lead to the vapour-pressure-controlled adsorption of water on soil surfaces (Lin and Horita, 2016; Lin et al., 2018).

Such detailed subscale processes are not incorporated specifically in the current version of our DIPMI approach. We focus on the abstraction of physical properties, which we think generally have the greatest influence on the diffusion process inside of soil domains, e.g. the diffusion coefficient of certain fluid or substance, distinct water pockets in soil pores of different sizes, and tortuosity of flow paths (Chou et al., 2012; Bowers et al., 2020). We lump these physical properties into the two main assets of our DIPMI approach, i.e. (i) the extent and characteristic flow length of the pore space L_D and the (ii) variable diffusivities D in different pore size classes, both of which can be directly derived from the soil water retention curve (cf. Sect. 2.2).

The results of Bowers et al. (2020) and our simulations both highlight that mixing processes in soils are by no means instantaneous or perfect, even in very small and fully saturated control volumes. Rather, the diffusive spreading of water depends on the pore size distribution and specific soil water retention properties. With this insight, it is of interest to examine, in the following virtual experiment (Sect. 4.3 and 4.4), how pore-size-dependent and non-instantaneous mixing affect simulations of water flow through a saturated soil column on a larger scale and to delineate their effects on solute breakthrough and redistribution in soil.

4.3 Simulation of the virtual experiment

Solute breakthrough curves exhibit a clear difference between simulations with the DIPMI approach and the perfect mixing assumption (Fig. 3). A simulation with the perfect mixing assumption results in a breakthrough curve that is shaped like an approximately normal distribution, with a concentration peak of 13 M L⁻³ after ∼42 h, followed by a sharp decrease converging in a zero solute concentration after ∼105 h. Thus, all solute stored initially in the surface soil layer is eluted completely from the entire soil column by ∼2.3 pore volumes when the simulation is performed with the perfect mixing assumption. In contrast, the simulation with the DIPMI approach and constant D values results in a breakthrough curve that exhibits a right-skewed distribution. The peak concentration of ∼7 M L⁻³ is reached after around 38 h and is followed by a long tailing of solute breakthrough, which converges the zero concentration only after 7 d, corresponding to 3.5 pore volumes. Hence, significantly more time is required to elute all solute from the soil column. The simulation with the DIPMI approach and pore-size-distributed D values generally results in a breakthrough curve with a similar shape, i.e. as the breakthrough curve with constant D values, with comparable peak concentration and long tailing behaviour. However, it needs a shorter time to peak (∼22 h), and the solute breakthrough tailing does not converge to zero concentration at all, not even after the simulation period of 7 d. Thus, more than 3.5 pore volumes are needed to leach all solute from soil in the case of pore-size-distributed D values.

Figure 3 Solute breakthrough curves over 7 d, simulated with (i) the DIPMI approach and constant D (const.) values (blue), (ii) DIPMI approach and pore-size-distributed D (distr.) values (red), and (iii) the perfect mixing assumption (black).

Figure 4 further shows the mean solute concentrations in the three tension areas per vertical soil layer at different points in time. Comparing the results of simulations with the DIPMI approach (with pore-size-distributed D values) and the perfect mixing assumption generally supports the previous insights of the breakthrough curves.

Figure 4 Vertical concentration profiles of a solute with mean concentrations in three tension areas over 7 d, simulated by the DIPMI approach with pore-size-distributed D values (red shades) and the perfect mixing assumption (black), respectively. Note that the black line and crosses illustrate the same mean concentration in all three tension areas per soil layer, as the perfect mixing assumption averages out pore-size-dependent differences.

There is no difference in concentration between tension areas when using the perfect mixing assumption, as this approach averages solute concentrations over the entire pore space of a soil layer (cf. Sect. 2.1). Consequently, the solute gradually propagates, in the form of a fast wave, through all soil layers, resulting in a complete elution of solute in all soil layers within the first 3–4 d. The results of the simulation with the DIPMI approach and pore-size-distributed D values exhibit a different and more heterogeneous picture (red shaded symbols and lines), especially regarding solute propagation behaviour and concentrations in different tension areas. Due to the imperfect mixing, a large fraction of solute preferentially propagates vertically through the low-tension areas (dark red) of the soil layers. This vertical leaching process is faster than the diffusive mixing over the tension areas within a soil layer, leading to the fast breakthrough of solute during the first 12 h (cf. Fig. 3). Hence, only smaller amounts of solute enter the mid- and high-tension areas (light red and orange) in deeper soil layers (0.1–1.0 m). In contrast, the initially solute-filled mid- and high-tension areas in the surface soil layer only release their solute slowly. This mechanism entails a retardation effect, which results in the persistence of solute in soil over the entire period of 7 d and the long breakthrough curve tailing with incomplete elution of solute (cf. Fig. 3).

4.4 Analysis of the virtual experiment

The use of the perfect mixing assumption shows two effects in our virtual experiment, namely (i) a longer time to the first breakthrough and peak (cf. Fig. 3) and (ii) a steeper and shorter tailing of the solute breakthrough concentration after the peak with a complete leaching of solute within 105 h (cf. Fig. 4). These effects can be explained by the fundamental assumption that solutes always perfectly and instantaneously mix over the entire pore space in each soil layer during passage of the soil domain. Consequently, solute spreads uniformly over the pore spaces of initially solute-free soil layers (0.1–1.0 m) before reaching the bottom boundary, which leads to dilution and the delayed breakthrough front arrival. Thereafter, the perfect and instantaneous mixing of pure infiltrating event water with solute-containing, pre-event water causes the rapid and complete elution of all solute, especially out of the surface soil layer, and also solute with a relatively high peak concentration. These dilution and elution effects are characteristic of the over-mixing phenomena (e.g. Green et al., 2002; Boso et al., 2013) and may be problematic for assessing the risk of contaminant leaching by potentially giving wrong predictions regarding breakthrough times and persistence in soil, for example.

Using the DIPMI approach with pore-size-distributed D values (cf. Eq. 4) for the simulation of the virtual experiment bears opposite effects, such as (i) a faster initial breakthrough (cf. Fig. 3) and (ii) a long tailing of solute breakthrough concentration with an incomplete elution of solute (cf. Fig. 4). The timescale for diffusive mixing over pore sizes is significantly larger than the timescale for vertical transport of solute. This leads to an early arrival of the breakthrough front, as solute travels downward mainly through the pores of low-tension areas without spreading uniformly over soil layers. The smaller peak concentration and the long tailing of the breakthrough curve are, thereafter, caused by a retardation effect of pores in the mid- and high-tension areas. These pores, with smaller diffusivities and higher capillary tensions, bind the solutes for a longer time (cf. Fig. 4) before they are mixed diffusively with pure infiltrating water and leached out. Regarding these effects, it is important to recall that the vertical displacement of fluid particles also depends on the different pore size classes/bins in our LAST-Model (cf. Sect. 2.1; Sternagel et al., 2019, 2021). Particles in coarse pores of low-tension areas are more mobile and displaced vertically by a higher advective velocity compared to particles in smaller pores. Consequently, water and solute flow mainly through these pores in low-tension areas, with a limited diffusive exchange with smaller pores, and thus, the mid- and high-tension areas are essentially bypassed. This behaviour implies that the saturated flow system is dominated by preferential bypassing flow through the low-tension area.

In the low-tension area, the mean D values of the constant D method ($\mathrm{2.272}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹) and the pore-size-distributed D method (1.243$\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹) are (relatively) similar, leading to breakthrough behaviours of both methods that are comparable in an early phase of breakthrough (cf. Fig. 3). However, the overall higher and constant D value causes a stronger leaching of retarded solute over time in later phases of breakthrough, especially from mid- and high-tension areas, with a complete elution after the period of 7 d. This implies that, during early phases of fast bulk leaching of solute, the influence of pore-size-dependent diffusive mixing is less significant due to preferential bypassing. Nevertheless, pore-size-dependent diffusive mixing becomes highly relevant in later phases of the breakthrough when residual amounts of solute, stored and retarded in small pores, are gradually moved back into coarser pores by diffusive mixing with infiltrating water and hence, remobilized. Note that we perform our simulations in saturated media because, in a fully saturated pore space, differences of the diffusivities between the largest saturated pore size class and the smallest saturated pore size class are more distinct than in a partially saturated pore space. Thus, simulations under saturated conditions are more suitable for highlighting the influence of diffusive mixing with pore-size-distributed D values (also in comparison to the use of constant D values and the perfect mixing assumption).

The breakthrough curves (cf. Fig. 3) also show the sensitivity of the DIPMI approach to variations in the diffusion coefficient value. The simulations are performed with significantly different D values in the set-ups with distributed (i.e. mean D of $\mathrm{9.7}\times {\mathrm{10}}^{-\mathrm{10}}$ m² s⁻¹) and constant (i.e. mean D of $\mathrm{2.272}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹) D values (i.e. range of 57 %), as well as with the instantaneous, perfect mixing assumption (i.e. infinite D of 10^∞ m² s⁻¹). We can infer that (i) higher D values result in breakthrough curves with shorter tailings, which gradually approach the shape of the curve of the perfect mixing assumption with higher D values, and (ii) smaller D values result in breakthrough curves with increasingly longer tailings. The latter situation is the case when simulating the diffusive spreading of common solute tracers (e.g. bromide and chloride) in porous media because diffusion coefficients of the dissolved solute molecules are usually smaller compared to those of pure water or stable water isotopes (Hasegawa et al., 2021).

Overall, the results of the virtual experiment reveal a different behaviour in the solute leaching and redistribution in soil for the simulation with the DIPMI approach, compared to simulations that invoke the perfect mixing assumption. In particular, the long tailing of breakthrough curves is in line with common observations of several studies (e.g. Zinn and Harvey, 2003; Willmann et al., 2008; Edery et al., 2014). Zinn and Harvey (2003) linked the long tailing of breakthrough curves to mass transfer between regions of different mobility, e.g. pore size classes in different tension areas. Edery et al. (2014) also simulated the solute breakthrough in saturated, porous media by a Lagrangian particle tracking approach. They showed that the broadness and tailing of breakthrough curves increase with a generally higher heterogeneity of pore space and flow paths. Our subdivision of the pore space into different pore size classes with pore-size-dependent diffusivities in the DIPMI approach is in line with this finding and adds, furthermore, an important aspect to account for imperfect subscale mixing in soil hydrological modelling.

Moreover, our results highlight that these typical shapes of breakthrough curves are not exclusively caused by explicit hydraulic (e.g. macropore flow) or chemical (e.g. adsorption and desorption) heterogeneity in soil, but that early breakthroughs and long tailing may also be a result of imperfect diffusive mixing within the pore space (Willmann et al., 2008), even when the flow domain is fully saturated and its soil properties are macroscopically homogeneous. Hence, we can verify our proposed hypothesis (cf. Sect. 1) and state that imperfect mixing across soil pores with different hydraulic properties within the soil matrix may entail an implicit heterogeneity of flow, causing typical shapes of breakthrough curves, even in the absence of explicitly defined preferential flow paths.

5 Conclusions and outlook

The main findings of this study are as follows:

• Simulation results show the feasibility of the DIPMI approach for reproducing the findings of the experiment of Bowers et al. (2020) by correctly capturing the measured diffusive mixing times and concentrations of water isotopes over the pore space of a fully saturated soil sample (cf. Sect. 4.1 and 4.2).

• A virtual experiment highlights that imperfect mixing in a macroscopically homogeneous soil matrix can produce breakthrough patterns that are usually associated with preferential flow in heterogeneous soil structures such as macropores (cf. Sect. 4.3 and 4.4).

• The DIPMI approach is physically based in the sense that its model parameters can be derived from the soil water retention curve, and no further calibration is needed (cf. Sect. 2.2).

In the future, in situ 1-D column experiments are planned to analyse the influence of the microstructure of partially saturated soils on the temporal and spatial mixing of isotopes over the pore space. The experimental results may provide a representative dataset serving as a reference to further extend and test the DIPMI approach, as comparable data on pore-scale mixing processes remain scarce.