# Influence of Clay on Solute Transport in Saturated Homogeneous Mixed Media

Authors: Albert Kwame Kwaw; Zhi Dou; Jinguo Wang; Yuting Zhang; Xueyi Zhang; Wenyuan Zhu; Portia Annabelle Opoku
Journal: Geofluids
 (2021)

---

## Abstract

In this study, four homogeneous porous media (HPM1-HPM4), consisting of distinct proportions of sand-sized and clay-sized solid beads, were prepared and used as single fracture infills. Flow and nonreactive solute transport experiments in HPM1-HPM4 under three flow rates were conducted, and the measured breakthrough curves (BTCs) were quantified using conventional advection-dispersion equation (ADE), mobile-immobile model (MIM), and continuous time random walk (CTRW) model with truncated power law transition time distribution. The measured BTCs showed stronger non-Fickian behaviour in HPM2-HPM4 (which had clay) than in HPM1 (which had no clay), implying that clay enhanced the non-Fickian transport. As the fraction of clay increased, the global error of ADE fits also increased, affirming the inefficiency of ADE in capturing the clay-induced non-Fickian behaviour. MIM and CTRW performed better in capturing the non-Fickian behaviour. Nonetheless, CTRW’s performance was robust. 12.5% and 25% of clay in HPM2 and HPM3, respectively, decreased the flowing fluid region and increased the solute exchange rate between the flowing and stagnant fluid regions in MIM. For CTRW, the power law exponent (βCTRW) values were 1.96, 1.75, and 1.63 in HPM1-HPM3, respectively, implying enhanced non-Fickian behaviour. However, for HPM4, whose clay fraction was 50%, the βCTRW value was 1.87, implying a deviation in the trend of non-Fickian enhancement with increasing clay fraction. This deviation indicated that non-Fickian behaviour enhancement depended on the fraction of clay present. Moreover, increasing flow rate enhanced the non-Fickian transport based on βCTRW.


---

## Body

## 1. Introduction

Filled fractures are ubiquitous in the subsurface. Despite the narrow nature of such fractures compared to the size of the host rock, the fracture infills, most often sediments, are regarded as porous media through which fluid flow and solute transport occur. The sediment infills are in different forms ranging from purely sand to complex compositions of sand, silt, and clay. Previous studies have shown that sediments that fill the void spaces in single fractures have significant implications for flow and transport [1, 2]. Despite numerous studies, the role of clay, as a proportion of homogeneous mixed sediment infill, on the transport of solutes in single fractures is poorly understood. Below are extensive researches about the significance of clay on dissolved contaminant transport in natural and engineered systems.At Sarnia in Ontario, Goodall and Quigley [3] detected pollution fronts beneath two young landfills built on dense silty clay formations. Crooks and Quigley [4] performed laboratory and field studies and compared the results of brine leachate through clay. In these studies, diffusion was identified as the main mechanism responsible for the migration of the dissolved chemical substances through the clays. Apart from diffusion, which is a physical process, clay also attenuates the transport of dissolved contaminants through chemical processes, for example, complexation and adsorption [5–8], ultimately causing the groundwater cleanup time to be longer than expected [9–11].To further understand the influence of clay on the contaminant transport behaviour in aquifer systems, researchers have paid much attention to scenarios where clays exist as lenses within aquifers [12, 13] and aquitards around aquifers [9, 14–20] in contact with dissolved plumes. These studies documented that clays occurring as lenses within aquifers and as aquitards near aquifers inhibited the transport of contaminants, yielding comparatively low concentrations of plumes at late times. In the context of groundwater quality, the occurrence of such low concentrations at times longer than expected is a problem and needs critical attention.Solute transport modelling is essential. Nonetheless, the question of which model can quantify the migration of solutes in porous media remains challenging. The conventional advection-dispersion equation (ADE) model has been proposed to be adequate in modelling solute migration through homogeneous porous media [21]. However, in homogeneous porous media, there exist inherent heterogeneities [22, 23] and this gives breakthrough curves (BTCs) characterised by early arrival times and long late temporal tails. Hence, using the ADE solution to model the travel of solutes in such media may not be sufficient.A mobile-immobile model (MIM), an improved version of the ADE model, was developed to resolve the nonequilibrium processes causing the anomalous transport behaviour of passive solutes. In the MIM model, long late-time tailing BTCs are ascribed to the trapping of portions of solutes in stagnant fluid regions. The transfer of solute between flowing and stagnant fluid regions is called the physical nonequilibrium process [24, 25]. Even though the MIM model is widely employed to model solute transport through homogeneous porous media [26–28], its validity is questionable in the sense that the transfer of solutes between the flowing and stagnant fluid regions, in reality, is a multiple-rate phenomenon and not a single rate suggested by the MIM model.To address the issue of anomalous transport, the continuous time random walk framework (CTRW), a type of time nonlocal transport model, has been developed and widely used in modelling solute transport through different porous media. In homogeneous porous media, CTRW has been compared with MIM and/or ADE in many studies. For example, in a sand medium, Levy and Berkowitz [29] observed that the CTRW model performed better than the ADE model in terms of capturing the quick arrival time and long late temporal tailed BTCs. An analogous finding was reported by Cortis et al. [30] when they explored the transport behaviour of a conservative tracer through homogeneous porous media and fitted the observed BTCs with both ADE and CTRW models. Bromly and Hinz [31] conducted repacked sand experiments and discovered that the MIM model could not fit the measured BTCs, especially the long late temporal tails, while the CTRW model did. Zhang and Lv [32] also reported similar diffusive-flux of conservative solute in homogeneous media. Recently, Zaheer et al. [33] conducted solute displacement experiments in clay and documented that the CTRW model performs better in accounting for early breakthrough time and long late temporal tails than the ADE and MIM models. This implies that it is insufficient to use the classical ADE model and its variants to describe the transport behaviour of a solute through homogeneous porous media.In order to choose a suitable model to quantify the transport of solutes in a sediment-filled fracture, the mechanism driving the solute in that fracture infill must be understood. However, in the case of homogeneous mixed media (with distinct fractions of clay and sand) filling a single fracture, the knowledge of solute transport behaviour is limited. To create well-controlled experimental conditions and have more access to the flow and transport data than in the field, we used an idealised one-dimensional physical model to represent the single fracture. It should be mentioned that the idealised single fracture employed in this study may not be the case in reality since natural filled single fractures have nonuniform apertures with rough surfaces. The objectives of this study are to (i) explore the effect of clay on solute transport characteristics in homogeneous mixed media as a single fracture infill, (ii) find an appropriate transport model that can describe the solute transport behaviour, and (iii) assess how distinct flow rates can impact the solute transport behaviour.


## 2. Materials and Methods

### 2.1. Experimental Set-Up and Materials

The experimental set-up (Figure1) comprised a deionised water tank, a tracer tank, a peristaltic pump (BT101L+YZ15, Lead Fluid, China), flow control valves, connecting tubes, a column, two digital pressure gauges (HD-100G, China), and an automatic effluent collector (SBS-100, Huxi, China). The two tanks were joined to the column through the peristaltic pump using tubes and flow control valves. The peristaltic pump, whose flow control accuracy is ±0.5%, has a head pump flow ranging from 0.006 to 2300 mL/min. Valve A changed the flow from the deionised water tank to the tracer tank to inject the solute and vice versa to flush the column. The column, made of acrylic glass, served as the idealised single fracture. The column was 51 cm in length with an inner diameter and a thickness of 5 cm and 1 cm, respectively, with both ends sealed by flanges and shims to prevent water leakage. Both ends of the column had three metal pipes purposely to get a uniform distribution of flow in the column. The pressure gauge measured the pressure at the inlet and outlet of the column.Figure 1
The experimental setup.To avoid the effect of chemical reaction in the solute transport process, we used synthetic glass beads with no intragranular pore space to represent the porous media. Four homogeneous porous media (HPM1-HPM4), comprising different fractions of sand-sized and clay-sized particles, were prepared and used as single fracture infills. The sand-sized and clay-sized beads were spherical. The estimated hydraulic conductivities of high- and low-permeability beads were1.94×10−2°cm/s and 3.05×10−5°cm/s. The range of sand and clay particle sizes, the proportions of sand and clay, and the porosity of the porous media are summarised in Table 1. Brilliant Blue FCF, commonly used as a colourant in foods and medications with a maximum optical absorption wavelength of 630 nm [34], was used as the nonreactive tracer.Table 1
Porous media properties.
Porous mediaParticle size range of sand and clay (mm)Mass proportions of sand and clay (%)Porosity (n)HPM10.355-0.51000.368HPM20.355-0.5, 0.071-0.0987.5 : 12.50.316HPM30.355-0.5, 0.071-0.0975 : 250.285HPM40.355-0.5, 0.071-0.0950 : 500.230


### 2.2. Experimental Procedure

