Technical Note: The use of an interrupted-flow centrifugation method to characterise preferential flow in low permeability media

We present an interrupted-flow centrifugation technique to characterise preferential flow in low permeability media. The method entails a minimum of three phases: centrifuge-induced flow, no flow and centrifuge-induced flow, which may be repeated several times in order to most effectively characterise multi-rate mass transfer behaviour. In addition, the method enables accurate simulation of relevant in situ total stress conditions during flow by selecting an appropriate centrifugal force. We demonstrate the utility of the technique for characterising the hydraulic properties of smectite-clay-dominated core samples. All core samples exhibited a non-Fickian tracer breakthrough (early tracer arrival), combined with a decrease in tracer concentration immediately after each period of interrupted flow. This is indicative of dual (or multi-)porosity behaviour, with solute migration predominately via advection during induced flow, and via molecular diffusion (between the preferential flow network(s) and the low hydraulic conductivity domain) during interrupted flow. Tracer breakthrough curves were simulated using a bespoke dual porosity model with excellent agreement between the data and model output (Nash–Sutcliffe model efficiency coefficient was > 0.97 for all samples). In combination, interrupted-flow centrifuge experiments and dual porosity transport modelling are shown to be a powerful method to characterise preferential flow in low permeability media.


Introduction
It is well known that heterogeneities, including biogenic pores/channels, desiccation cracks, fissures, fractures, nonuniform particle size distributions and inter-aggregate pores, are widespread in the subsurface and lead to a range of preferential flow phenomena (Beven and Germann, 1982;Cuthbert et al., 2013;Cuthbert and Tindimugaya, 2010;Flury et al., 1994). The coexistence of a relatively high hydraulic conductivity (K) domain(s) and an impermeable one, often termed dual porosity, results in a non-Fickian breakthrough curve. Solute transport in such systems is often characterised by an early arrival of solutes originating from the more mobile domain (macropores) and a slow approach to the final concentration caused by diffusion into the immobile domain (matrix or microporous network). When fitting breakthrough curves, therefore, it is often difficult to differentiate between contributions from the micro-and macropore transport mechanisms. As a consequence, in recent years there has been much research into the development of effective empirical and modelling techniques to characterise solute transport processes for dual porosity systems. One method investigated has been the use of interruptedflow solute-breakthrough experiments. Amongst the original work on this topic Murali and Aylmore (1980) discussed the influence of nonconstant flow on solute transport in aggregated soil. Brusseau et al. (1989) developed a flowinterruption method for use in measuring rate-controlled sorption processes in soil systems, which was subsequently applied by Koch and Fluhler (1993) to investigate advec-tion and diffusion phenomena occurring for nonreactive solute transport in aggregated media. The idea proposed was that, by interrupting flow during nonreactive tracer breakthrough, the degree of nonequilibrium between any fast-and slow-flow pathways can be determined. Central to this hypothesis is that the magnitude of the change in nonreactive tracer concentration in effluent samples taken immediately after a no-flow period is indicative of such nonequilibrium. Subsequent work within this field has included determination of physical (e.g. diffusive mass transfer between advective and nonadvective water) and chemical (e.g. nonlinear sorption) nonequilibrium processes in soil (Brusseau et al., 1997), determination of nonreactive solute exchange between the matrix porosity and preferential flow paths in fractured shale (Reedy et al., 1996), quantifying the effect of aggregate radius on diffusive timescales in dual porosity media (Cote et al., 1999), numerical modelling of aqueous contaminant release in nonequilibrium flow conditions (Wehrer and Totsche, 2003), empirical modelling of the release of dissolved organic species (Guimont et al., 2005;Ma and Selim, 1996;Totsche et al., 2006;Totsche, 2005, 2009) and heavy metals (Buczko et al., 2004), increasing the efficiency of solute leaching (Cote et al., 2000), empirical modelling of conservative tracer transport in a laminated sandstone core sample (Bashar and Tellam, 2006), and characterising in situ aquifer heterogeneity (Gong et al., 2010). One area where comparatively few studies exist, however, is in characterising the hydraulic properties of aquitards (e.g. clay-dominated soils and sediments, shales, and mudstones). Such research is of particular interest because preferential flow paths, by their intrinsic nature, can significantly compromise the integrity of aquitard units as local and regional barriers to the movement of groundwater contaminants. There are significant technical difficulties at present, however, in characterising such features at appropriate scales . For example, it is well known that the K of glacial till is scale dependent, with laboratory permeability measurements often yielding values lower than field-based measurements and modelling . As a consequence, a key requirement of laboratory-scale aquitard characterisation is that the core sample must be of sufficient volume in order to incorporate the key dual porosity features which govern the overall formation. A second technical challenge is that laboratory testing typically requires generation of flow through the sample whilst maintaining relevant in situ hydro-geotechnical conditions. One method which has been demonstrated as effective for this purpose is centrifugation, which is increasingly being used for hydraulic and geotechnical testing of low K materials (Hensley and Schofield, 1991;Nimmo and Mello, 1991;Timms et al., 2009;Timms and Hendry, 2008). Moreover, experiments using geotechnical centrifuges with payload capacities exceeding several kilograms can provide the additional benefit of being able to use core samples of representative scale for the overall formation. Here we present, for the first time, an interrupted-flow methodology using a centrifuge permeameter (CP) to characterise possible dual porosity behaviour of low permeability porous media. A novel dual domain model is also described which has been used to guide physical interpretation of the experimental tracer breakthrough curves.

