Effects of large-scale heterogeneity and temporally varying hydrologic processes on estimating immobile pore space: A mesoscale-laboratory experimental and numerical modeling investigation

Highlights

• Impermeable barriers and density dependence increase estimated immobile porosity.

• Decreasing dispersivity leads to increased estimated immobile porosity.

• Immobile porosity estimates were insensitive to changes in system flow rates.

Abstract

The advection-dispersion equation (ADE) often fails to predict solute transport, in part due to incomplete mixing in the subsurface, which the development of non-local models has attempted to deal with. One such model is dual-domain mass transfer (DDMT); one parameter that exists within this model type is called immobile porosity. Here, we explore the complexity of estimating immobile porosity under varying flow rates and density dependencies in a large-scale heterogeneous system. Immobile porosity is estimated experimentally and using numerical models in 3-D flow systems, and is defined by domains of comparatively low advective velocity instead of truly immobile regions at the pore scale. Tracer experiments were conducted in a mesoscale 3-D tank system with embedded large impermeable zones and the generated data were analyzed using a numerical model. The impermeable zones were used to explore how large-scale structure and heterogeneity affect parameter estimation of immobile porosity, assuming a dual-porosity model, and resultant characterization of the aquifer system. Spatially and temporally co-located fluid electrical conductivity (σ_f) and bulk apparent electrical conductivity (σ_b)—using geophysical methods—were measured to estimate immobile porosity, and numerical modeling (i.e., SEAWAT and R3t) was conducted to explore controls of the immobile zones on the experimentally observed flow and transport. Results showed that density-dependent flow increased the hysteresis between measured fluid and bulk electrical conductivity, resulting in larger interpreted immobile pore-space estimates. Increasing the dispersivity in the model simulations decreased the estimated immobile porosity; flow rate had no impact. Overall, the results of this study highlight the difficulty faced in determining immobile porosity values in field settings, where hydrogeologic processes may vary temporally. Our results also highlight that immobile porosity is an effective parameter in an upscaled model whose physical meaning is not necessarily clear and that may not align with intuitive interpretations of a porosity.

Introduction

Non-Fickian or “anomalous” solute transport, characterized by features such as early arrival times, elongated tailing, and concentration rebounding, are not adequately described by the classic advection-dispersion equation (ADE) (e.g., Sudicky et al., 1985; Benson et al., 2000; Kosakowski et al., 2001; Berkowitz et al., 2006; Culkin et al., 2008). Some studies attribute the advection-dispersion equation's (ADE) inability to capture such experimental observations to the presence of small-scale zones, or domains (i.e., relative to the measurement scale) that are “immobile” with respect to bulk flow where advection controls the transport (e.g., Ceriotti et al., 2019). These immobile domains are the result of the natural physical- and chemical-property variability of porous media across a range of scales and flow rates (Harvey et al., 1994; Haggerty and Gorelick, 1995; Liu and Kitanidis, 2012; Bouquain et al., 2012). Their importance can be explicitly quantified in terms of a local Peclet number, a dimensionless number that represents the ratio of advective and diffusive transport rates. The accessibility of solutes to immobile domains is thought to be velocity dependent and can be inferred from the analysis of the Damkohler (DaI) number, a dimensionless number that relates the physical mass transfer rate of solutes between the mobile and immobile domains to the advective transport rate. Mass transfer between the two domains can be assumed to be instantaneous for DaI >> 1 and the mobile and immobile domains treated as single domain. For values of DaI << 1, advection dominates and the contribution of the immobile pore space to total transport may be small, even though long tails may still exist (Bahr and Rubin, 1987). Quantifying the contributions and volume of the immobile pore space is important as they likely control the presence of secondary contaminant sources that render remediation processes inefficient and slow (e.g., Haggerty and Gorelick, 1994; Harvey et al., 1994; Wilson, 1997).

One of the simplest groups of models capable of representing immobile pathways are dual-domain mass transfer (DDMT) models, developed as a modification of the ADE. These types of models represent the non-equilibrium conditions between two domains (Coats and Smith, 1964; van Genuchten and Wierenga, 2010; Ma and Zheng, 2011), which cannot be captured through the use of a retardation coefficient, which assumes instantaneous equilibrium between mobile and immobile regions in the ADE. Lab-scale (mm to m) (e.g., Coats and Smith, 1964; Swanson et al., 2012, Swanson et al., 2015), field-scale (m to km) (e.g., Feehley et al., 2000; Zheng et al., 2011; Briggs et al., 2018) and numerical studies (e.g., Haggerty and Gorelick, 1995; Wheaton and Singha, 2010; Liu and Kitanidis, 2012; Baratelli et al., 2014) have shown that DDMT models can match observations better than the traditional ADE under a variety of conditions. DDMT models have even been successfully applied in systems where velocity differences between the mobile and immobile domains are small (e.g., Baratelli et al., 2011) and in fractured rock where diffusion was considered negligible (e.g., Becker and Shapiro, 2000).