For each column experiment, the beads were washed and soaked in deionised water for 24 h purposely to remove any impurities on the surface of the beads. In filling the column, water was first poured into the column placed in a vertical position followed by the wet beads. Small quantities of the wet beads were poured and pressed firmly with a rod, and this was to obtain dense uniform media. Water was always above the beads throughout the filling, and the reason was to prevent air bubbles in the column. The column was laid horizontally on two metal supports and connected to the feed tanks, pump, and pressure gauges using the tubes and valves. Filters were placed at both ends of the column to prevent the beads from flowing out of the column. Deionised water was injected into the column until steady-state saturated flow was achieved.An accurate description of a flow field is essential for studying the transport behaviour of solutes in porous media. In this study, flow experiments, with different pump rates (0.48, 0.72, 0.96, 1.20, 1.44, and 1.68 mL/min), were performed. To verify Darcian flow, since the groundwater flow is limited to Darcy’s law, the hydraulic gradient (J) and specific discharge (q=Q/A) were estimated, where Q and A are the discharge rate (L3T-1) and inner cross-sectional area (L2) of the column, respectively. The hydraulic gradient (J) was computed using
(1)h=Pρwg,(2)J=dhl,where h is the hydraulic head, p is the inlet or outlet pressure (ML-1 T-2), ρw is the density of water (ML-3), g denotes the acceleration due to gravity (LT-2), dh is the difference between the inlet and outlet hydraulic heads, and l is the length of the column (L). Out of the six tested pump rates, three (i.e., 0.96, 1.44, and 1.68 mL/min) were used in the solute displacement experiments. Table 2 summarises the discharge rates and average flow velocities for the three pump rates in HPM1-HPM4.Table 2
Flow conditions in HPM1-HPM4 at 0.96, 1.44, and 1.68 mL/min.
Porous mediaPump rate (mL/min)Discharge rateQ (mL/min)Average flow velocityv¯∗ (cm/min)HPM10.960.79950.1111.441.29910.1801.681.57930.219HPM20.960.89980.1451.441.34920.2181.681.61890.261HPM30.960.89970.1611.441.23260.2201.681.49900.268HPM40.960.84960.1881.441.24920.2771.681.49910.332v¯∗=Q/An, where n and A are the porosity and cross-sectional area of the column, respectively. A is 19.63 cm2.A dimensionless number, Reynolds number (Re), was employed to assess the flow regime in HPM1-HPM4. The Reynolds number is simply a ratio of inertial force to viscous force and is given as(3)Re=ρqdμ,where ρ, q, d, and μ denote the density of the fluid (ML-3), specific discharge (LT-1), the characteristic length of flow (L), and the dynamic viscosity of the fluid (ML-1 T-1), respectively. The properties of water at 25°C (i.e., ρ=997.04 kg/m3, μ=0.89×10-3 Pa·s) and d=5 cm were used. In this study, d is the diameter of the column, representing the aperture of the filled single fracture. The aperture of a filled single fracture has also been used in a previous study by Qian et al. [1]. In their study, they documented that using the average grain size as the characteristic length may not be always convenient when the size of grains varies significantly, which was exactly the case in this study. The estimated Re values were less than some values within the range 1 and 10, signifying the Darcian flow regime [21] in HPM1-HPM4 for the six tested pump rates. Only Re values for 0.96, 1.44, and 1.68 mL/min are shown in Table 3. Moreover, the relation between hydraulic gradient and flow velocity was linear with the coefficient of determination (R2) above 0.990 as shown in Figures 2(a)–2(d). The linear relationship implied Darcian flow in HPM1-HPM4. The verification of Darcian flow regime was done to ensure that there was no chaotic flow phenomenon in the packed media since such flow can induce a non-Fickian transport behaviour [35] even in homogeneous porous media [36].Table 3
Estimated Reynolds number (Re) for HPM1-HPM4 at 0.96, 1.44, and 1.68 mL/min.
Pump rate (mL/min)HPM1HPM2HPM3HPM4ReReReRe0.960.3800.4280.4280.4041.440.6180.6420.5860.5941.680.7510.7700.7130.713Figure 2
Relationship between hydraulic gradient and specific discharge for HPM1 (a), HPM2 (b), HPM3 (c), and HPM4 (d).
(a)(b)(c)(d)After the verification of Darcian flow, a Brilliant Blue FCF solution, with a concentration of 4 mg/L, was injected continuously into the column, representing a step input boundary condition. The tracer injection continued until the effluent concentration was uniform and equivalent to the input concentration (C0). The 4 mg/L solution was replaced with deionised water for the flushing of the tracer. The tracer injection and flushing were carried out using 0.96, 1.44, and 1.68 mL/min at a constant room temperature of 25°C. At the outlet of the column, effluents were collected in 10 mL-sized glass tubes placed in an automatic part collector (SBS-100, Huxi, China) over the entire period of the experiment. An ultraviolet-visible light spectrophotometer (HACH DR5000) was used to measure the absorbance (Ab) of the effluents at a wavelength of 630 nm. At this wavelength, a concentration above 4 mg/L was not detected by the spectrophotometer, hence the choice of 4 mg/L as the initial concentration in this study. The absorbance (Ab) values were then converted to concentrations (C) using the relation Ab=0.1397C–0.0006, with the coefficient of determination (R2) above 0.999.


## 2.1. Experimental Set-Up and Materials

The experimental set-up (Figure1) comprised a deionised water tank, a tracer tank, a peristaltic pump (BT101L+YZ15, Lead Fluid, China), flow control valves, connecting tubes, a column, two digital pressure gauges (HD-100G, China), and an automatic effluent collector (SBS-100, Huxi, China). The two tanks were joined to the column through the peristaltic pump using tubes and flow control valves. The peristaltic pump, whose flow control accuracy is ±0.5%, has a head pump flow ranging from 0.006 to 2300 mL/min. Valve A changed the flow from the deionised water tank to the tracer tank to inject the solute and vice versa to flush the column. The column, made of acrylic glass, served as the idealised single fracture. The column was 51 cm in length with an inner diameter and a thickness of 5 cm and 1 cm, respectively, with both ends sealed by flanges and shims to prevent water leakage. Both ends of the column had three metal pipes purposely to get a uniform distribution of flow in the column. The pressure gauge measured the pressure at the inlet and outlet of the column.Figure 1
The experimental setup.To avoid the effect of chemical reaction in the solute transport process, we used synthetic glass beads with no intragranular pore space to represent the porous media. Four homogeneous porous media (HPM1-HPM4), comprising different fractions of sand-sized and clay-sized particles, were prepared and used as single fracture infills. The sand-sized and clay-sized beads were spherical. The estimated hydraulic conductivities of high- and low-permeability beads were1.94×10−2°cm/s and 3.05×10−5°cm/s. The range of sand and clay particle sizes, the proportions of sand and clay, and the porosity of the porous media are summarised in Table 1. Brilliant Blue FCF, commonly used as a colourant in foods and medications with a maximum optical absorption wavelength of 630 nm [34], was used as the nonreactive tracer.Table 1
Porous media properties.
Porous mediaParticle size range of sand and clay (mm)Mass proportions of sand and clay (%)Porosity (n)HPM10.355-0.51000.368HPM20.355-0.5, 0.071-0.0987.5 : 12.50.316HPM30.355-0.5, 0.071-0.0975 : 250.285HPM40.355-0.5, 0.071-0.0950 : 500.230


## 2.2. Experimental Procedure

For each column experiment, the beads were washed and soaked in deionised water for 24 h purposely to remove any impurities on the surface of the beads. In filling the column, water was first poured into the column placed in a vertical position followed by the wet beads. Small quantities of the wet beads were poured and pressed firmly with a rod, and this was to obtain dense uniform media. Water was always above the beads throughout the filling, and the reason was to prevent air bubbles in the column. The column was laid horizontally on two metal supports and connected to the feed tanks, pump, and pressure gauges using the tubes and valves. Filters were placed at both ends of the column to prevent the beads from flowing out of the column. Deionised water was injected into the column until steady-state saturated flow was achieved.An accurate description of a flow field is essential for studying the transport behaviour of solutes in porous media. In this study, flow experiments, with different pump rates (0.48, 0.72, 0.96, 1.20, 1.44, and 1.68 mL/min), were performed. To verify Darcian flow, since the groundwater flow is limited to Darcy’s law, the hydraulic gradient (J) and specific discharge (q=Q/A) were estimated, where Q and A are the discharge rate (L3T-1) and inner cross-sectional area (L2) of the column, respectively. The hydraulic gradient (J) was computed using
(1)h=Pρwg,(2)J=dhl,where h is the hydraulic head, p is the inlet or outlet pressure (ML-1 T-2), ρw is the density of water (ML-3), g denotes the acceleration due to gravity (LT-2), dh is the difference between the inlet and outlet hydraulic heads, and l is the length of the column (L). Out of the six tested pump rates, three (i.e., 0.96, 1.44, and 1.68 mL/min) were used in the solute displacement experiments. Table 2 summarises the discharge rates and average flow velocities for the three pump rates in HPM1-HPM4.Table 2
Flow conditions in HPM1-HPM4 at 0.96, 1.44, and 1.68 mL/min.
Porous mediaPump rate (mL/min)Discharge rateQ (mL/min)Average flow velocityv¯∗ (cm/min)HPM10.960.79950.1111.441.29910.1801.681.57930.219HPM20.960.89980.1451.441.34920.2181.681.61890.261HPM30.960.89970.1611.441.23260.2201.681.49900.268HPM40.960.84960.1881.441.24920.2771.681.49910.332v¯∗=Q/An, where n and A are the porosity and cross-sectional area of the column, respectively. A is 19.63 cm2.A dimensionless number, Reynolds number (Re), was employed to assess the flow regime in HPM1-HPM4. The Reynolds number is simply a ratio of inertial force to viscous force and is given as(3)Re=ρqdμ,where ρ, q, d, and μ denote the density of the fluid (ML-3), specific discharge (LT-1), the characteristic length of flow (L), and the dynamic viscosity of the fluid (ML-1 T-1), respectively. The properties of water at 25°C (i.e., ρ=997.04 kg/m3, μ=0.89×10-3 Pa·s) and d=5 cm were used. In this study, d is the diameter of the column, representing the aperture of the filled single fracture. The aperture of a filled single fracture has also been used in a previous study by Qian et al. [1]. In their study, they documented that using the average grain size as the characteristic length may not be always convenient when the size of grains varies significantly, which was exactly the case in this study. The estimated Re values were less than some values within the range 1 and 10, signifying the Darcian flow regime [21] in HPM1-HPM4 for the six tested pump rates. Only Re values for 0.96, 1.44, and 1.68 mL/min are shown in Table 3. Moreover, the relation between hydraulic gradient and flow velocity was linear with the coefficient of determination (R2) above 0.990 as shown in Figures 2(a)–2(d). The linear relationship implied Darcian flow in HPM1-HPM4. The verification of Darcian flow regime was done to ensure that there was no chaotic flow phenomenon in the packed media since such flow can induce a non-Fickian transport behaviour [35] even in homogeneous porous media [36].Table 3
Estimated Reynolds number (Re) for HPM1-HPM4 at 0.96, 1.44, and 1.68 mL/min.
Pump rate (mL/min)HPM1HPM2HPM3HPM4ReReReRe0.960.3800.4280.4280.4041.440.6180.6420.5860.5941.680.7510.7700.7130.713Figure 2
Relationship between hydraulic gradient and specific discharge for HPM1 (a), HPM2 (b), HPM3 (c), and HPM4 (d).
(a)(b)(c)(d)After the verification of Darcian flow, a Brilliant Blue FCF solution, with a concentration of 4 mg/L, was injected continuously into the column, representing a step input boundary condition. The tracer injection continued until the effluent concentration was uniform and equivalent to the input concentration (C0). The 4 mg/L solution was replaced with deionised water for the flushing of the tracer. The tracer injection and flushing were carried out using 0.96, 1.44, and 1.68 mL/min at a constant room temperature of 25°C. At the outlet of the column, effluents were collected in 10 mL-sized glass tubes placed in an automatic part collector (SBS-100, Huxi, China) over the entire period of the experiment. An ultraviolet-visible light spectrophotometer (HACH DR5000) was used to measure the absorbance (Ab) of the effluents at a wavelength of 630 nm. At this wavelength, a concentration above 4 mg/L was not detected by the spectrophotometer, hence the choice of 4 mg/L as the initial concentration in this study. The absorbance (Ab) values were then converted to concentrations (C) using the relation Ab=0.1397C–0.0006, with the coefficient of determination (R2) above 0.999.