Core and groundwater sampling methodology
The clay core (101.6 mm in diameter, Treifus core barrel, nonstandard C size) and groundwater were sourced from a 40 m thick, semi-consolidated, clay-rich alluvium deposit located approximately 100 km south of Gunnedah, New South Wales, Australia (31 • 31 9 S, 150 • 28 7 E). Equipment and procedures for obtaining minimally disturbed cores were compliant with ASTM (2012). See Timms et al. (2014) for a review of the procedure. Groundwater samples were taken from piezometers using standard groundwater quality sampling techniques (Sundaram et al., 2009). A 240 V electric submersible pump (GRUNDFOS MP1) and a surface flow cell were used to obtain representative samples after purging stagnant water to achieve constant field measurements of electrical conductivity (EC), pH, dissolved oxygen (DO) and reduction potential (Eh).

Centrifuge permeameter theory
During centrifugation, increased centrifugal force generates a body force which accelerates both solid and fluid phases within the sample. Centrifugal acceleration at any point within a centrifuge sample is calculated as follows: where a is the centrifugal acceleration (m s −2 ), ω is the angular velocity (rad s −1 ), and r is the radius from the axis of rotation (m). The g level is the scaling factor (a/g) for accelerated gravity, where g is gravity at Earth's surface. Vertical hydraulic conductivity, K v (m s −1 ), is calculated using ASTM (2000) (Eq. 2), where Q is the steady-state fluid flux (mL h −1 ), A is the sample flow area (cm 2 ), r m is the radial distance at the midpoint of the core sample (cm), and RPM is revolutions per minute.
The estimated in situ stress applied at the base of the core samples was calculated according to Eq. (3) and assumes that the overlaying formations were fully saturated and of a similar density to the core samples.
where σ i is the in situ stress (kPa), ρ s is saturated density of core (kg m −3 ), d is the depth to the base of the core sample (m below ground level (b.g.l.)); and g is the gravitational acceleration (m s −2 ). The applied stress at the base of the core (σ g , kPa) during the centrifuge experiments was calculated according to Eq. (4) (Timms et al., 2014).
where ρ b is the core bulk density (kg m −3 ), L c is the length of the CP core specimen (mm), ρ w is the influent density (kg m −3 ), h w is the height of influent water above the CP core specimen (mm), and a b is the centrifugal acceleration at the base of the CP core specimen.