DDMT conceptually treats one domain as immobile with respect to flow; the mass transfer between mobile and immobile domains is generally modeled as a first-order exchange process mimicking diffusive mass transfer, although more complex higher-order exchange processes can also readily be represented (e.g., Haggerty and Gorelick, 1995; Carrera et al., 1998). The resulting transport of a solute through the mobile and immobile domains can be described according to:

$${\theta }_m\frac{{\partial C}_m}{\partial t}=\nabla \bullet \left({\theta }_m{D}_m\nabla C_m\right)-\nabla \bullet \left({\theta }_{m}v\nabla C_m\right)-R$$

where θ_m is the porosity of the mobile domain [−], C_m is the dissolved species concentration in the mobile domain [kg/m³], v is the velocity in the mobile domain [m/s], and D_m is the hydrodynamic dispersion tensor [m²/s]. R is a source/sink term, where R_im is an exchange (reaction) term [kg/m³s] that is a component of R, which for a single-rate DDMT model is expressed as:

$$R_{\mathit{im}}={\theta }_{\mathit{im}}\alpha \left(C_m-C_{\mathit{im}}\right)$$

where θ_im is the porosity of the immobile domain [−], C_im is the dissolved species concentration in the immobile domain [kg/m³], and α is a first-order single-rate constant [1/s] and, along with other parameters within these transport models, is largely treated as a fitting parameter. In reality, concentration within “immobile” pore space may just have low advective velocities compared to the primary flowpaths. While some previous work has shown that these mass transfer parameters may be correlated to hydraulic conductivity or fracture length scales (e.g., Benson et al., 2001; Reeves et al., 2008), the physical meanings are generally considered to be difficult to correlate to anything meaningful in the field (e.g., Zhang et al., 2007; Willmann et al., 2008). Additionally, it has been recognized that velocity may control the estimated mass-transfer parameters (e.g., Haggerty and Gorelick, 1995) and they therefore display local Peclet number dependency (Ceriotti et al., 2019). While immobile zones are often conceptualized at the pore scale, field observations have indicated that they can be caused by large-scale heterogeneity as well (e.g., Briggs et al., 2018).

Recently, geophysical methods have been used to help constrain the mass-transfer parameters in Eq. (2)—specifically θ_im and α—during ionic tracer tests in the field given co-located fluid (σ_f) and bulk electrical conductivity (σ_b) measurements (e.g., Singha et al., 2007; Hampton et al., 2019; MahmoodPoor Dehkordy et al., 2018, MahmoodPoor Dehkordy et al., 2019). The σ_f measurements preferentially sample well-connected pore spaces (i.e., the mobile domain), whereas σ_b measurements are thought to collectively capture both the mobile and immobile domains (Singha et al., 2007); the σ_b measurements provide a way to “see” into the immobile domain and estimate the parameters from Eq. (2) at particular locations. In a single-domain medium, co-located σ_f and σ_b measurements are expected to respond similarly to an ionic tracer. In a dual-domain medium on the other hand, co-located σ_f and σ_b will deviate from one another due to differences in the solute loading and unloading rate(s) between the mobile and immobile domains; this difference allows these key mass-transfer parameters to be estimated.

While many column experiments have been used to calculate the parameters from Eq. (2) (e.g., Swanson et al., 2012, Swanson et al., 2015; Briggs et al., 2014), there are currently no larger-scale, controlled experiments reported in the literature that robustly explore the estimation of immobile porosity in the presence of large-scale heterogeneity. Here, we similarly look to explore estimates of θ_im in a mesoscale 3-D sand tank, with large impermeable blocks emplaced that produce a nonuniform flow field. These barriers produce low local velocities that result in tank-scale immobile zones. This work is motivated in part by that of Briggs et al. (2018), who hypothesized that their immobile porosity estimates were controlled by largely stagnant zones downgradient of large cobbles within glacially deposited lake sediments at their field site. We additionally explore the controls of density-dependent tracers, typically neglected in field tracer tests but clearly a process that affects transport in many settings (e.g., Beinhorn et al., 2005). Co-located measurements of σ_f and σ_b were collected at a variety of locations within the tank, and those data were then used with graphical and numerical models to estimate θ_im throughout the tank; inverted maps of σ_b were also used to image the transport of the tracer. A numerical model of the soil tank was also created as part of this work to support the experimental efforts. Flow velocities were varied in the model to control the global Peclet number to provide insight into how θ_im estimates and the contributions of advective and diffusive transport change in the presence of large-scale heterogeneities. As the generation of accurate data needed under controlled conditions is not feasible in the field, this study used an existing intermediate scale test tank where the boundary and initial conditions can be controlled and the soil conditions are accurately known.