## 3. Inverse Models for Solute Transport

### 3.1. ADE Model

The ADE model is based on two assumptions. First, the ADE model treats the transition of solute particles as a Markovian process, meaning the present trait of the solute not dependent on its history [37]. Second, the ADE model adopts a Fickian-driven dispersive flux of a solute, with the coefficient of dispersion assumed to be constant in spatial and temporal scales [38]. For one-dimensional nonreactive solute transport through a homogeneous, isotropic, steady flow porous medium, ADE is given by [28]
(4)∂c∂t=DADE∂2c∂x2−v¯∂c∂x,where C is the solute concentration (ML-3), v¯ is the average pore-water velocity (LT-1), and Dx is the coefficient of hydrodynamic dispersion (L2T-1), t is time (T), and x is a distance (L). The analytical solution of ADE for flux type upstream boundary condition is given as [39]
(5)cx,t=c02erfcx−v¯t2DADEt+c0v¯2tπDADEexp−x−v¯t24DADEt−c021+v¯xDADE+v¯2tDADEexpv¯xDADEerfcx+v¯t2DADEt,where erfc means the complementary error function.


### 3.2. MIM Model

For the MIM model, which describes the physical nonequilibrium part applied in this study, the transport of a solute plume is modelled on the assumption that the porous media through which the plume migrates and has two different liquid phases: mobile (flowing) and immobile (stagnant). The migration of a solute in the flowing region is driven by an advection-dispersion process, while the immobile region exchanges its solute content with the mobile region through a rate-limited diffusion process and is described as a first-order process [25, 39]. For one-dimensional transport of a solute in the absence of chemical reaction, MIM is given as
(6)θm∂cm∂t+θim∂cim∂t=θmDm∂2cm∂x2−θmvm∂cm∂x,(7)ωcm−cim=θim∂cim∂t,where θm, θim, cm, and cim, respectively, denote the amounts of water and concentration in the flowing and stagnant liquid regions; vm and Dm, respectively, define the mean flow velocity and coefficient of hydrodynamic dispersion in the flowing fluid region; and ω is the first-order coefficient of mass transfer between the two flowing and stagnant liquid regions (T-1). The flowing liquid fraction is expressed by a dimensionless β=θm/θm+θim.


### 3.3. CTRW Model

The idea of CTRW was introduced in electronics to explore transport in semiconductors [40] and later applied in geological media to study solute transport [41]. In porous media, a solute is transported by a fluid in which it dissolves through tortuous paths. This transport process is conceptualised by CTRW as particles randomly undergoing a series of “hops” and each “hop” is described by a displacement and waiting time. Here, the displacement of the particles is defined by the physical and chemical heterogeneous nature of the media, and the waiting times describe the durations of hopping events. The displacements and the waiting times are independent random variables, and their distributions are governed by a joint probability density function (PDF), represented by ψs,t. For simplicity, the PDF is decoupled as ψs,t=Psψt, where Ps defines the probability distribution of the distance travelled by the particles between sites and ψt is the probability rate for the transition time to cover the transition displacement [42].As a solute migrates in porous media, a portion of the solute is held in immobile regions and later diffuses slowly out of the region to produce late temporal tailed BTCs, a mark of non-Fickian behaviour. The capturing of the waiting times of the solute held in the immobile region is the core of the CTRW model. This means that CTRW places more emphasis on the transition time distribution than the transition length since the residence times of fluids and solutes quantify transport in porous media. In this study, the concept of CTRW is summarised and the detail of the importance of waiting time distribution, which the CTRW model is all about, is documented in the publications [7, 43].Assuming a Gaussian distribution for the transition displacement, the one-dimensional solute concentration defined in Laplace space for a uniform system is expressed as [44]
(8)uc~x,u−c0x=−M~uvψ∂∂xc~x,u−Dψ∂2∂x2c~x,u,where
(9)M~u≡t¯uψ~u1−ψ~u.In the above equations,M~u is a memory function and is governed by the waiting time distribution [7]. M~u accounts for undetectable heterogeneities that exist at all scales, and u denotes the Laplace variable. c~x,u is Laplace transformed ensemble-average, normalised concentration, whereas c0 and x denote the initial concentration and location in space, respectively. The characteristic time, transport velocity, and generalised coefficient of dispersion are defined by t¯, vψ, and Dψ, respectively. The transport velocity, vψ, is not equivalent to the average pore velocity, v, as both are reported to be similar in the classical ADE. Likewise, Dψ is different from D in ADE.Considering a unit constant flux boundary conditionj~=u−1 at the inlet (x=0) of a finite system of length 1 (i.e., x=0,1), a zero concentration gradient boundary condition (i.e., ∂c~/∂x=0) at the outlet (x=1) and initial condition c0x=0, the resident concentration in the Laplace transform is given as
(10)c~u,x≡1uw−1eq+w+1erw+1M~vψ+2u∝ψew/∝ψ+w−1M~vψ−2u∝ψ,where
(11)w=1+4u∝ψM~vψ,(12)q=1+wx2∝ψ,(13)r=2w−xw2∝ψ.The behaviour of solute transport is quantified by the distribution of transition time, which can take on three different expressions, namely, asymptotic, truncated power law, and modified exponential models. In this study, only the truncated power model is discussed. The comprehensive description of the remaining models is documented in a publication by Cortis and Berkowitz [45]. To quantify the temporal regimes of transport in disordered media, characterised by heterogeneity at all scales, the transition time distribution of solute particles is defined by a truncated power law (TPL) [42]:
(14)ψt=t1τ2−βexpτ2−1Г−β,τ2−1−1exp−t/t21+t/t11+β,which transforms to Fickian behaviour at extended times. In Equation (14), τ2=t2/t1 and Г denotes the incomplete Gamma function. t1 denotes the characteristic median transition time and sets the lower limit for the power law transport behaviour such that t>t1. t2 is a cutoff time for the truncated power-law transition time distribution function that corresponds to the largest heterogeneity length scale. For the time regime t1≪t≪t2, ψt behaves in line with power law whereas ψt decreases exponentially for t≫t2 [42]. To find the memory function, TPL in Laplace transform (Equation (15)) [42] is substituted into Equation (9), with t¯=t1.
(15)ψ~u=1+τ2ut1βexpt1uГ−β,τ2−1+t1uГ−β,τ2−1.For the truncated power law transition time distribution, the nature of solute transport (i.e., non-Fickian or Fickian) is dictated by the parameterβ, where β is a constant exponent. The value of β is in three ranges, 0<β<1, 1<β<2, and β>2. The nature of transport for the interval 0<β<1 is extremely non-Fickian with early and long late-temporal tailed BTCs. For 1<β<2, the mean of the plume migrates with steady velocity and yields non-Gaussian-shaped BTCs with long late-time tailings. When β>2, the movement of the dissolved chemical species is already Fickian for t≫t1, and t2 becomes insignificant in the solute migration process. In this case, the solute plume mass centre migrates with the average fluid velocity and yields a dispersion coefficient that is constant spatially and temporally [42].


### 3.4. Estimation of Parameters from Measured BTCs