Centrifuge permeameter sample preparation
A Broadbent geotechnical centrifuge (GMT GT 18/0.7 F) with a custom-built permeameter module (Timms et al., 2014) was used for this study. Prior to mounting into the CP, the outer 5 mm of the clay cores were trimmed and the trimmed cores were then inserted into Teflon cylindrical core holders (100 mm internal diameter, 220 mm length) using a custom-built mechanical cutting and loading device.
The cores were trimmed in order to remove any physical and chemical disturbance associated with the core extraction (drilling) process. A 5 mm thick A14 Geofabrics Bidim geofabric filter (100 micron, K = 33 m s −1 ) was placed above and below the sample in order to prevent clogging of the effluent drainage plate with colloid material from the sample. The geofabric filter was held in position above the sample using a plastic clamp. The core holders (with the core sample held within) were placed into 3000 mL glass beakers containing 1000 mL of groundwater derived from the piezometer at the closest depth to the core sample (see Table 1) and allowed to saturate from the base upwards. In total three core samples were analysed, which were taken from depths of 5.03, 9.52, and 21.75 m b.g.l. Saturation was performed by immersing the core holder into a reservoir of groundwater with the level of the water 5 cm higher than the top of the core sample. The mass of each core was then monitored every 24 h until no further increase in mass was recorded, saturation was then assumed to have occurred. The core holders (containing the saturated core samples) were mounted to the CP system via double O-ring seals. An influent head was added to all samples (see Table 1), which was maintained during centrifugation by a custom-built automated influent level monitoring and pumping system. The system comprises a carbon fibre EC electrode array which is connected via a fibre optic rotary joint to a peristaltic pump that supplies influent from an external 100 mL burette. Effluent samples were collected in an effluent reservoir and extracted using a 50 mL syringe. All experiments were conducted under steady-state flow, which is defined as a < 10 % difference between influent and effluent flow rates. The influent volume was determined by manual measurements of the water level in the external burette and effluent volumes were measured by multiplying their mass by their density.

Interrupted-flow experiment methodology
The idea of interrupting the flow during a breakthrough experiment is to differentiate between advection and diffusion processes. The method comprises a minimum of three phases.
1. Flow is induced at a constant centrifugal force for a fixed time period with effluent samples collected at multiple periodic intervals. The g level and influent reservoir height are selected so that the maximum total stress on the core approached the estimated in situ stress of the material at the given depth in the formation (Eqs. 3, 4). The time period between each effluent sampling interval is selected in order to gain sufficient effluent volume (namely > 1 mL) for accurate volume and nonreactive tracer concentration measurement.
2. Flow is interrupted (stopped) for a fixed time period during which time the permeameters are disconnected from the centrifuge module and positioned upright, the influent reservoir is also removed to limit any downward migration of solutes. A relatively long interrupted-flow period (> 12 h) is selected so that slow mass transfer processes can be identified.
All phases can be repeated multiple times in order to record sufficient nonreactive tracer breakthrough which enables the mass transport behaviour to be accurately characterised. Deuterium oxide (D 2 O) (Acros Organics, 99.8 % concentration) was used as a nonreactive tracer. A concentration of 3.12 mL L −1 was used, which raised the concentration of D 2 O to approximately 200 %. This was selected as sufficiently high in order to result in accurately measurable mass transfer changes. Effluent samples were filtered using a 0.2 µm cellulose acetate filter, stored at 4 • C and analysed for δD within 7 days of testing. δD was determined by measuring the 1 H/ 2 H ratio to an accuracy of 0.1 % using a Los Gatos DLT100 isotope analyser.