The inverse models of ADE, MIM (from the nonlinear least-squares optimisation software CXTFIT2.0 [39]), and CTRW (from the CTRW Matlab Toolbox v.4.0 [44]), subject to an initial condition, Robin boundary condition at the inlet, and Neumann boundary condition at the outlet,
(16)Cx,⋅t=0=0,(17)vxC−D∂C∂xx=0=vxC0,(18)∂C∂xx=L=0,were fitted to the measured BTCs to obtain the transport parameters: vADE, DADE for the ADE model; vMIM, DMIM, βMIM, and ωMIM for the MIM model; and vCTRW, DCTRW, and βCTRW for the CTRW framework. The values of vADE and DADE were used as the initial guesses for the CTRW inverse modelling. The choice of the Robin boundary condition at the upstream was based on mass conservation [46]. Besides, the infinite outlet condition applied to the finite system of length L was based on the assumption that the solute concentrations at the upstream were not altered by the downstream boundary [46]. The measured effluent concentrations (C) were normalised by the input concentration (C0).To assess the goodness of the fitted models, the coefficient of determination (R2) and global error (Ej) were computed as
(19)R2=1−∑i=1NCio−Cif2∑i=1NCio−Co¯2,(20)Ej=1N∑i=1NCio−Cif2,wherej is ADE, MIM, or CTRW model; Cio and Cif denote the observed and fitted tracer concentrations, respectively; Co¯ is the mean of the observed concentrations; and N represents the total number of observed concentrations obtained in a transport experiment.


## 3.1. ADE Model

The ADE model is based on two assumptions. First, the ADE model treats the transition of solute particles as a Markovian process, meaning the present trait of the solute not dependent on its history [37]. Second, the ADE model adopts a Fickian-driven dispersive flux of a solute, with the coefficient of dispersion assumed to be constant in spatial and temporal scales [38]. For one-dimensional nonreactive solute transport through a homogeneous, isotropic, steady flow porous medium, ADE is given by [28]
(4)∂c∂t=DADE∂2c∂x2−v¯∂c∂x,where C is the solute concentration (ML-3), v¯ is the average pore-water velocity (LT-1), and Dx is the coefficient of hydrodynamic dispersion (L2T-1), t is time (T), and x is a distance (L). The analytical solution of ADE for flux type upstream boundary condition is given as [39]
(5)cx,t=c02erfcx−v¯t2DADEt+c0v¯2tπDADEexp−x−v¯t24DADEt−c021+v¯xDADE+v¯2tDADEexpv¯xDADEerfcx+v¯t2DADEt,where erfc means the complementary error function.


## 3.2. MIM Model

For the MIM model, which describes the physical nonequilibrium part applied in this study, the transport of a solute plume is modelled on the assumption that the porous media through which the plume migrates and has two different liquid phases: mobile (flowing) and immobile (stagnant). The migration of a solute in the flowing region is driven by an advection-dispersion process, while the immobile region exchanges its solute content with the mobile region through a rate-limited diffusion process and is described as a first-order process [25, 39]. For one-dimensional transport of a solute in the absence of chemical reaction, MIM is given as
(6)θm∂cm∂t+θim∂cim∂t=θmDm∂2cm∂x2−θmvm∂cm∂x,(7)ωcm−cim=θim∂cim∂t,where θm, θim, cm, and cim, respectively, denote the amounts of water and concentration in the flowing and stagnant liquid regions; vm and Dm, respectively, define the mean flow velocity and coefficient of hydrodynamic dispersion in the flowing fluid region; and ω is the first-order coefficient of mass transfer between the two flowing and stagnant liquid regions (T-1). The flowing liquid fraction is expressed by a dimensionless β=θm/θm+θim.


## 3.3. CTRW Model

The idea of CTRW was introduced in electronics to explore transport in semiconductors [40] and later applied in geological media to study solute transport [41]. In porous media, a solute is transported by a fluid in which it dissolves through tortuous paths. This transport process is conceptualised by CTRW as particles randomly undergoing a series of “hops” and each “hop” is described by a displacement and waiting time. Here, the displacement of the particles is defined by the physical and chemical heterogeneous nature of the media, and the waiting times describe the durations of hopping events. The displacements and the waiting times are independent random variables, and their distributions are governed by a joint probability density function (PDF), represented by ψs,t. For simplicity, the PDF is decoupled as ψs,t=Psψt, where Ps defines the probability distribution of the distance travelled by the particles between sites and ψt is the probability rate for the transition time to cover the transition displacement [42].As a solute migrates in porous media, a portion of the solute is held in immobile regions and later diffuses slowly out of the region to produce late temporal tailed BTCs, a mark of non-Fickian behaviour. The capturing of the waiting times of the solute held in the immobile region is the core of the CTRW model. This means that CTRW places more emphasis on the transition time distribution than the transition length since the residence times of fluids and solutes quantify transport in porous media. In this study, the concept of CTRW is summarised and the detail of the importance of waiting time distribution, which the CTRW model is all about, is documented in the publications [7, 43].Assuming a Gaussian distribution for the transition displacement, the one-dimensional solute concentration defined in Laplace space for a uniform system is expressed as [44]
(8)uc~x,u−c0x=−M~uvψ∂∂xc~x,u−Dψ∂2∂x2c~x,u,where
(9)M~u≡t¯uψ~u1−ψ~u.In the above equations,M~u is a memory function and is governed by the waiting time distribution [7]. M~u accounts for undetectable heterogeneities that exist at all scales, and u denotes the Laplace variable. c~x,u is Laplace transformed ensemble-average, normalised concentration, whereas c0 and x denote the initial concentration and location in space, respectively. The characteristic time, transport velocity, and generalised coefficient of dispersion are defined by t¯, vψ, and Dψ, respectively. The transport velocity, vψ, is not equivalent to the average pore velocity, v, as both are reported to be similar in the classical ADE. Likewise, Dψ is different from D in ADE.Considering a unit constant flux boundary conditionj~=u−1 at the inlet (x=0) of a finite system of length 1 (i.e., x=0,1), a zero concentration gradient boundary condition (i.e., ∂c~/∂x=0) at the outlet (x=1) and initial condition c0x=0, the resident concentration in the Laplace transform is given as
(10)c~u,x≡1uw−1eq+w+1erw+1M~vψ+2u∝ψew/∝ψ+w−1M~vψ−2u∝ψ,where
(11)w=1+4u∝ψM~vψ,(12)q=1+wx2∝ψ,(13)r=2w−xw2∝ψ.The behaviour of solute transport is quantified by the distribution of transition time, which can take on three different expressions, namely, asymptotic, truncated power law, and modified exponential models. In this study, only the truncated power model is discussed. The comprehensive description of the remaining models is documented in a publication by Cortis and Berkowitz [45]. To quantify the temporal regimes of transport in disordered media, characterised by heterogeneity at all scales, the transition time distribution of solute particles is defined by a truncated power law (TPL) [42]:
(14)ψt=t1τ2−βexpτ2−1Г−β,τ2−1−1exp−t/t21+t/t11+β,which transforms to Fickian behaviour at extended times. In Equation (14), τ2=t2/t1 and Г denotes the incomplete Gamma function. t1 denotes the characteristic median transition time and sets the lower limit for the power law transport behaviour such that t>t1. t2 is a cutoff time for the truncated power-law transition time distribution function that corresponds to the largest heterogeneity length scale. For the time regime t1≪t≪t2, ψt behaves in line with power law whereas ψt decreases exponentially for t≫t2 [42]. To find the memory function, TPL in Laplace transform (Equation (15)) [42] is substituted into Equation (9), with t¯=t1.
(15)ψ~u=1+τ2ut1βexpt1uГ−β,τ2−1+t1uГ−β,τ2−1.For the truncated power law transition time distribution, the nature of solute transport (i.e., non-Fickian or Fickian) is dictated by the parameterβ, where β is a constant exponent. The value of β is in three ranges, 0<β<1, 1<β<2, and β>2. The nature of transport for the interval 0<β<1 is extremely non-Fickian with early and long late-temporal tailed BTCs. For 1<β<2, the mean of the plume migrates with steady velocity and yields non-Gaussian-shaped BTCs with long late-time tailings. When β>2, the movement of the dissolved chemical species is already Fickian for t≫t1, and t2 becomes insignificant in the solute migration process. In this case, the solute plume mass centre migrates with the average fluid velocity and yields a dispersion coefficient that is constant spatially and temporally [42].


## 3.4. Estimation of Parameters from Measured BTCs

The inverse models of ADE, MIM (from the nonlinear least-squares optimisation software CXTFIT2.0 [39]), and CTRW (from the CTRW Matlab Toolbox v.4.0 [44]), subject to an initial condition, Robin boundary condition at the inlet, and Neumann boundary condition at the outlet,
(16)Cx,⋅t=0=0,(17)vxC−D∂C∂xx=0=vxC0,(18)∂C∂xx=L=0,were fitted to the measured BTCs to obtain the transport parameters: vADE, DADE for the ADE model; vMIM, DMIM, βMIM, and ωMIM for the MIM model; and vCTRW, DCTRW, and βCTRW for the CTRW framework. The values of vADE and DADE were used as the initial guesses for the CTRW inverse modelling. The choice of the Robin boundary condition at the upstream was based on mass conservation [46]. Besides, the infinite outlet condition applied to the finite system of length L was based on the assumption that the solute concentrations at the upstream were not altered by the downstream boundary [46]. The measured effluent concentrations (C) were normalised by the input concentration (C0).To assess the goodness of the fitted models, the coefficient of determination (R2) and global error (Ej) were computed as
(19)R2=1−∑i=1NCio−Cif2∑i=1NCio−Co¯2,(20)Ej=1N∑i=1NCio−Cif2,wherej is ADE, MIM, or CTRW model; Cio and Cif denote the observed and fitted tracer concentrations, respectively; Co¯ is the mean of the observed concentrations; and N represents the total number of observed concentrations obtained in a transport experiment.


## 4. Results and Discussion

### 4.1. Effect of Clay on Breakthrough Curves

The measured BTCs were analysed to explore the effect of clay on the nature of arrival and elution portions of the measured BTCs. The BTCs measured for HPM1-HPM4 at 1.68 mL/min are shown in Figure3. The measured BTCs showed a non-Fickian transport behaviour (i.e., early breakthrough and long late temporal tails of the solute plumes). The distinct shaped BTCs, indicating a varying degree of non-Fickian characteristics, were due to microscale heterogeneity induced by the increasing fraction of clay. The poorly or highly interconnected pores created by the clay led to the evolution of high-velocity and low-velocity fluid regions [47].Figure 3
Measured BTCs in HPM1-HPM4 for 1.68 mL/min.In HPM1, the BTC was near Gaussian in shape with its elution limb exhibiting a weak late-temporal tail. The result was consistent with previous studies [22, 23, 29] that the transport behaviour of solutes in even uniformly packed porous media will not be Fickian as generally assumed for a classical advection-dispersion equation model. The near Gaussian-shaped BTC could be a result of the narrow distribution of solute advective velocities induced by the uniformly packed spherical-shaped solid grains. Moreover, the weak solute plume tail at late times implied the occurrence of a small diffusion-dominated fluid region. An analogous weak tail characteristic was reported in a recent study by Crevacore et al. [48] who employed a comparable shape of solid particles as the component of the porous medium. Compared to HPM1, the solute transport in HPM2 also displayed a near Gaussian-shaped BTC with a low degree of late temporal tail. However, the solute particles arrived at the measurement plane earlier than observed in HPM1, defining the spatial memory of the medium. This breakthrough time was enhanced by the presence of preferential pathways, characterised by high advective velocities, along which the solute particles migrated quickly with small residence time. The preferential pathways could result from the effect of bridging [29] induced by 12.5% of clay in HPM2.In contrast, the non-Fickian characteristics (i.e., early breakthrough time and heavy late temporal tail) were more pronounced in HPM3 than in HPM2. The solute particles in HPM2 and HPM3 arrived at the measurement point earlier than that in HPM1. 12.5% and 25% fractions of clay in HPM2 and HPM3, respectively, were not enough to fill all primary pores between the sand particles; meaning, some pores were less filled or not filled at all. The unfilled pore spaces formed connected channels with enhanced hydraulic conductivity. The solute particles migrated along the preferential paths with relatively high velocities, explaining the earlier arrival times of solute in HPM2 and HPM3 compared to that in HPM1. The closeness of arrival times of solute in HPM2 and HPM3 could be attributed to the similarity in their preferential pathways along which the solute particles migrated. An increase in the proportion of clay from 12.5% to 25% in HPM3 did not only yield the earliest solute breakthrough time but also yield a solute elution limb, characterised by a heavy late temporal tail, compared to HPM1 and HPM2. This strong solute plume tailing was due to the increased stagnant fluid region, causing solute mass exchange between the flowing fluid regions. At early times, the stagnant fluid region served as a sink for the solute. The retention of the solute in the stagnant fluid region and subsequent release of held-up solute particles into the flowing fluid region delayed the mixing of the solute [27]. It is worth mentioning that the proportion of clay had a significant impact on the pace at which the solute diffused into and out of the stagnant fluid region.The solute transport in HPM1-HPM3 revealed an interesting trend of non-Fickian characteristics enhancement with an increasing proportion of clay. However, this trend did not continue in HPM4, consisting of equal fractions of sand and clay. The breakthrough time and late-temporal tail characteristics of HPM4 were almost the same as those of HPM1 even though HPM4 contained both sand and clay. Compared to HPM2 and HPM3, 50% of clay in HPM4 was enough to fill the primary pores between the sand particles, yielding well-connected and uniformly distributed pores with less immobile fluid regions. The well-connected and uniformly distributed pores explained the closeness of non-Fickian characteristics in HPM1 and HPM4. Our results indicated that in homogeneous mixed media, comprising sand- and clay-sized particles, the degree of non-Fickian transport of the solute was dependent on the proportion of clay present.


### 4.2. Inverse Modelling of Measured BTCs