Dual domain transport modelling
Dual porosity models were created using COMSOL Multiphysics v. 4.4 (http://www.comsol.com) modified from wellknown formulations described, for example, by Coats and Smith (1964) and Bear and Bachmat (1990). The purpose of the modelling was to aid physical interpretation of the tracer breakthrough curves and validate the hypothesis that the step changes in tracer concentrations observed during no-flow periods could be explained by the presence of dual porosity in the samples. The models comprised a classical advectiondispersion equation for a mobile zone (subscript m) representing preferential flow pathways with a source/sink term representing exchange of solute with an immobile zone (subscript im). Solute transport in the immobile zone was by diffusion only. The exchanged flux between the immobile and mobile zones was modelled as being proportional to the concentration difference between the zones. The governing equations are as follows: where C is the δD isotope ratio (1), t is time (T), z is distance along the column (L), q is fluid flux (L T −1 ), α is hydrodynamic dispersivity (L), µ is the coefficient of molecular diffusion (L 2 T −1 ). The porosity, ∅, of the mobile and immobile domain is defined as where V p,m is the pore volume of the mobile domain (L), V p,im is the pore volume of the immobile domain (L) and V T is the total volume of the saturated core (L). The mass transfer coefficient, γ (T −1 ), is defined as where β is the dimensionless geometry coefficient, which typically ranges from 3 for rectangular slabs to 15 for spherical aggregates, and a is the characteristic half width of the matrix block (L) (Gerke and van Genuchten, 1993). The initial concentration conditions were set to zero for both domains for all model runs. During centrifugation periods, a variable solute flux upper boundary condition was used for the mobile domain and varied according to the product of the measured fluid flux and input concentration (C 0 ) during each experiment as follows: A Dirichlet (constant concentration) upper boundary condition was used for the immobile domain during times of centrifugation. A novel aspect of the models, facilitated by the flexibility of model structure variations possible in COM-SOL Multiphysics, was that the upstream transport boundary for both domains was switched to a zero flux condition during the interrupted-flow phases. The downstream transport boundary conditions for both domains were given by at z = L b , where L b was sufficiently large to ensure the results at the column outlet distance (at z = L c , L c L b ) were not sensitive to the position of the boundary. The total mass flux at the distance from the upstream boundary corresponding with the length of the experimental column was output from the models and integrated over the sampling periods for comparisons to the observed breakthrough curves. µ was calculated as 3.43×10 −5 m 2 d −1 which is the diffusion coefficient of D 2 O in H 2 O at 25.0 • C (Orr and Butler, 1935) multiplied by the average tortuosity of 0.15 reported by Barnes and Allison (1988) for clay bearing media. Model output was fitted to the observed data by varying the unconstrained parameters: α and γ . Note that ∅ m and ∅ im were also considered unconstrained parameters but their sum was constrained to equal total ∅ measured for each sample by oven drying at 105 • C for 24 h. In order to quantify the deviation between the recorded data and the dual porosity model, the normalised root mean square error (NRMSE) and the Nash-Sutcliffe model efficiency coefficient (NSMEC) were calculated (Nash and Sutcliffe, 1970). The mesh size and model tolerance were set sufficiently small so that the results were no longer sensitive to further reduction, to ensure the accuracy of the model output. The models runs presented were all executed using an extra fine mesh size and a relative tolerance of 0.00001.

Dual domain model sensitivity testing
Sensitivity analysis of the dual domain model (for the core taken from 5.03 m) was conducted in order to determine how sensitive the model was to changes in the constrained (L c , ∅ and µ) and unconstrained (∅ m , α and γ ) parameters. Sensitivity factors for constrained parameters were determined according to the estimated percentage error associated with each parameter, whilst ±50 % was selected for the unconstrained parameters in order to determine their influence on the NSMEC. The percentage error for L c was calculated to be ±2.78 % due to the core length being 36 mm, and the error associated with measurement at each end was ±0.5 mm. The percentage error for ∅ was calculated to be ±2.79 %, which comprises the L c measurement error plus 0.0026 % which is the calculated error associated with the two mass measurements. The percentage error for µ was determined to be ±50 % due to the range in tortuosity of 0.1-0.2 documented by Barnes and Allison (1988) and references therein.

D 2 O breakthrough
D 2 O breakthrough data and best-fit dual porosity model output for the interrupted-flow experiments conducted using core samples taken from 5.03, 9.52, and 21.75 m b.g.l. are displayed in Fig. 1. A close fit was achieved between the dual porosity model output and the original data, with a NS-MEC of 0.97, 0.99 and 0.97 and a NRMSE of 5, 3, and 5 % recorded for D 2 O breakthrough data from core samples taken from 5.03, 9.52, and 21.75 m b.g.l., respectively. The D 2 O breakthrough curves for all core samples exhibited a relatively elongated shape, with 100 % breakthrough not recorded for any of the timescales tested. This was expected given that a "long tailing" is a common feature of dual (or multi-)porosity materials, i.e. systems where the mobile domain is coupled to a less mobile, or immobile, domain. In such instances the dominant solute transport mechanism during imposed flow in the mobile domain(s) is typically advection; however, solute exchange also occurs in parallel with the immobile domain(s), typically via molecular diffusion. Following each interrupted-flow (no-flow) period a decrease in δD was recorded for all samples, and attributed to the diffusion of D 2 O from the preferential flow domain(s) into the low-flow, or immobile-flow, domain(s). The shape of the D 2 O breakthrough curves and the magnitude of the δD decrease following the interrupted-flow periods are different for all samples, with a 42.6, 18.5, and 28.4 % decrease recorded for the core samples taken from 5.03, 9.52, and 21.75 m b.g.l., respectively, after the first interrupted-flow period. In addition, the K v of each sample was recorded as different (Fig. 2), with average values of 1.4×10 −8 , 3.9×10 −9 , and 2.7 × 10 −9 m s −1 for the core samples taken from 5.03, 9.52, and 21.75 m b.g.l., respectively. The K v was recorded to decrease during the initial stages of each centrifugation period, attributed to the partial consolidation of the clay due to the stress applied by the centrifugal force. Following this initial consolidation period a more constant K v as a function of time was recorded for all cores, indicating that relative equilibrium had been achieved between stress applied by the centrifugal force and the compaction state of the core.

Dual domain model
The close model fits confirm that preferential flow through a dual porosity structure is a plausible hypothesis to explain the shape of the observed breakthrough curves. The unconstrained (∅ m , α and γ ) parameters that yielded the best dual domain model output fit to the D 2 O breakthrough data are displayed in Table 2. It is noted that the pore volume of the mobile domain per total volume of the core, ∅ m , was modelled to be 0.06, 0.04, and 0.08 for core taken from 5.03, 9.52, and 21.75 m b.g.l., respectively. With total porosity, ∅, measured as 0.44, 0.47, and 0.43, this equates to 13.6, 8.5, and 18.6 % of the total pore volume, respectively, suggesting that preferential flow features comprise a relatively large proportion of the total pore porosity in each sample. Hydrodynamic dispersivity, α, for best-fit model output for all core samples was L c /2, which is larger than typically reported for laboratory-scale column experiments (e.g. Shukla et al., 2003). It can be noted that all of the core samples were assumed to have remained saturated throughout the breakthrough experiments because all influent and effluent flow rates were recorded at steady state. Whilst dispersion is known to increase substantially as moisture content decreases from saturation (e.g. Wilson and Gelhar, 1981), it is therefore unlikely that this could have been a factor. The mass transfer coefficient, γ , was also modelled as different for each core sample with 0.65, 1.50, and 1.20 yielding the best model fit for the core samples taken from 5.03, 9.52, and 21.75 m b.g.l., respectively. Using Eq. (10), the half width of the matrix block (using a β range of 3-15 (3 for parallel slabs and 15 for spherical aggregates after Gerke and van Genuchten, 1993)), a, is calculated as within the range of 8.0-17.8, 5.4-12.1, and 5.5-12.3 mm for the core samples taken from 5.03, 9.52, and 21.75 m b.g.l., respectively. This suggests that the preferential flow channels present are likely to be separated by distances in the order of several millimetres from each other within the cores. With the dimensions of the cores significantly greater than these values, the model output therefore suggests that several preferential flow features are present in each core sample. Model output for the mobile and immobile domains at the top, middle, and base of the core samples is displayed in Fig. 3. It is noted that, for all core samples, diffusion into the immobile domain during the induced-flow periods is relatively significant, with δD im /δD m at the end of the first centrifugation (induced-flow) period recorded as 0.16, (right). The data points represent the concentration averaged over each sampling period and the dashed line for the model output represents the raw model output time series. In the empirical experiment it was therefore not possible to measure the concentration of the effluent during the no-flow phase because there was no effluent to collect for analysis. Thus, due to this averaging, in the rising limb of the breakthrough curve, the first point obtained by measurement during each flow phase can be observed as consistently greater than the "starting concentration" for the raw model output.  0.32, and 0.34 for the base of the core samples taken from 5. 03, 9.52, and 21.75 m b.g.l., respectively. With respective average flow rates recorded as 0.017, 0.007, and 0.015 m d −1 this behaviour is not obviously related to the variation in flow rates between the samples but more likely to the in-trinsic properties of the preferential flow domain (namely ∅ m , γ , and α). It is also noted that for all core samples full equilibration between the mobile and immobile domains occurred (δD im = δD m ) during each no-flow period. For example, δD im and δD m were modelled to be within ±1 % of each other after 7.0, 2.6, and 6.1 h during the first no-flow period for the core samples taken from 5.03, 9.52, and 21.75 m b.g.l., respectively.

Sensitivity analysis
Sensitivity analysis plots for a ±50 % change in unconstrained parameters (α, γ , and ∅ m ) for the core sample taken from 5.03 m b.g.l. are displayed in Fig. 4, with corresponding NSMEC data displayed in Table 3. The model fitting efficiency is relatively insensitive to all three unconstrained parameters in the range tested, with a less than 12 % change in the NSMEC compared to the NSMEC recorded for the best fit (Table 3). Sensitivity for the estimated percentage error associated with constrained parameters (∅, L c , and µ) are displayed in Fig. 5, with corresponding NSMEC data displayed in Table 3. The model fitting efficiency is also relatively insensitive, with a less than 1 % change in the NSMEC compared to the NSMEC recorded for the best fit (Table 3). For the data presented, the relatively low sensitivity to the parameters indicates that further testing, such as by dye tracing or geophysical tomography, is necessary to resolve more precisely the nature of the preferential flow paths. Nevertheless, the modelling has supported the preferential flow conceptual model we have used to explain the step changes in concentration observed after resting periods. It has also provided a first-order approximation of the likely geometry of the flow paths.

Comparison of dual and single domain modelling
In order to further demonstrate the practicality of the interrupted-flow methodology, a numerical experiment was carried out using the dual domain model developed above. Using the best-fit parameters from the core from 9.52 m b.g.l., an equivalent simulation to the laboratory experiment described above was run but without interrupted-flow phases. The breakthrough curve produced was then fit to the Ogata-Banks equation (Ogata and Banks, 1961) on the assumption that flow was occurring only through a single domain. The resulting fit was good (NRMSE = 3 %) with just one fitting parameter being the dispersion term which yielded a reasonable value of 1.27×10 −8 m 2 s −1 . This illustrates that, without the use of interrupted-flow phases to reveal the disequilibrium between two or more flow domains, a false assumption could easily be made with regard to the structure and associated transport properties of the core on the basis of a simple 1-D analytical model. This could have very significant consequences for the prediction and management of solute migration through such deposits. An additional numerical experiment was also undertaken to attempt to match the observed data to a single domain model which included resting phases, since no analytical solution is known for such a simulation. This was accomplished using COMSOL Multiphysics with identical settings to the dual domain models described above but with a disabled immobile domain. Calibrating to the δD breakthrough data recorded for the core from 9.52 m b.g.l. by just varying dispersivity, but using the measured porosity, we were unable to achieve a better fit than a NRMSE of 46 %, even with an unrealistically high dispersivity. A better fit is possible (NRMSE = 9 %, NSMEC = 0.9) if porosity is decreased to 0.1 but, again, only with an unrealistically high value for dispersivity of 1000L c (see Fig. 6). While such a model may be useful to suggest that the effective porosity of the core through which solute is moving is much less than the total porosity, it is only possible to fit the early time data (e.g. only the first flow stage) very accurately at the expense of the later time data. Perhaps more importantly than the lower NSMEC (or higher NRMSE) compared to the dual domain models, the single domain model also misses a key feature of the observed breakthrough curves: the decrease in concentration during resting phases. Instead, modelled concentrations   Table 3. NSMEC for the core sample taken from 5.03 m b.g.l. due to changes in constrained (L c , ∅, µ) and unconstrained (∅ m , α and γ ) model parameters. Changes in constrained parameters comprised the estimated percentage error per each parameter, which was 2.78, 2.79, and 50 % for L c , ∅, and µ, respectively. Changes in unconstrained parameters were ±50 %. The NSMEC for the best fit was 0.972.

Model
Pore increase during resting phases as would be expected in a single domain model due to redistribution of the solute along the core by diffusion. This additional numerical experiment thus strengthens the conclusions of the study, which are that dual domain behaviour is indicated by our interrupted-flow experiment observations, and that single domain models are inappropriate as a means of analysis.

Conclusions and outlook
Solute transport in the subsurface can be influenced by multiple nonlinear, rate-limited processes, and it is often difficult to determine which processes predominate for any given system. In this work we demonstrate the utility of interruptedflow solute transport experiments using a centrifuge permeameter to quantify the relative contributions of preferential flow pathways and surrounding matrix porosity to mass transfer processes in low permeability dual porosity materials. Dual domain transport modelling was used to validate the hypothesis that the step changes in tracer concentrations observed during no-flow periods could be explained by the presence of dual porosity in the samples. The modelling also enabled a first-order approximation of the physical properties of the two domains to be inferred. Smectite clay core samples were used (101.6 mm in diameter) as an example low K dual porosity media; however, it is anticipated that the methodology would also be suitable for the characterisation of any dual porosity material where mass transfer occurs via both advection and diffusion (e.g. fractured rock, heterogeneous soils, mine tailings). The methodology entails a minimum of three phases: induced flow, no flow, and induced flow; however, this may be repeated several times in order to most effectively characterise the multi-rate mass transfer behaviour. In addition, it is necessary to tailor the inducedflow rate, interrupted-flow timescales and nonreactive tracer concentrations in order to most effectively identify different mass transfer processes whilst also simulating realistic total stress conditions. Future work will seek to further investigate the structure of the clay samples studied using quantitative tomography techniques (e.g. X-ray computed tomography and magnetic resonance imaging) and how these physical features can be integrated into site-scale numerical flow modelling.