To quantify the effect of clay on the solute transport in the filled single fracture, the BTCs measured in HPM1-HPM4 were analysed and interpreted using the inverse models of ADE, MIM, and CTRW. The representative fits to the BTCs for 1.68 mL/min are shown in Figures4–7. The coefficient of determination (R2) and global error (Ej) and the optimised transport parameters are presented in Tables 4 and 5, respectively.Figure 4
Measured BTCs fitted with the ADE, MIM, and CTRW models for HPM1 at 1.68 mL/min.Figure 5
Measured BTCs fitted with the ADE, MIM, and CTRW models for HPM2 at 1.68 mL/min.Figure 6
Measured BTCs fitted with the ADE, MIM, and CTRW models for HPM3 at 1.68 mL/min.Figure 7
Measured BTCs fitted with the ADE, MIM, and CTRW models for HPM4 at 1.68 mL/min.Table 4
Coefficient of determination (R2) and global error (Ej) for assessing the goodness of the fitted ADE, MIM, and CTRW inverse models for 1.68 mL/min.
Porous mediaADEMIMCTRWR2EjR2EjR2EjHPM10.99670.02730.99670.02720.99810.0206HPM20.99370.03680.99710.02510.99730.0240HPM30.99250.03740.99830.01780.99870.0159HPM40.99100.04450.99130.04370.99830.0195Table 5
Estimated parameters for the ADE, MIM, and CTRW models for HPM1-HPM4 at 1.68 mL/min.
ModelsParametersHPM1HPM2HPM3HPM4ADEvADE (cm/min)0.2040.2190.2200.225DADE (cm2/min)0.04700.2040.8520.0718MIMvMIM(cm/min)0.1730.1890.2260.212DMIM (cm2/min)0.03920.1530.1770.0626βMIM0.8490.8150.6720.9400ωMIM (min-1)0.004190.005091.5800.0178CTRWvCTRW (cm/min)0.2080.3030.3650.248DCTRW (cm2/min)0.04790.1761.7520.0722βCTRW1.961.751.631.87t1 (min)8.75×10−43.86×10−35.22×10−31.51×10−3t2 (min)3.42×1098.63×1061.52×1019.53×108The model of ADE performed relatively poorly in fitting the non-Fickian characteristics of the measured BTCs. In contrast, the models of MIM and CTRW showed a better performance in describing the measured BTCs, especially early arrival and late temporal tails, than the ADE model. Compared to MIM, the fitted results of CTRW were better than MIM. While CTRW allows multiple rates of waiting times of solute particles, MIM assumes a single rate for all diffusion-driven solute mass exchanges between the mobile and immobile fluid regions, which explains MIM’s relatively poor performance as discussed by Gao et al. [49]. According to Gao et al. [49], an alternative reason for the CTRW’s best performance was its five fitting parameters, which enhanced the degree of freedom in fitting the measured BTCs compared to the two and four fitting parameters of ADE and MIM models, respectively. The largest and smallest values of coefficient of determination (R2) and global error (Ej), respectively, revealed the best performance of the CTRW model as shown in Table 4.The ADE model produced an increasing global error. Adding 12.5%, 25%, and 50% of clay in HPM2-HPM4 led to increasing global error by 34.8%, 37%, and 63%, respectively. This showed that the classical ADE model could not quantify the characteristics of non-Fickian behaviour as the proportion of clay increased. For HPM1, which had no clay, the errors of the ADE and MIM models were very close with an error difference of 0.0001. However, as clay was considered in an increasing proportion in HPM2-HPM4, the error difference also increased (i.e., 0.0117, 0.0196, and 0.0008, respectively). HPM2 and HPM3 recorded large values, indicating the inefficiency of ADE to capture the enhanced non-Fickian characteristics. The results were consistent with the results by Scheibe et al. [50], who employed ADE and diffusion-driven mobile-immobile mass exchange model to quantify solute transport and found reasonable ADE fits for single particle-sized medium and poorer ADE fits for mixed media. The relatively low errors displayed by the MIM inverse model explained the presence of clay-induced physical nonequilibrium process, which was not accounted for in the classical ADE model.The difference between the average flow velocities estimated from the flow experiments (Table2) and the optimised transport velocities of the ADE model (Table 5) widened as the proportion of clay increased, and this could be explained by increasing microscale heterogeneity and the variation in flow velocities induced by the fractions of clay. As the clay was considered in varying fractions, the tortuosity of the porous media was enhanced and led to variable local velocities, which could be the reason for the widening of the difference between the average flow velocities and transport velocities. Moreover, comparing the velocity values obtained from the three inverse models (i.e., vADE, vMIM, and vCTRW), the vCTRW values were higher. This contrast was expected because the average velocities of solute particles and water molecules in the porous media were not the same [7]. vCTRW described the average solute particle whereas vADE and vMIM expressed the average water velocity, hence the reason for the contrast.The dispersion coefficient values (which quantify the spatial extent of solute through spreading) of ADE, MIM, and CTRW, as shown in Table5, increased with a growing proportion of clay except for equal proportions of clay and sand. An increase in dispersion coefficient values with a growing proportion of clay except for equal proportions of clay and sand could be explained by the tortuosity of the media and a variation in velocity distribution. For HPM1, which contained only spherical sand-sized particles, there was a narrow distribution of local velocities due to less tortuous flow paths. Because the velocity distribution was narrow, the spatial extent of solute through spreading was also narrow, explaining its smaller dispersion coefficient value. 12.5% and 25% of clay in HPM2 and HPM3, respectively, aided in the evolution of more tortuous pathways, producing a wide range of local velocities. As the velocity distribution became broader, the spread of solute through the pores also widened, explaining the larger dispersion coefficient values of HPM2 and HPM3. However, 50% of clay was large to fill the pore spaces between the larger spherical sand-sized particles in HMP4, yielding a well-connected and even distribution of pores (with a narrow velocity distribution). The narrow velocity distribution resulted in less solute spreading, hence the smaller dispersion coefficient value of HMP4 compared to that of HPM2 and HPM3. In all, the dispersion coefficient values of the MIM model were less than those of the ADE model. The dispersion coefficient values of the MIM model only described the spreading of the solute particles in the flowing fluid region, whereas the ADE model merged both solute dispersion in the flowing fluid region and diffusive solute mass exchange between the flowing and stagnant fluid regions, hence the reason for the relatively low dispersion coefficient values of the MIM model. There was a wide range of stagnant fluid regions in the porous media. For HPM1-HPM4, the stagnant fluid regions were 15.1%, 18.5%, 32.8%, and 6%, respectively. An increasing spread of solute was reflected in the enhancement of the early breakthrough of BTCs. According to Dou et al. [51], the optimised dispersion coefficients might not be enough to describe non-Fickian transport since the process of solute spreading only quantifies the early arrival characteristic of BTCs and not the late temporal tail, which is mainly controlled by solute mixing. In other words, the spreading of the solute plume did not necessarily mean the true mixing of the solute plume.The additional parameters (i.e.,βMIM and ωMIM) in diffusion-driven mobile-immobile mass exchange models are not only for improving the degrees of freedom in fitting measured BTCs but also to explain the physical connection between microscale geometrical configuration and diffusion rates of solute mass exchange [52]. Therefore, to further understand the effect of clay fractions on the physical nonequilibrium solute transport that resulted in late temporal tailed BTCs, the optimised values of βMIM and ωMIM were analysed. As presented in Table 5, fitting the measured BTCs with the MIM model revealed a variation in the flowing fluid region βMIM. In HPM1, the value of βMIM was 0.849, which implied that 15.1% of the medium’s fluid phase was stagnant. The occurrence of such a stagnant fluid regions was consistent with a previous study by Scheibe et al. [50] who observed trapped solute particles in the intergranular pore spaces in single particle-sized media when the simulated spatial solute distribution was visualised. The addition of clay (i.e., 12.5% and 25%) in HPM2 and HPM3, respectively, further reduced the fraction of the flowing region (βMIM) and increased the fractions of immobile fluid region (1-βMIM) in these two mixed media. For HPM2 and HPM3, the stagnant fluid region increased to 18.5% and 32.8%, respectively. However, for HPM4, which had 50% clay, the fraction of mobile region was 0.940, which implied that 6% of the medium’s fluid phase was stagnant. The largest mobile fluid region of HPM4 resulted from well-connected and even distributed pores caused by 50% of clay which filled the primary pore spaces between the sand-sized particles. Compared to 1-βMIM in HPM1, this 6% implied a reduction. The occurrence of stagnant fluid region showed a deviation from the Fickian solute transport behaviour. In addition, the variation in the stagnant fluid region indicated the significance of clay fractions in the development of stagnant fluid regions in mixed media consisting of sand and clay.The presence of stagnant fluid regions, according to Padilla et al. [27], has a considerable effect on solute mixing in terms of contact areas between flowing and stagnant regions. Therefore, we analysed the interactions between stagnant and flowing fluid regions in terms of solute mass exchange expressed by first-order mass transfer coefficient (ωMIM). In HPM1, which had no clay, the ωMIM value was small, connoting a low diffusion-driven solute mass exchange between the flowing and stagnant fluid regions. The small rate of solute diffusion explained the weak tailed BTC. Because the synthetic glass beads were impermeable [53], the low solute mass exchange resulted from solute particles trapped in potential dead-end pore spaces and narrow pores caused by the packing process [54].Compared to HPM1, the fraction of clay in HPM2-HPM4 enhanced the trapping of solute particles in the stagnant regions and increased the diffusive solute mass transfer between the flowing and stagnant fluid regions by 21%, 3761%, and 325%, respectively. These percentage increase values implied a higher degree of physical nonequilibrium process (i.e., mass exchange between mobile and immobile fluid regions), which reflected in the heavy late temporal tails, especially in HPM3. The heavy late-temporal tailed BTCs, a typical non-Fickian characteristic, implied an incomplete solute mixing process. The increasing trend ofωMIM values with growing clay fraction in HPM2 and HPM3 created a notion that the larger the clay fraction in mixed media (consisting of sand- and clay-sized particles), the higher the diffusion-driven solute mass transfer between flowing and stagnant fluid regions. However, in HPM4, which had 50% of clay, the ωMIM value was smaller compared to that of HPM3. This could be explained by the relatively less pore-scale heterogeneity in HPM4. An increasing rate of solute mass exchange between the two regions (with clay fraction less than the sand fraction) and a sudden drop in the solute mass transfer rate when both fractions were equal did show not only the significance of clay occurrence but also the importance of its proportion on non-Fickian characteristics enhancement in mixed media.The study also analysed the non-Fickian behaviour of the solute in HPM1-HPM4, using the remaining optimised parameters (t1, t2, and βCTRW) of the CTRW inverse model. The cut-off time values (t2), which served as a limit of the transition from a power-law model to an exponential model, were large, implying that the solute transport during the entire period of the laboratory experiments was non-Fickian. Moreover, compared to the time scale of the experiments, the characteristic median transition time values (t1) were small, and therefore, the transition time distribution was controlled mainly by the exponent of the power law (βCTRW). As presented in Table 5, the βCTRW values were in the range 1<βCTRW<2, implying a less non-Fickian transport behaviour. In this case, the solute transport velocity was steady, yet the coefficient of dispersion varied with the time scale. Even though the non-Fickian transport was less, the degree of non-Fickian characteristics (i.e., early arrival and late temporal tails) varied in HPM1-HPM4 as shown in Table 5. The variation in the degree of non-Fickian behaviour was due to a wide range of advective velocities and slow diffusion of solute particles into and out of stagnant fluid regions [29, 55, 56] induced by the fractions of clay.This study has strengthened the understanding that the contribution of clay-induced physical heterogeneity is significant and cannot be neglected in solute transport through mixed media as also suggested by Liu et al. [57].


### 4.3. Effect of Flow Rate on Solute Transport Behaviour

Since the power law exponent (βCTRW) expresses the nature of solute transport in porous media, the estimated βCTRW values of the CTRW fits to the measured BTCs for 0.96, 1.44, and 1.68 mL/min were employed to assess the impact of flow rate on the solute transport behaviour. Figures 8(a)–8(f) show the measured BTCs fitted with the ADE, MIM, and CTRW models for HPM1 and HPM3. In HPM1, increasing the pump rate (0.96, 1.44, and 1.68 mL/min) yielded the power-law exponent (βCTRW) values of 1.97, 1.974, and 1.96, respectively. The same trend was observed in HPM3. However, the βCTRW values (1.74, 1.71, and 1.63) were smaller than those in HPM1. The closeness of βCTRW values, especially in HPM1, indicated minor heterogeneity created by the flow rates.Figure 8
Measured BTCs fitted with the ADE, MIM, and CTRW models for HPM1: 0.96 mL/min (a), 1.44 mL/min (b), and 1.68 mL/min (c), and HPM3: 0.96 mL/min (d), 1.44 mL/min (e), and 1.68 mL/min (f).
(a)(b)(c)(d)(e)(f)For 0.96 mL/min, theβCTRW value in HPM1 was close to the threshold defining the onset of Fickian transport. The reason could be the greater interaction between the advective dominated region and diffusion dominated region, leading to greater solute mixing. As the rate of flow increases, the interaction between the two regions lessens, which amplifies non-Fickian behaviour [58]. This mechanism could be the reason for the relatively more non-Fickian transport as the flow rate increased, which was consistent with previous studies [29, 56]. The small yet significant distinction in the solute transport behaviour in the porous media at 0.96, 1.44, and 1.68 mL/min indicated the relevance of flow rate in enhancing non-Fickian solute transport, especially in the presence of clay, which was well quantified by the CTRW model with truncated power law transition time distribution.


## 4.1. Effect of Clay on Breakthrough Curves

The measured BTCs were analysed to explore the effect of clay on the nature of arrival and elution portions of the measured BTCs. The BTCs measured for HPM1-HPM4 at 1.68 mL/min are shown in Figure3. The measured BTCs showed a non-Fickian transport behaviour (i.e., early breakthrough and long late temporal tails of the solute plumes). The distinct shaped BTCs, indicating a varying degree of non-Fickian characteristics, were due to microscale heterogeneity induced by the increasing fraction of clay. The poorly or highly interconnected pores created by the clay led to the evolution of high-velocity and low-velocity fluid regions [47].Figure 3
Measured BTCs in HPM1-HPM4 for 1.68 mL/min.In HPM1, the BTC was near Gaussian in shape with its elution limb exhibiting a weak late-temporal tail. The result was consistent with previous studies [22, 23, 29] that the transport behaviour of solutes in even uniformly packed porous media will not be Fickian as generally assumed for a classical advection-dispersion equation model. The near Gaussian-shaped BTC could be a result of the narrow distribution of solute advective velocities induced by the uniformly packed spherical-shaped solid grains. Moreover, the weak solute plume tail at late times implied the occurrence of a small diffusion-dominated fluid region. An analogous weak tail characteristic was reported in a recent study by Crevacore et al. [48] who employed a comparable shape of solid particles as the component of the porous medium. Compared to HPM1, the solute transport in HPM2 also displayed a near Gaussian-shaped BTC with a low degree of late temporal tail. However, the solute particles arrived at the measurement plane earlier than observed in HPM1, defining the spatial memory of the medium. This breakthrough time was enhanced by the presence of preferential pathways, characterised by high advective velocities, along which the solute particles migrated quickly with small residence time. The preferential pathways could result from the effect of bridging [29] induced by 12.5% of clay in HPM2.In contrast, the non-Fickian characteristics (i.e., early breakthrough time and heavy late temporal tail) were more pronounced in HPM3 than in HPM2. The solute particles in HPM2 and HPM3 arrived at the measurement point earlier than that in HPM1. 12.5% and 25% fractions of clay in HPM2 and HPM3, respectively, were not enough to fill all primary pores between the sand particles; meaning, some pores were less filled or not filled at all. The unfilled pore spaces formed connected channels with enhanced hydraulic conductivity. The solute particles migrated along the preferential paths with relatively high velocities, explaining the earlier arrival times of solute in HPM2 and HPM3 compared to that in HPM1. The closeness of arrival times of solute in HPM2 and HPM3 could be attributed to the similarity in their preferential pathways along which the solute particles migrated. An increase in the proportion of clay from 12.5% to 25% in HPM3 did not only yield the earliest solute breakthrough time but also yield a solute elution limb, characterised by a heavy late temporal tail, compared to HPM1 and HPM2. This strong solute plume tailing was due to the increased stagnant fluid region, causing solute mass exchange between the flowing fluid regions. At early times, the stagnant fluid region served as a sink for the solute. The retention of the solute in the stagnant fluid region and subsequent release of held-up solute particles into the flowing fluid region delayed the mixing of the solute [27]. It is worth mentioning that the proportion of clay had a significant impact on the pace at which the solute diffused into and out of the stagnant fluid region.The solute transport in HPM1-HPM3 revealed an interesting trend of non-Fickian characteristics enhancement with an increasing proportion of clay. However, this trend did not continue in HPM4, consisting of equal fractions of sand and clay. The breakthrough time and late-temporal tail characteristics of HPM4 were almost the same as those of HPM1 even though HPM4 contained both sand and clay. Compared to HPM2 and HPM3, 50% of clay in HPM4 was enough to fill the primary pores between the sand particles, yielding well-connected and uniformly distributed pores with less immobile fluid regions. The well-connected and uniformly distributed pores explained the closeness of non-Fickian characteristics in HPM1 and HPM4. Our results indicated that in homogeneous mixed media, comprising sand- and clay-sized particles, the degree of non-Fickian transport of the solute was dependent on the proportion of clay present.


## 4.2. Inverse Modelling of Measured BTCs

To quantify the effect of clay on the solute transport in the filled single fracture, the BTCs measured in HPM1-HPM4 were analysed and interpreted using the inverse models of ADE, MIM, and CTRW. The representative fits to the BTCs for 1.68 mL/min are shown in Figures4–7. The coefficient of determination (R2) and global error (Ej) and the optimised transport parameters are presented in Tables 4 and 5, respectively.Figure 4
Measured BTCs fitted with the ADE, MIM, and CTRW models for HPM1 at 1.68 mL/min.Figure 5
Measured BTCs fitted with the ADE, MIM, and CTRW models for HPM2 at 1.68 mL/min.Figure 6
Measured BTCs fitted with the ADE, MIM, and CTRW models for HPM3 at 1.68 mL/min.Figure 7
Measured BTCs fitted with the ADE, MIM, and CTRW models for HPM4 at 1.68 mL/min.Table 4
Coefficient of determination (R2) and global error (Ej) for assessing the goodness of the fitted ADE, MIM, and CTRW inverse models for 1.68 mL/min.
Porous mediaADEMIMCTRWR2EjR2EjR2EjHPM10.99670.02730.99670.02720.99810.0206HPM20.99370.03680.99710.02510.99730.0240HPM30.99250.03740.99830.01780.99870.0159HPM40.99100.04450.99130.04370.99830.0195Table 5
Estimated parameters for the ADE, MIM, and CTRW models for HPM1-HPM4 at 1.68 mL/min.
ModelsParametersHPM1HPM2HPM3HPM4ADEvADE (cm/min)0.2040.2190.2200.225DADE (cm2/min)0.04700.2040.8520.0718MIMvMIM(cm/min)0.1730.1890.2260.212DMIM (cm2/min)0.03920.1530.1770.0626βMIM0.8490.8150.6720.9400ωMIM (min-1)0.004190.005091.5800.0178CTRWvCTRW (cm/min)0.2080.3030.3650.248DCTRW (cm2/min)0.04790.1761.7520.0722βCTRW1.961.751.631.87t1 (min)8.75×10−43.86×10−35.22×10−31.51×10−3t2 (min)3.42×1098.63×1061.52×1019.53×108The model of ADE performed relatively poorly in fitting the non-Fickian characteristics of the measured BTCs. In contrast, the models of MIM and CTRW showed a better performance in describing the measured BTCs, especially early arrival and late temporal tails, than the ADE model. Compared to MIM, the fitted results of CTRW were better than MIM. While CTRW allows multiple rates of waiting times of solute particles, MIM assumes a single rate for all diffusion-driven solute mass exchanges between the mobile and immobile fluid regions, which explains MIM’s relatively poor performance as discussed by Gao et al. [49]. According to Gao et al. [49], an alternative reason for the CTRW’s best performance was its five fitting parameters, which enhanced the degree of freedom in fitting the measured BTCs compared to the two and four fitting parameters of ADE and MIM models, respectively. The largest and smallest values of coefficient of determination (R2) and global error (Ej), respectively, revealed the best performance of the CTRW model as shown in Table 4.The ADE model produced an increasing global error. Adding 12.5%, 25%, and 50% of clay in HPM2-HPM4 led to increasing global error by 34.8%, 37%, and 63%, respectively. This showed that the classical ADE model could not quantify the characteristics of non-Fickian behaviour as the proportion of clay increased. For HPM1, which had no clay, the errors of the ADE and MIM models were very close with an error difference of 0.0001. However, as clay was considered in an increasing proportion in HPM2-HPM4, the error difference also increased (i.e., 0.0117, 0.0196, and 0.0008, respectively). HPM2 and HPM3 recorded large values, indicating the inefficiency of ADE to capture the enhanced non-Fickian characteristics. The results were consistent with the results by Scheibe et al. [50], who employed ADE and diffusion-driven mobile-immobile mass exchange model to quantify solute transport and found reasonable ADE fits for single particle-sized medium and poorer ADE fits for mixed media. The relatively low errors displayed by the MIM inverse model explained the presence of clay-induced physical nonequilibrium process, which was not accounted for in the classical ADE model.The difference between the average flow velocities estimated from the flow experiments (Table2) and the optimised transport velocities of the ADE model (Table 5) widened as the proportion of clay increased, and this could be explained by increasing microscale heterogeneity and the variation in flow velocities induced by the fractions of clay. As the clay was considered in varying fractions, the tortuosity of the porous media was enhanced and led to variable local velocities, which could be the reason for the widening of the difference between the average flow velocities and transport velocities. Moreover, comparing the velocity values obtained from the three inverse models (i.e., vADE, vMIM, and vCTRW), the vCTRW values were higher. This contrast was expected because the average velocities of solute particles and water molecules in the porous media were not the same [7]. vCTRW described the average solute particle whereas vADE and vMIM expressed the average water velocity, hence the reason for the contrast.The dispersion coefficient values (which quantify the spatial extent of solute through spreading) of ADE, MIM, and CTRW, as shown in Table5, increased with a growing proportion of clay except for equal proportions of clay and sand. An increase in dispersion coefficient values with a growing proportion of clay except for equal proportions of clay and sand could be explained by the tortuosity of the media and a variation in velocity distribution. For HPM1, which contained only spherical sand-sized particles, there was a narrow distribution of local velocities due to less tortuous flow paths. Because the velocity distribution was narrow, the spatial extent of solute through spreading was also narrow, explaining its smaller dispersion coefficient value. 12.5% and 25% of clay in HPM2 and HPM3, respectively, aided in the evolution of more tortuous pathways, producing a wide range of local velocities. As the velocity distribution became broader, the spread of solute through the pores also widened, explaining the larger dispersion coefficient values of HPM2 and HPM3. However, 50% of clay was large to fill the pore spaces between the larger spherical sand-sized particles in HMP4, yielding a well-connected and even distribution of pores (with a narrow velocity distribution). The narrow velocity distribution resulted in less solute spreading, hence the smaller dispersion coefficient value of HMP4 compared to that of HPM2 and HPM3. In all, the dispersion coefficient values of the MIM model were less than those of the ADE model. The dispersion coefficient values of the MIM model only described the spreading of the solute particles in the flowing fluid region, whereas the ADE model merged both solute dispersion in the flowing fluid region and diffusive solute mass exchange between the flowing and stagnant fluid regions, hence the reason for the relatively low dispersion coefficient values of the MIM model. There was a wide range of stagnant fluid regions in the porous media. For HPM1-HPM4, the stagnant fluid regions were 15.1%, 18.5%, 32.8%, and 6%, respectively. An increasing spread of solute was reflected in the enhancement of the early breakthrough of BTCs. According to Dou et al. [51], the optimised dispersion coefficients might not be enough to describe non-Fickian transport since the process of solute spreading only quantifies the early arrival characteristic of BTCs and not the late temporal tail, which is mainly controlled by solute mixing. In other words, the spreading of the solute plume did not necessarily mean the true mixing of the solute plume.The additional parameters (i.e.,βMIM and ωMIM) in diffusion-driven mobile-immobile mass exchange models are not only for improving the degrees of freedom in fitting measured BTCs but also to explain the physical connection between microscale geometrical configuration and diffusion rates of solute mass exchange [52]. Therefore, to further understand the effect of clay fractions on the physical nonequilibrium solute transport that resulted in late temporal tailed BTCs, the optimised values of βMIM and ωMIM were analysed. As presented in Table 5, fitting the measured BTCs with the MIM model revealed a variation in the flowing fluid region βMIM. In HPM1, the value of βMIM was 0.849, which implied that 15.1% of the medium’s fluid phase was stagnant. The occurrence of such a stagnant fluid regions was consistent with a previous study by Scheibe et al. [50] who observed trapped solute particles in the intergranular pore spaces in single particle-sized media when the simulated spatial solute distribution was visualised. The addition of clay (i.e., 12.5% and 25%) in HPM2 and HPM3, respectively, further reduced the fraction of the flowing region (βMIM) and increased the fractions of immobile fluid region (1-βMIM) in these two mixed media. For HPM2 and HPM3, the stagnant fluid region increased to 18.5% and 32.8%, respectively. However, for HPM4, which had 50% clay, the fraction of mobile region was 0.940, which implied that 6% of the medium’s fluid phase was stagnant. The largest mobile fluid region of HPM4 resulted from well-connected and even distributed pores caused by 50% of clay which filled the primary pore spaces between the sand-sized particles. Compared to 1-βMIM in HPM1, this 6% implied a reduction. The occurrence of stagnant fluid region showed a deviation from the Fickian solute transport behaviour. In addition, the variation in the stagnant fluid region indicated the significance of clay fractions in the development of stagnant fluid regions in mixed media consisting of sand and clay.The presence of stagnant fluid regions, according to Padilla et al. [27], has a considerable effect on solute mixing in terms of contact areas between flowing and stagnant regions. Therefore, we analysed the interactions between stagnant and flowing fluid regions in terms of solute mass exchange expressed by first-order mass transfer coefficient (ωMIM). In HPM1, which had no clay, the ωMIM value was small, connoting a low diffusion-driven solute mass exchange between the flowing and stagnant fluid regions. The small rate of solute diffusion explained the weak tailed BTC. Because the synthetic glass beads were impermeable [53], the low solute mass exchange resulted from solute particles trapped in potential dead-end pore spaces and narrow pores caused by the packing process [54].Compared to HPM1, the fraction of clay in HPM2-HPM4 enhanced the trapping of solute particles in the stagnant regions and increased the diffusive solute mass transfer between the flowing and stagnant fluid regions by 21%, 3761%, and 325%, respectively. These percentage increase values implied a higher degree of physical nonequilibrium process (i.e., mass exchange between mobile and immobile fluid regions), which reflected in the heavy late temporal tails, especially in HPM3. The heavy late-temporal tailed BTCs, a typical non-Fickian characteristic, implied an incomplete solute mixing process. The increasing trend ofωMIM values with growing clay fraction in HPM2 and HPM3 created a notion that the larger the clay fraction in mixed media (consisting of sand- and clay-sized particles), the higher the diffusion-driven solute mass transfer between flowing and stagnant fluid regions. However, in HPM4, which had 50% of clay, the ωMIM value was smaller compared to that of HPM3. This could be explained by the relatively less pore-scale heterogeneity in HPM4. An increasing rate of solute mass exchange between the two regions (with clay fraction less than the sand fraction) and a sudden drop in the solute mass transfer rate when both fractions were equal did show not only the significance of clay occurrence but also the importance of its proportion on non-Fickian characteristics enhancement in mixed media.The study also analysed the non-Fickian behaviour of the solute in HPM1-HPM4, using the remaining optimised parameters (t1, t2, and βCTRW) of the CTRW inverse model. The cut-off time values (t2), which served as a limit of the transition from a power-law model to an exponential model, were large, implying that the solute transport during the entire period of the laboratory experiments was non-Fickian. Moreover, compared to the time scale of the experiments, the characteristic median transition time values (t1) were small, and therefore, the transition time distribution was controlled mainly by the exponent of the power law (βCTRW). As presented in Table 5, the βCTRW values were in the range 1<βCTRW<2, implying a less non-Fickian transport behaviour. In this case, the solute transport velocity was steady, yet the coefficient of dispersion varied with the time scale. Even though the non-Fickian transport was less, the degree of non-Fickian characteristics (i.e., early arrival and late temporal tails) varied in HPM1-HPM4 as shown in Table 5. The variation in the degree of non-Fickian behaviour was due to a wide range of advective velocities and slow diffusion of solute particles into and out of stagnant fluid regions [29, 55, 56] induced by the fractions of clay.This study has strengthened the understanding that the contribution of clay-induced physical heterogeneity is significant and cannot be neglected in solute transport through mixed media as also suggested by Liu et al. [57].


## 4.3. Effect of Flow Rate on Solute Transport Behaviour

Since the power law exponent (βCTRW) expresses the nature of solute transport in porous media, the estimated βCTRW values of the CTRW fits to the measured BTCs for 0.96, 1.44, and 1.68 mL/min were employed to assess the impact of flow rate on the solute transport behaviour. Figures 8(a)–8(f) show the measured BTCs fitted with the ADE, MIM, and CTRW models for HPM1 and HPM3. In HPM1, increasing the pump rate (0.96, 1.44, and 1.68 mL/min) yielded the power-law exponent (βCTRW) values of 1.97, 1.974, and 1.96, respectively. The same trend was observed in HPM3. However, the βCTRW values (1.74, 1.71, and 1.63) were smaller than those in HPM1. The closeness of βCTRW values, especially in HPM1, indicated minor heterogeneity created by the flow rates.Figure 8
Measured BTCs fitted with the ADE, MIM, and CTRW models for HPM1: 0.96 mL/min (a), 1.44 mL/min (b), and 1.68 mL/min (c), and HPM3: 0.96 mL/min (d), 1.44 mL/min (e), and 1.68 mL/min (f).
(a)(b)(c)(d)(e)(f)For 0.96 mL/min, theβCTRW value in HPM1 was close to the threshold defining the onset of Fickian transport. The reason could be the greater interaction between the advective dominated region and diffusion dominated region, leading to greater solute mixing. As the rate of flow increases, the interaction between the two regions lessens, which amplifies non-Fickian behaviour [58]. This mechanism could be the reason for the relatively more non-Fickian transport as the flow rate increased, which was consistent with previous studies [29, 56]. The small yet significant distinction in the solute transport behaviour in the porous media at 0.96, 1.44, and 1.68 mL/min indicated the relevance of flow rate in enhancing non-Fickian solute transport, especially in the presence of clay, which was well quantified by the CTRW model with truncated power law transition time distribution.


## 5. Summary and Conclusions

To explore the influence of clay on solute transport in homogeneous mixed media as a single fracture infill, flow and nonreactive solute transport experiments in HPM1-HPM4 under three flow rates were conducted. The conventional ADE, MIM model, and CTRW framework with truncated power law transition time distribution were fitted to the measured BTCs to quantify the solute transport based on the optimised parameters. Our results allowed the following conclusions to be drawn.The measured BTCs showed non-Fickian behaviour (i.e., early arrival times and long late temporal tails) in HPM1-HPM4. However, the non-Fickian behaviour was stronger in HPM2-HPM4 (which contained some fractions of clay) than in HPM1 (which contained no clay). This implied that the presence of clay enhanced the non-Fickian transport.Fitting the classical ADE model to the measured BTCs yielded an increasing global error as the fraction of clay increased in HPM2-HPM4, affirming the inefficiency of classical ADE in capturing clay-induced non-Fickian transport. The models of MIM and CTRW performed better in describing the non-Fickian characteristics. However, CTRW’s performance was robust.Compared to HPM1, 12.5% and 25% of clay in HPM2 and HPM3, respectively, decreased the flowing fluid region and increased the solute transfer rate between the flowing and stagnant fluid regions in the MIM model. In addition, the CTRW model yielded power law exponent (βCTRW) values of 1.96, 1.75, and 1.63 in HPM1-HPM3, respectively, implying enhanced non-Fickian behaviour with increasing clay fraction. However, for HPM4, whose clay fraction was 50%, the βCTRW value was 1.87, indicating a deviation in the trend of non-Fickian enhancement with increasing clay fraction. This deviation showed that the enhancement of non-Fickian characteristics in homogeneous mixed media was dependent on the fraction of clay present.Moreover, increasing flow rate enhanced the non-Fickian behaviour based on the optimised values ofβCTRW.Our study has presented interesting results about the impact of clay and flow rate on nonreactive solute transport in homogeneous mixed media as an infill of a single fracture, which may be useful in groundwater contamination and remediation.



---

*Source: 1207971-2021-08-25.xml*