AN ASSESSMENT OF SCALING AND INFORMATION REQUIREMENTS OF THE CLASSICAL ADVECTION DISPERSION EQUATION AND A
TEMPORAL-NONLOCAL ADVECTION DISPERSION EQUATION AT THE MADE SITE
by
c
Copyright by Savannah R. Miller, 2016 All Rights Reserved
A thesis submitted to the Faculty and the Board of Trustees of the Colorado School of Mines in partial fulfillment of the requirements for the degree of Master of Science (Hydrol-ogy). Golden, Colorado Date Signed: Savannah R. Miller Signed: Dr. David A. Benson Thesis Advisor Golden, Colorado Date Signed: Dr. Terri Hogue Professor and Department Head Department of Civil and Environmental Engineering
ABSTRACT
The MacroDispersion Experiment (MADE) Site in Columbus, MS, is a research site where tracer tests have been performed to further understand transport processes in heterogeneous material. Tracer tests have exhibited anomalous plume behavior, including early arrival times and heavy, power-law late-time tails that the classical Advection-Dispersion Equation (ADE) is unable to match. The heavy tailing is due to mass transfer into low velocity zones where dissolved solute becomes relatively immobile for certain periods of time. Transport modeling insufficiencies have led to two additional areas of research: collecting more hydraulic conductivity (K) data to better define subsurface heterogeneities and modifying the ADE to account for anomalous transport. This study seeks to further both research areas through the investigation of a single-well injection withdrawal test (SWIW) at the MADE Site with a skewed breakthrough curve (BTC). High resolution hydraulic conductivity and facies data were collected in the area of the SWIW test. This data was used to statistically generate a high resolution K domain of the aquifer material and was spatially upscaled (averaged) to examine transport model capabilities while destroying K information. Two models were compared at a range of K domain resolutions, the ADE and a temporal-nonlocal advection dispersion equation (t-fADE). The t-fADE model was used as an alternative to the classical ADE because it is able to simulate the power-law tailing observed in the SWIW BTC. Water table fluctuations throughout the SWIW test were also included by solving Richard’s equation.
When the K domain was defined at a high-resolution the ADE was able to match the heavy, power-law, late-time tailing. By explicitly defining small-scale heterogeneities, the ADE can simulate solute trapping in low K zones where solute becomes immobile for a period of time. As K information was destroyed, the ADE model’s accuracy decreased and simulating the SWIW test with the t-fADE model became more advantageous. Both trans-port model simulations exhibited their largest decline in performance at the scaling point
where facies boundaries were averaged out, indicating the importance of mass transfer be-tween facies rather than within facies boundaries. The vadose zone proved to be a significant contributor of solute trapping.
To explore the predictive capabilities of both transport models at such a highly hetero-geneous site, many equally probable realizations at the finest K resolution were created. The variability between realizations as well as simulated BTCs were analyzed to assess if the best-fit parameter values, which represent hydraulic properties of the subsurface, are characteristic of the aquifer or of the particular K field realization. While both transport models were able to fit the measured data with enough K information, the predictive capabil-ities were poor. Additional K field realizations could not repeatedly simulate the same BTC. Aquifer heterogeneity was too complex for a predictive transport model to be obtained.
TABLE OF CONTENTS
ABSTRACT . . . iii
LIST OF FIGURES . . . vii
LIST OF TABLES . . . x
CHAPTER 1 INTRODUCTION . . . 1
CHAPTER 2 BACKGROUND . . . 8
2.1 Groundwater Flow . . . 8
2.2 Solute Transport, ADE . . . 10
2.3 Anomalous Solute Transport . . . 10
2.3.1 Non-local Transport Models, t-fADE . . . 11
2.4 MADE Site . . . 14
2.5 SWIW Test . . . 15
2.5.1 Previous studies modeling the SWIW Test with the t-fADE . . . 18
CHAPTER 3 METHODOLOGY . . . 24
3.1 Hydraulic Conductivity Data . . . 24
3.2 Non-stationary Conditioned K Field with Hydrofacies . . . 24
3.3 Upscaling . . . 29
3.4 Fine K Realizations . . . 30
3.5 Model Simulations . . . 32
3.5.1 Parflow . . . 32
3.5.2 Slim-FAST . . . 33
3.6 Error Statistics and Model Comparison . . . 35
CHAPTER 4 RESULTS AND DISCUSSION . . . 36
4.1 Upscaling . . . 36
4.1.1 Upscaling of the ADE Model . . . 38
4.1.2 Upscaling of the t-fADE Model . . . 42
4.1.3 Comparison and Contrast of the ADE and t-fADE Models . . . 45
4.2 Analysis of Equally Probable K fields at the Finest Resolution . . . 52
CHAPTER 5 CONCLUSIONS . . . 56
CHAPTER 6 REFERENCES . . . 58
APPENDIX A - SUPPLEMENTARY FIGURES . . . 63
A.1 ADE versus t-fADE BTCs . . . 63
A.2 ADE and t-fADE Parameter Sensitivites . . . 65
LIST OF FIGURES
Figure 2.1 MADE site and location of the Single Well Injection Withdrawal (SWIW) site highlighted (Liu et al., 2010) . . . 16 Figure 2.2 Liu et al.’s (2010) simulated results for the ADE (AD) and SRMT (DDMT)
models. The best-fit model was the SRMT model with a zonal K domain. 18 Figure 2.3 Ronayne et al.’s (2010) final simulated BTCs for the a) ADE model and
b) SRMT (DDMT) model. The best simulated to measured match was the SRMT heterogeneous model. . . 19 Figure 2.4 Ibrahim’s (2013) log-log plot of the best-fit ADE and t-fADE BTCs. . . 20 Figure 2.5 Illustration of the three K model domains used in Stumb’s (2013) research:
a homogeneous domain (Case 1), layered domain (Case 2), and the highly detailed K domain (Case 3). KAand KGare the arithmetic and geometric
means respectfully. . . 22 Figure 2.6 Best-fit BTCs from the study conducted by Stumb (2013). Top is the ADE
simulated BTCs in the Layered Domain. Bottom is the t-fADE simulated BTCs in the Layered Domain. . . 23 Figure 3.1 MADE Site with the ICA highlighted and locations of GPR lines and DP
profiles from Dogan et al., 2011. . . 25 Figure 3.2 Four hydrofacies distinguished from the combined analysis of the GPR
data and DP measurements from Dogan et al., 2011. . . 26 Figure 3.3 Locations of the 9 DP measurements, domain boundary conditioning points,
and distinguished GPR facies within the 50 m by 50 m model domain. The extent of the facies represents a 10 by 10 meter area within the ICA cube. . 27 Figure 3.4 Illustration of creating a fBm K-field in 1-dimension. The first step is
creating a smooth trend of the original dataset. The second step is gen-erating an unconditional fBm by convolving Gaussian noise and the filter kernel. The third step is generating a smooth trend through the uncon-ditional fBm at the same locations as the actual K data. The fourth and final step is calculating the residuals of the unconditional fBm and the smoothed trend through the unconditional fBm and adding them to the original smoothed trend in Step 1. (Monnig, 2012). . . 28
Figure 3.5 osfBm hydraulic conductivity (K) field. This K field will be denoted f Bm K. The units of K are in meters per hour. . . 29 Figure 3.6 Slice through each upscaled hydraulic conductivity (K) domain from f bm
K to Homogeneous K (5). This illustration is to visualize the size of each grid block rather than the K value of each grid block so a color legend is not provided. . . 30 Figure 4.1 Time-lapse of particle injection, from left to right: a) finest (f Bm) K
field, b) middle (K2), and c) coarsest (Homogeneous) K field. Particle locations are recorded at hours 1, 4, 6, 10, 20, and 31. . . 37 Figure 4.2 ADE simulated BTCs for each hydraulic conductivity (K) domain
reso-lution (standard and log-log form): F Bm K to Homogeneous K. The measured BTC is also shown. . . 40 Figure 4.3 Particle locations for the fBm K field at the end of the extraction period.
Particles remaining in the model domain account for unrecovered mass. The well is located in the center of the figure and the black line represents the steady state water table. During the extraction period, the water table lowers to 8 m because of drawdown. . . 41 Figure 4.4 t-fADE simulated BTCs for each hydraulic conductivity (K) domain
res-olution (standard and log-log form): F Bm K to Homogeneous K. The measured BTC is also shown . . . 43 Figure 4.5 ADE and t-fADE simulated BTCs at the finest hydraulic conductivity
resolution (F bm K) in log-log and standard form. . . 46 Figure 4.6 ADE and t-fADE simulated BTCs at the coarsest hydraulic conductivity
resolution (Homogeneous K) in log-log and standard form. . . 47 Figure 4.7 Percent Mass Recovery versus hydraulic conductivity domain
discretiza-tion for the measured, ADE model, and t-fADE model. . . 49 Figure 4.8 Weighted Sum of Squared Residuals (WSSR) versus hydraulic conductivity
domain discretization for the ADE and t-fADE models. . . 50 Figure 4.9 Akaike information criterion (AICc) for the ADE and t-fADE models. . . . 50 Figure 4.10 Optimized parameter value versus discretization for both the ADE and
t-fADE models. a) Top left: longitudinal dispersivity, αT, b) Top right:
transverse dispersivity, αT, c) Bottom left: porosity, φ, d) Bottom right:
Figure 4.11 Top a) ADE simulated BTCs for all fBm hydraulic conductivity field re-alizations. Bottom b) t-fADE simulated BTCs for all fBm hydraulic con-ductivity field realizations. The dashed line represents the orignial fBm BTC. . . 54 Figure 4.12 Functional boxplots of the Top a) ADE and Bottom b) t-fADE simulated
BTCs for all fine-scale K realizations in standard and log form. The blue curves are the envelopes and the black line is the median curve. Dark ma-genta denotes the 25% central region, mama-genta represents the 50% central region, and the pink areas are the 75% central region. . . 55 Figure A.1 ADE and t-fADE simulated BTCs from the K0 hydraulic conductivity
domain in log-log and standard form. . . 63 Figure A.2 ADE and t-fADE simulated BTCs from the K1 hydraulic conductivity
domain in log-log and standard form. . . 63 Figure A.3 ADE and t-fADE simulated BTCs from the K2 hydraulic conductivity
domain in log-log and standard form. . . 64 Figure A.4 ADE and t-fADE simulated BTCs from the K3 hydraulic conductivity
domain in log-log and standard form. . . 64 Figure A.5 ADE and t-fADE simulated BTCs from the K4 hydraulic conductivity
domain in log-log and standard form. . . 65 Figure A.6 Composite Scaled Sensitivites of the optimized parameters in the ADE:
porosity, longitudinal dispersivity, and transverse dispersivity. . . 65 Figure A.7 Composite Scaled Sensitivities of the optimized parameters in the t-fADE:
porosity, longitudinal dispersivity, transverse dispersivity, and the capacity coefficient, beta. . . 66
LIST OF TABLES
Table 2.1 Upscaling process in Ibrahim’s (2013) study. . . 20 Table 3.2 The mean, variance, and standard deviation of ln(K) for each fractional
Brownian motion (f Bm) hydraulic conductivity (K) field realization. . . . 31 Table 3.3 Parameter values needed for Parflow simulations . . . 32 Table 4.1 Summary of Upscaling model runs, their Weighted Sum of Squared
Resid-uals (WSSR), and % Mass Recovery. . . 36 Table 4.2 Optimized parameter values for the ADE and t-fADE model for all
hy-draulic conductivity (K) domains. Parameters: Porosity (φ), Longitu-dinal Dispersivity (αL), Transverse Dispersivity (αT), and the Capacity
CHAPTER 1 INTRODUCTION
Accurately predicting solute movement through an aquifer is an important component of groundwater remediation processes. After a contaminant has entered the groundwater system, predicting the time at which peak concentrations will arrive at a particular loca-tion, or the time it takes for background concentrations to return, is crucial for developing cost-effective remediation strategies and estimating cleanup costs. The advection-dispersion equation (ADE) is commonly used to model solute movement through a porous medium. The ADE defines transport by advective and dispersive fluxes. Advective transport occurs along with the mean velocity of groundwater and dispersive transport accounts for diffusive and mechanical dispersion. Mechanical dispersion is defined by velocity distributions and variability of pore sizes at the pore scale as well as differences in hydraulic conductivity (K) at larger scales (Gelhar and Axness, 1981). Similar to many equations used to represent natural systems, the ADE makes several simplifying assumptions. One essential assump-tion of the ADE is that mechanical dispersion is analogous to Fick’s second law of diffusion (Anderson,1984). This assumption does not hold true for many field-scale systems with a large variance in groundwater velocity. In particular, the assumption appears to break down in highly heterogeneous geologic settings in which either the variance of velocity or corre-lation length of velocity (proportional to hydraulic conductivity (K) to the first order) is large. When the geology is highly variable and a field-scale model is needed, parameters are adjusted to account for areas where there are no data measurements. For groundwater trans-port models, the Fickian dispersion parameter is forced to increase and account for unknown K (Neuman, 1990). Solute dispersion, however, typically cannot be defined as Fickian until the plume has traveled several correlation lengths, i.e., a large distance from the source or a long period of time has elapsed (Anderson, 1984, Matheron and DeMarsily, 1980; Gelhar and Axness, 1981). When modeling solute migration over a short time period and/or a small
distance, Fick’s law can represent mechanical dispersion only when the heterogeneity and mixing serve as randomizing processes. If this is not the case, non-Fickian transport can oc-cur and result in non-Gaussian breakthrough oc-curves (BTCs). Numerous studies have shown that the ADE model is an insufficient method for predicting solute movement when non-Fickian dispersion is observed (Boggs et al, 1992; Harvey 1996; Eggleston and Rojstaczer, 1998; Benson et al., 2001; Schumer et al., 2003).
The idea that inadequacies of the ADE to model anomalous transport are due to unde-fined small-scale heterogeneities gives rise to the speculation that the ADE model may be sufficient if the K field is described at a fine enough scale. If this is true, at what scale might the ADE be valid? Hypothetically, if minuscule water fluctuations caused by aquifer hetero-geneity could be precisely defined, then a model could accurately predict solute movement based solely on a combination of advection and molecular diffusion (Zheng et al., 2011). Due to the insufficiencies of the ADE at the field-scale it has been suggested that if high-resolution K data can be obtained to populate a highly detailed velocity field, the predictive capabilities of the ADE in heterogeneous aquifers would improve. For any model domain, known parameters are averaged over the volume of a grid block. If there are multiple points of measurement within a grid block volume then that information is degraded by averaging the point measurements into one value to characterize a larger volume. A highly detailed domain, one with small grid block volumes at the same scale of measurement spacing, can directly account for small-scale heterogeneities (low permeability zones, preferential flow-paths, etc) without destroying information. To construct a detailed model domain, many grid blocks and data measurements are needed. It is important to be aware that this type of model domain will increase computation time and collecting this level of detailed K data is impractical for many contaminated sites.
Other researchers have proposed modifying existing governing equations to account for upscaled non-Fickian behavior through spatial and/or temporal nonlocality (Benson et al, 2001; Major et al., 2011). Spatial nonlocality defines concentration changes at a point as a
function of the concentrations at all upstream points in the aquifer (Schumer et al., 2001). Solute mass in temporal nonlocalities have memory of the times arriving and duration at a location (Schumer et al., 2003). These non-local equations contain additional parameters that replace detailed velocity information in, presumably, a more representative way. This study aims to further investigate both of these methods of improving solute transport mod-eling by: comparing the classical transport equation to a modified version that represents anomalous transport and obtaining high-resolution K data to create a highly detailed model domain.
Previous research using alternative transport models and collecting high-resolution K data have been conducted at laboratory and field scales. In a laboratory setting Klise et al. (2008) examined whether the ADE can be a valid model in an advection-dominated system if the K field is represented at a high enough resolution. Klise et al. (2008) exhaustively sampled a cross-bedded sandstone slab, 30.5 cm by 30.5 cm by 2.1 cm to obtain detailed permeability and porosity data at the sub-cm scale and coupled this with high-resolution imaging of solute transport. A numerical transport model was used to reproduce the solute tracer test experiment with the high-resolution data. They found that while the ADE could predict the bulk characteristics of solute transport, the intricate fingering observed was not adequately modeled (Klise et al., 2008). Experimental errors, such as boundary conditions, initial conditions, etc., were said to be minor and not the cause for poor model predictions from the ADE. Klise et al. (2008) believe that even though measurements were collected on the sub-cm scale, the major source of error was caused by the loss of material connectivity below the scale of measurement.
Major et al. (2011) continued this study by statistically downscaling the K data of the sandstone slab to incorporate with Klise et al.’s (2008) data. With this additional data they were able to create an even higher resolution K field. Major et al. (2011) then coarsened this K field (upscaled) by averaging grid block values to analyze scaling effects on transport parameters. They found that the ADE still did not fit the late-time tail of the BTC with very
fine-scale K data (Major et al., 2011). The discrepancy between simulated and measured BTCs is not well explained by subgrid K heterogeneity (Major et al., 2011). They also found that the longitudinal dispersivity parameter increased as the K field was coarsened. Major et al. (2011) simulated this same BTC with a temporal nonlocal equation, the time-fractional advection-dispersion equation (t-fADE), and found that it better predicted late-time concentrations even without detailed K information. The t-fADE model is based on the classical ADE but includes an additional fractional time derivative that accounts for mass transfer between mobile and immobile phases. At the small-scale, mass transfer occurs between connected and dead-end pores and at larger scales, mass transfer occurs between coarse material (sand and/or gravel) and clay layers or lenses (Benson et al, 2000). Other models have also included this type of mass transfer. Single-rate mass transfer (SRMT) and multi-rate mass transfer (MRMT) models also split K domains into mobile and immobile phases. The SRMT model, however, simulates an exponential tail rather than a power law tail. Coupling the use of high-resolution K data with nonlocal models has been performed at the field scale (Ronayne et al., 2010; Liu et al., 2007), but the benefits of using the t-fADE model in particular will be assessed.
The Macrodispersion Experiment (MADE) Site is an extensively studied field site with high levels of heterogeneity. The MADE site contains gravel, sand, and clay lenses all within close proximity to one another (Zheng et al., 2011). The original intent of the MADE site was to study large-scale transport theories. The existence of preferential flow paths and low permeability zones creates a highly complex aquifer that is difficult to characterize. Research endeavors to simulate BTCs at the MADE site with the ADE proved to be ineffec-tive in modeling plume migration and because of this further experimentation with revised methodologies was performed (Liu et al., 2010 and Ronayne et al., 2010).
A single-well injection withdrawal (SWIW) test was conducted by Liu et al. (2010) to analyze the predictive capabilities of the ADE as compared to a non-local single-rate dual-domain mass transfer (SRMT) model. A SWIW test or a push-pull test is an experiment
where a tracer is injected into a well and then extracted from the same well after a pe-riod of time. Liu et al. (2010) used a conservative bromide tracer in the SWIW test. A conservative tracer will not sorb or undergo any reactions during transport, which means impedances to transport are due to subsurface heterogeneity, such as low permeability zones and/or dead-end pores. The SRMT model divides the groundwater domain into mobile and immobile regions where mass transfer between the two is defined by an exponential function with a single mass transfer rate. Liu et al. (2010) found that the SRMT model always produced better results than the ADE, which is expected because it contains two additional parameters. The simulated BTCs, however, were unable to match the heavy-tail observed in the SWIW BTC. Both models resulted in less error when applied to a zonal K domain, rather than a homogeneous one emphasizing the importance of defining heterogeneity when simulating solute transport. Ronayne et al. (2010) performed a follow-up study on the SWIW test comparing the ADE to the SRMT model. They created a higher resolution K field that included connected channels generated using a multivariate Gaussian distribution. Ronayne et al. (2010) found that although the ADE model was improved with small-scale heterogeneity, it was still unable to reproduce the observed BTC, suggesting that even finer resolution K data may be needed for the ADE to be a sufficient model for this test (Ronayne et al., 2010). Both Liu et al. (2010) and Ronayne et al. (2010) modeled the SWIW test under saturated conditions and reported problems during the injection period. Water table mounding was observed, a process that can cause solute entrapment in the vadose zone. Dis-regarding solute entrapment in the vadose zone can lead to an inflated mass recovery. Since the SWIW test was simulated under saturated conditions they removed the upper portion of the screened interval during extraction to represent water table drawdown. Artificially re-moving a portion of the well causes a disconnect between particle pathways during injection and extraction moving the peak concentration to a later time. While this is a way to have better simulated fit, this is not a realistic way to replicate the SWIW test.
et al.’s (2010) work have led to further discussion on the proper methods and level of data collection needed to reproduce the SWIW test. This study will continue the efforts of Liu et al. (2010) and Ronayne et al. (2010) to simulate the SWIW test by comparing the ADE to a nonlocal model, but will replace the SRMT model with a more appropriate representation of the power law tail observed in the SWIW BTC. The t-fADE model can simulate a power law tail rather than an exponential tail making it superior to the SRMT with the same number of parameters. An additional improvement will be the incorporation of the vadose zone to account for water table fluctuations. These additions provide a more realistic approach and minimize the number of assumptions when simulating the SWIW BTC.
The need for higher resolution K data was addressed by Dogan et al. (2011) who collected high-resolution K data as well as GPR data that distinguishes different subsurface layers (facies) in the area of the SWIW test. This information was then used to generate multiple high-resolution K fields of equal probability. Multiple K fields were produced to analyze the variability in model simulations on a realization to realization basis. If the differences are large between equally probable K fields, then can one actually simulate contaminant transport in an aquifer with this degree of complexity? While we hope to improve ADE and t-fADE results by explicitly defining small-scale heterogeneities, this level of data collection is not likely in other field areas. To provide further information on the robustness of the two models the high-resolution K field will be coarsened (upscaled) multiple times to analyze the feasibility of accurately simulating contaminant transport in a heterogeneous aquifer with different levels of velocity information.
To reiterate, the research goals are to analyze:
1. if a high-resolution K domain improves the late-time ADE predictions.
2. how the ADE and t-fADE models compare while coarsening (upscaling) the K domain. 3. how simulated BTCs change when equally probable K fields are generated at the finest
While analyzing model simulations, an important component of relating this experiment to real world applications is considering the level of model and numerical complexity needed to make competent predictions.
CHAPTER 2 BACKGROUND
2.1 Groundwater Flow
The groundwater flow equation combines Darcy’s Law and the conservation of mass equation to describe the movement of groundwater through an aquifer. Darcy’s Law deter-mines the volumetric flux rate or specific discharge, q (LT−1), through porous media by the
equation,
q = −K∇h (2.1)
where K is the hydraulic conductivity tensor (LT−1) and h is the hydraulic head (L) (Fetter,
2001). For most groundwater flow the hydraulic head can be represented as h = pγ+z, where p is the pressure, γ is the unit weight of water, and z is the elevation head or potential energy associated with change in water height above some datum (Fetter, 2001). Although the specific discharge, q, has units of velocity, to calculate the average linear water velocity (v), the specific discharge is divided by the effective porosity (φ) to obtain (Fetter, 2001):
v = q
φ (2.2)
Darcy’s law is within the groundwater flow equation to solve for velocity. The conservative transient groundwater flow equation is:
∇ · (K∇h) = Ss
∂h
∂t (2.3)
Ss is the specific storage (L−1) and ∂h∂t is the change in hydraulic head per change in time
(LT−1) (Fetter, 2001). The Darcy flux and Groundwater Flow equations are applied to solve
for groundwater velocity within a discretized domain once the appropriate hydraulic head boundary and initial conditions are applied.
The groundwater flow equation is only valid for saturated conditions, meaning all pores are completely filled with water. During the SWIW test variably saturated conditions exist in certain portions of the test area because of water table fluctuations. Previously unsaturated soils can fill during injection and saturated soils may drain during extraction. These water
table fluctuations can result in possible solute entrapment in the vadose zone. In order to account for solute entrapment in variably saturated soils, models must include a portion of the vadose zone where pores are partially filled with water and air. The pressure head, which is the pressure of fluid within the pore space, in the vadose zone is less than atmospheric pressure and suction is the dominant force in fluid movement. Pressure head and hydraulic conductivity are both functions of water content (the ratio of fluid volume to total volume). Moisture retention curves are used to relate the pressure head, volumetric water content, and hydraulic conductivity for different types of soils. This relationship can then be plugged into Richard’s equation to solve for groundwater flow in variably saturated conditions. Richard’s equation is a nonlinear partial differential equation that combines the groundwater flow equation with the continuity equation (Delleur, 2007):
C(h)∂h ∂t = ∇ · K(h)∇( p γ + z) (2.4) C(h) = ∂θ ∂h (2.5)
where θ is the water content (volume of water, Vw, divided by the total volume, VT) (Fetter,
2001). The shape of the moisture retention curve is characterized by the van Genuchten model. Four independent parameters (α, n, θr, and θs) within the model are estimated from
soil-water retention data and used to relate water content (θ) and hydraulic head (h) by the following equation,
θ(ψ) = θr+
θs− θr
[1 + (αψ)n]1−n1
(2.6) where ψ is the suction pressure [L], θr is the residual water content [-], θs is the saturated
water content, α is related to the inverse of the air entry suction [L−1], and n is a measure of
the pore-size distribution (van Genuchten, 1980). No soil-water retention data were obtained for this study so van Genuchten parameter values for a similar soil were taken from Carsel and Parrish (1988). Richard’s equation could then be solved with these parameter values to determine hydraulic head values in unsaturated soils.
2.2 Solute Transport, ADE
The advection dispersion equation, based on the conservation of mass equation, is given by,
∂C
∂t = −∇ · (vC − ∇C) (2.7)
C(x, t) is the solute concentration [M L−3], v is the average pore water linear velocity [LT−1],
and D is the hydrodynamic dispersion tensor [L2T−1]. The ADE is a combination of
ad-vective flux (vC) and dispersive flux (D∇C). The adad-vective component of the ADE is the transport of solute by the bulk groundwater flow. The dispersion component combines me-chanical dispersion and molecular diffusion where Dij = (αL − αT)
vivj
|v| + (αT|v| + D ∗)I
ij
(Bear, 2004). Parameter α [L] is the dispersivity coefficient in the longitudinal (αL) and
transverse (αT) directions and D∗ [L2T−1] is the diffusion coefficient. Laboratory and field
experiments have found that dispersivity values are much greater in field-scale experiments than lab-scale experiments for the same material (Gelhar et al., 1992). The general reason-ing for this phenomenon is the influence of natural heterogeneities that cause irregular flow patterns at the field scale (Gelhar et al., 1992). Because heterogeneities are encountered with greater frequency as scale increases, longitudinal dispersivity is scale-dependent in the ADE (Hiscock, 2005). Diffusive processes are important when modeling at the pore scale but at larger scales they become negligible in comparison to mechanical dispersion.
2.3 Anomalous Solute Transport
The SWIW test produced non-Gaussian concentration-time BTCs with skewed early-time peaks and long heavy tails (Liu et al., 2010 and Ronayne et al., 2010). When simulating these BTCs, the ADE model underestimates measured concentrations at late-times. This tailing is attributed to the slow release of mass from immobile zones over time that causes low concentrations to persist (Zhang and Benson, 2008 and Major et al., 2011). While this could be due to undefined heterogeneity, heavy tailing has also been observed in relatively homogeneous material (Major et al., 2011; Benson et al., 2000; Schumer et al.,2001). In cases
where small-scale heterogeneity is not the cause for heavy tailing, solute entrapment in dead-end pores is likely the culprit. Because defining velocity at the pore-scale is extremely difficult or impossible, Darcy velocity often replaces pore-scale velocity as an upscaled average. By averaging velocity, information is lost and is transferred to dispersion parameters. In some cases increasing dispersion in the ADE is not enough to capture heavy tailing. Considering the poor predictive capability of the ADE to simulate the SWIW test and many others, alternative nonlocal models have been proposed in hopes of capturing anomalous plume behavior and asymmetric BTCs.
2.3.1 Non-local Transport Models, t-fADE
Nonlocal models define the concentration at a particular location based on concentrations at previous times and/or surrounding space (Benson et al, 2013). Non-locality arises due to loss of fine-scale velocity information when model domains are upscaled (coarsened) and grid blocks represent larger volumes. As models are upscaled and velocity information is lost, this information is transferred to the non-local operators (Benson et al. 2001; Morales-Casique et al., 2006). Temporal non-localities, which are dependent upon prior concentrations or loadings, are physical manifestations of the presence of mobile and immobile zones within the aquifer (Benson et al, 2001). Spatial non-locality arises when water from different upstream locations move into an upscaled measuring “point.” The classic case is a monitoring well (point) that intersects multiple stratigraphic layers of facies in the vertical direction and receives solute mass from many locations upstream (Benson et al., 2000). There are various types of nonlocal models, including single- and multi-rate mass transfer models, continuous time random walk (CTRW) methods, and fractional-order forms of the ADE.
This research will compare the ADE to the temporally nonlocal, time fractional-order form of the advection dispersion equation (t-fADE). The t-fADE model represents mobile zones, where solute movement is dominated by advection and mechanical dispersion, and immobile zones, where molecular diffusion processes dominate solute movement. Models
that define these two phases separately are called mobile/immobile (MIM) models and mass transfer between phases is described by one or more rate coefficients. MIM models parse ∂C
∂t
in the original conservation of mass equation (Equation 2.7) into a mobile (Cm(x, t)) phase
and an immobile (Cim(x, t)) phase (Schumer et al., 2001):
∂Cm
∂t + β ∂Cim
∂t = −∇ · (vCm− D∇Cm) (2.8)
Here, β = θim
θm is the mobile/immobile capacity coefficient and is equal to the ratio of immobile
porosity to mobile porosity (Coats and Smith, 1964). The total concentration is equal to Ctot = θmCm + θimCim. A simple constitutive relationship between Cm and Cim is the
single-rate mass transfer (SRMT) model, ∂Cim
∂t = ω(Cm− Cim), where ω is a first order rate
coefficient. Solving for ∂Cim
∂t as a function of Cm using Laplace transforms gives (Haggerty
and Gorelick, 1995): ∂Cim ∂t = ωe −ωt∗ ∂Cm ∂t + ωe −ωt(C m(x, 0) − Cim(x, 0)) (2.9)
For Equation 2.9, ∗ denotes a convolution, C(x, 0) are the initial concentrations in both phases, and f (t) = ωe−ωt is an exponential filter function or memory function. Equation 2.9
with the exponential memory function is the SRMT model and can be easily transitioned into a MRMT model. The memory function can take many forms, and assuming the initial concentration in the immobile phase equals zero gives a more general mass transfer equation of (Schumer et al., 2001),
∂Cim
∂t = f (t) ∗ ∂Cm
∂t + f (t)(Cm(x, 0)) (2.10)
The general mass transfer Equation 2.8 can then be plugged into the original conservation of mass Equation 2.7 to determine the mobile solute transport equation or multi-rate mass transfer equation,
∂Cm
∂t + β ∂Cm
∂t ∗ f (t) = −∇ · (vCm− D∇Cm) − βCm(x, 0)f (t) (2.11) Benson et al. (2004) state that the functional form of the BTC tail is a characteristic tool for determining underlying transport/retention processes. For example, the SRMT model predicts late-time concentrations at a fixed distance with an exponential tail, but often
power-law BTCs are observed (Major et al., 2011). Schumer et al. (2003) introduced the power-law memory function to replicate the observed tailing. Power-law heavy tails are a result of the immobilization of solute in heterogeneous media (Schumer et al., 2003).This memory function defines the time it takes for a particle to reenter the mobile zone from the immobile zone based upon a power law probability density distribution (Schumer et al, 2003). BTCs with power law tails have been observed at the MADE site and are observed from the SWIW test. This aspect of the BTC can be clearly seen by plotting concentration versus time on a log-log plot emphasizing the low concentrations that persist for long periods of time. As a result, the t-fADE better simulates late time mass decay than the exponentially decaying SRMT model. The exponential memory function is replaced with a power law function, f (t) = t−γ
Γ(1−γ), where Γ(·) is the Gamma function (Schumer et al., 2003; Dentz et
al., 2003), ∂Cm(x, t) ∂t ∗ f (t) = ∂γC m(x, t) ∂tγ (2.12)
Equation 2.12 is a Caputo fractional derivative of order γ (0<γ<1). Parameter γ (-) is the time-fractional exponent. The closer γ is to 1 the steeper the tailing slope and the faster the solute concentration falls off in the late-time tail of the BTC. Plugging Equation 2.12 into Equation 2.11 gives the t-fADE for mobile solute concentration,
∂Cm ∂t + β ∂γC m ∂tγ = −∇ · (vCm− D∇Cm) − βCm(x, 0) t−γ Γ(1 − γ) (2.13) Parameters β and γ are somewhat inversely related; larger β corresponds to higher immobile to mobile porosity causing more retention in solute transport while larger γ lessens retention by increasing the decay rate of solute concentration at late-times. Since the t-fADE defines immobile zones and slow moving mass, it is hypothesized that the values in the hydrodynamic dispersion tensor, D, might remain relatively constant at multiple scales (Major et al., 2011). How D changes with the scale of measurement is a key difference between the t-fADE and the ADE. We suspect the t-fADE model results will show similar values of D while coarsening the K field.
experiment to compare the classical ADE and the t-fADE models. Based on Klise et al.’s conclusions that discrepancy between simulated and measured BTCs was due to unaccounted for heterogeneity below the scale of measurement, Major et al. (2011) generated finer-scale K realizations that statistically honored the true K data. At the finest scale, the ADE could still not fit the late-time tail of the BTC. They determined that the mismatch between the modeled BTC and the observed BTC is not well explained by additional unsampled K heterogeneity within the sandstone slab (Major et al., 2011). The t-fADE model however was able to accurately simulate the late time tail of the BTC without sacrificing fit elsewhere. Major et al. (2011) concludes that the nonlocal model is necessary to simulate an advective-dominated system at all scales. They also challenged the previously conceived notion that the ADE could be used to model solute transport in highly heterogeneous media if finer detailed K data was obtained (Major et al., 2011). In this case, the level at which the K field was defined was very thorough and the likelihood that more information could be obtained is slight. There was no apparent small-scale representative elementary volume (REV) at which the transport is Fickian within the sandstone slab (Major et al., 2011). For the ADE model, the hydrodynamic dispersion coefficient, D, scaled linearly with the scale of K averaging whereas there were no apparent scaling effects of D for the t-fADE model (Major et al., 2011).
2.4 MADE Site
The Macrodispersion Experiment (MADE) site in Columbus, Mississippi has been an ongoing research field site for analyzing transport theories in heterogeneous media for almost 30 years (Zheng et al., 2011). The MADE site contains a shallow unconfined braided stream and fluvially deposited aquifer that has decimeter to centimeter scale preferential flow paths and clay lenses (Zheng et al., 2011). The aquifer is on average 10 m thick and is underlain by a continuous clay layer (Zheng et al., 2011). The overall variance of the natural logarithm of hydraulic conductivity (ln(K)) of 6.6, is reported to be much higher than previous tracer test
research sites (Bohling et al, 2012). Large- and small-scale tracer tests have been performed at the MADE site and high resolution K data has been collected throughout the area. ADE results from the tracer tests show large discrepancies between modeled and observed concentrations (Zheng et al., 2011; Bianchi et al., 2010). Preferential flow paths likely lead to early arrival time concentrations whereas slow release of contaminant from low permeability zones are responsible for the late time heavy-tailed concentrations seen in BTCs.
2.5 SWIW Test
An injection well was drilled in the Intensively Cored Area (ICA) at the MADE site for the SWIW test (Figure 2.1). The well was drilled to a depth of 9.75 m and has a screened interval of 3.66-9.75 m below the ground surface (Liu et al., 2010). The water table was measured by Dogan et al. (2011) to be 3.6 m below land surface and a confining aquitard was measured at a depth of approximately 12.2 m (Liu et al., 2010). The SWIW test conducted by Liu et al. (2010) consisted of five steps: (1) native water was injected into the well, (2) a known mass of bromide tracer was injected, (3) native water was injected once again to push the tracer further into the aquifer away from the well, (4) an inactive period of no injection or extraction, allowing tracer to enter low permeability zones by diffusion processes, and (5) water was pumped out of the well in an extraction phase, during which samples were simultaneously collected to measure bromide concentrations (Liu et al., 2010). The duration of the test was 410.3 hours. During the extraction phase the effluent concentrations were recorded to obtain a BTC and only 78.4% of mass was recovered (Liu et al, 2010).
Liu et al.’s (2010) analysis of the SWIW test showed improved fit when the SRMT model was used and also when heterogeneity was incorporated (Figure 2.2). These results are unsurprising because a model with two additional parameters and a more complicated K domain should improve the measured to simulated fit. Liu et al. (2010) created two K model domains of the SWIW test area. The first was a homogeneous K domain with a single K value of 5.27 (md−1) and the second was a zonal K domain with three facies: a high K
Figure 2.1: MADE site and location of the Single Well Injection Withdrawal (SWIW) site highlighted (Liu et al., 2010)
facies, low K facies, and an average K facies equal to the homogeneous case (Liu et al., 2010). The heterogeneous domain or zonal K domain used in Liu et al.’s study was fairly coarse; therefore, Ronayne et al. (2010) extended this study by collecting more detailed K data of the area. They created a sub-meter scale K field based on a correlated multivariate Gaussian distribution to input into the ADE and SRMT models. This K field represented subfacies heterogeneity i.e., heterogeneity within each facies or subsurface layer, along with highly connected channels for explicit simulation of preferential flows and slow advection (Ronayne
et al., 2010). Ronayne et al. (2010) found that when the small scale heterogeneities were accounted for, the ADE could better model the BTC (Figure 2.3). Parameters within the SRMT also decreased, revealing that small-scale heterogeneity between facies is responsible for local scale mass transfer. Although the ADE model improved, it still generated higher error than the SRMT model. Ronayne et al.’s (2010) overall conclusion was that extremely detailed data are needed to appropriately model the SWIW test contrary to Major et al.’s (2011) findings with the sandstone slab experiment. Both Ronayne et al. (2010) at the field scale and Klise et al.’s at the laboratory scale state the necessity of obtaining higher resolution K data for the ADE to produce accurate results. The gauge on how high this resolution must be is still uncertain however and is likely dependent on a site-by-site basis.
Liu et al. (2010) and Ronayne et al. (2010) both present normalized concentration versus time BTCs. By doing this they do not show differences in the true mass recovery and simulated mass recovery and it appears that their simulated peak concentrations are equivalent to the measured peak concentration when in actuality they are not. Distinguishing mass recovery is important because this accounts for mass that has been recovered and mass that remains in the aquifer in low permeability zones, dead-end pores, and/or the vadose zone. Liu et al. (2010) and Ronayne et al. (2010) also simulated the SWIW test under saturated conditions. Since the vadose zone was not included in model simulations they had to artificially remove an upper portion of the model domain during extraction to account for drawdown. This resulted in two problems; first, they were unable to capture mass trapped in the vadose zone and secondly, removing a portion of the model domain causes a disconnect in particle pathways and shifts the peak concentration to a later time. Normalized concentrations and simulating the SWIW test under saturated conditions are not a realistic ways to simulate the SWIW test, but rather a way to fit simulations to measured data.
Figure 2.2: Liu et al.’s (2010) simulated results for the ADE (AD) and SRMT (DDMT) models. The best-fit model was the SRMT model with a zonal K domain.
2.5.1 Previous studies modeling the SWIW Test with the t-fADE
Two additional studies were performed on the SWIW test after the initial results from Liu et al. (2010) and Ronayne et al (2010). Masters theses by Ibrahim (2013) and Stumb (2013) are the basis of the research presented in this study. The first, conducted by Ibrahim (2013), examined how upscaling influences a model’s capability to simulate the SWIW test under saturated conditions. Ibrahim compared two transport models while coarsening the K field, the classical advection dispersion equation (ADE) and the time-fractional advec-tion dispersion equaadvec-tion (t-fADE). He obtained high-resoluadvec-tion K data from Dogan et al. (2011) to create a model domain that defined small-scale heterogeneity in the location of the SWIW test. The high-resolution K domain was upscaled until a homogeneous domain was achieved (Table 2.1). He found that high-resolution hydraulic conductivity data did not im-prove results of the ADE model at late-times. The ADE simulated BTC consistently under
Figure 2.3: Ronayne et al.’s (2010) final simulated BTCs for the a) ADE model and b) SRMT (DDMT) model. The best simulated to measured match was the SRMT heterogeneous model.
predicted concentrations at late-times and the t-fADE BTC had better late-time predictions at all scales (Ibrahim, 2013) (Figure 2.4). Both models simulated peak concentrations earlier than the measured peak concentration (Ibrahim, 2013).
Table 2.1: Upscaling process in Ibrahim’s (2013) study. Grid Discretization (x, y, z) Grid block [m3]
128 x 128 x 640 0.00298 128 x 128 x 256 0.00597 128 x 128 x 128 0.01192 64 x 64 x 64 0.09541 32 x 32 x 32 0.76294 16 x 16 x 16 6.10352 8 x 8 x 8 48.8281 4 x 4 x 4 390.625 2 x 2 x 2 3125 1 x 1 x 1 25000 R e la ti v e C o n c e n tr a ti o n [ ] Withdrawal time [hr] Field data ADE: 128x128x512 t-fADE: 128x128x512 10-3 10-1 10-2 103 10-4 101 100 102
Figure 2.4: Ibrahim’s (2013) log-log plot of the best-fit ADE and t-fADE BTCs.
(2013) to create a fine scaled K domain, but added the vadose zone to capture water table fluctuations throughout the simulation (Stumb, 2013). She hypothesized that water table mounding and sequential drawdown trapped mass in the vadose zone and partly accounted for the low mass recovery measured by Liu et al. (2010) (Stumb, 2013). The addition of the vadose zone creates a more realistic model, but it also increases computation time be-cause it is necessary to solve the nonlinear Richards equation for each time step. Instead of upscaling the highly detailed K domain many times, Stumb only examined three K do-mains: the high-resolution K domain, a four-layered K domain that represents four distinct hydrofacies, and a homogeneous K domain (Figure 2.5) (Stumb, 2013). Stumb (2013) found that including the vadose zone improved the ADE results by lowering mass recovery closer to what was measured in the field (Figure 2.6). Although the ADE model was improved, the t-fADE model still produced better results for all K domains when the unsaturated zone was included (Stumb, 2013). The high-resolution K field improved model results, but she found that distinguishing facies boundaries was of more importance (Figure 2.6) (Stumb, 2013). Mass transfer between facies was of a higher order than between low K zones within facies boundaries.
Figure 2.5: Illustration of the three K model domains used in Stumb’s (2013) research: a homogeneous domain (Case 1), layered domain (Case 2), and the highly detailed K domain (Case 3). KA and KG are the arithmetic and geometric means respectfully.
0 50 100 150 200 250 300 350 400 450 0 10 20 30 40 50 60 70 Time [hours] Conc ent ration (ppm) 100 101 102 103 101 100 101 102 Time [hours]
Concentration (ppm) Saturated Geometric ADE Saturated Arithmetic ADE Variably Saturated Geometric ADE Variably Saturated Arithmetic ADE Field Data BTC 450 3 0 50 100 150 200 250 300 350 400 0 5 10 15 20 25 30 35 40 45 50 Time [hours] Concentration (ppm) 100 101 102 10 101 100 101 102 Time [hours] Concentration (ppm)
Saturated Geometric t fADE Saturated Arithmetic t fADE Variably Saturated Geometric t fADE Variably Saturated Arithmetic t fADE Field Data BTC
Figure 2.6: Best-fit BTCs from the study conducted by Stumb (2013). Top is the ADE simulated BTCs in the Layered Domain. Bottom is the t-fADE simulated BTCs in the Layered Domain.
CHAPTER 3 METHODOLOGY
3.1 Hydraulic Conductivity Data
The Intensively Cored Area (ICA) is a small portion of the MADE site where the SWIW test was conducted (Figure 3.1). Dogan et al. (2011) collected high-resolution 3D ground-penetrating radar (GPR) and combined this with high-resolution hydraulic conductivity (HRK) data from cm-scale Direct Push (DP) measurements to characterize the ICA. The DP measurements were collected using a DP HRK probe (Dogan et al., 2011). As the probe is pushed into the subsurface water is injected out of a small screen positioned near the tip of the probe (Dogan et al., 2011). The ratio of the injection rate to injection-induced back pressure can be transformed into hydraulic conductivity (K) following an approach by Liu et al. (2009). A total of 9 DP vertical profiles (Figure 3.3) and 3.8 km of GPR lines were collected in the ICA (Figure 3.1) (Dogan et al., 2011). The analysis of both datasets provided high-resolution estimates of the 3D hydrostratigraphic facies (Figure 3.2) A hydrostratigraphic facie is a sedimentary layer with composition and depositional environment distinct from the surrounding sediments (Reading, 1996). Four GPR facies were distinguished and show good correlation with the DP measurements (Dogan et al., 2011).
3.2 Non-stationary Conditioned K Field with Hydrofacies
Lognormal K fields were stochastically generated from the data collected by Dogan et al. (2011). The code to generate the K fields uses operator-scaling fractional Brownian motion (osfBm) to condition the recorded K points and interpolate a K domain (Appendix B). The osfBm code is different from typical fBm because scaling varies with direction (Monnig et al., 2008). This is important because material can have stronger correlations in certain directions. An example of this is material deposited by a braided stream. The material will have a stronger correlation in the direction of the braided stream.
Figure 3.1: MADE Site with the ICA highlighted and locations of GPR lines and DP profiles from Dogan et al., 2011.
An f Bm K field was created for each of the four facies and then spliced together to create a K field for the entire domain. Statistical parameters (mean and variance) are cal-culated for the DP measurements within each GPR defined hydrofacies and then used for K interpolation. The conditioning points for each hydrofacies are from the 9 DP measurements shown in Figure 3.3. The osfBm code can describe heterogeneity at multiple scales because the variance scales as a power law. To create a fBm K field, an initial trend surface, Kavg,
is obtained by ordinary kriging methods from the conditioning data, K. Convolving ran-dom Gaussian noise (B(~x)) to an appropriate filter kernel (ϕ(~x)) creates a three dimensional unconditional fBm field (Monnig et al, 2008):
Bϕ(~x) = B(~x) ∗ ϕ(~x) (3.1)
A Fast Fourier Transform (FFT) is used for the convolution and then the inverse FFT (iFFT) to obtain Bϕ(~x). The filter kernel, ϕ(~x), describes the scaling relationship and correlation
Figure 3.2: Four hydrofacies distinguished from the combined analysis of the GPR data and DP measurements from Dogan et al., 2011.
ϕ(x) = x−H+d2 (3.2)
where H is the Hurst coefficient(s) and d is dimension. The increments of Bϕ(~x) are
sta-tionary. These increments are described by fractional Gaussian noise in any given direction (f Gn), G(x, h) = Bϕ(x + h) − Bϕ(x). The Hurst coefficient describes the nature of these
increments and is bounded from 0 to 1. When H equals 0.5 then the increments are uncor-related and a positive K increment is equally likely to be followed by a negative or positive increment. When H>0.5 there is a positive correlation structure and when H<0.5 there is a negative correlation structure. A positive correlation structure means a positive K in-crement is more likely to be followed by a positive inin-crement and vice versa. A negative correlation structure means a K increment is more likely to be followed by an increment of the opposite sign. Dispersional analysis and rescaled range analysis are used to estimate H
in the longitudinal and transverse directions (Monnig et al., 2008). The slope of the rescaled range statistic versus partition size on a log log plot is equal to H (Monnig et al, 2008).
Y-Dir ec tion [m] X-Direction [m] Z -Dir ec tion [m]
Figure 3.3: Locations of the 9 DP measurements, domain boundary conditioning points, and distinguished GPR facies within the 50 m by 50 m model domain. The extent of the facies represents a 10 by 10 meter area within the ICA cube.
The unconditional fBm, Bϕ, is conditioned at each DP location to create a smoothed
trend through Bϕ denoted Bavg. The residuals of the fBm smoothed trend surface (Bavg)
and the unconditional fBm (Bϕ) are added to the initial trend surface (Kavg) to obtain the
final conditioned fBm K field (Figure 3.5, Equation 3.3).
Kout = Kavg + (Bϕ− Bavg) (3.3)
For further understanding, Figure 3.4 is a one-dimensional conceptual diagram of the process described above (Monnig, 2012).
0 200 400 600 800 1000 1200 −1 −0.5 0 0.5 1 1.5
Step 1: Interpolate Conditioning Points
Kavg 0 200 400 600 800 1000 1200 −1.5 −1 −0.5 0
Step 3: Interpolate fBm at Conditioning Points
fBm (Conditioning Points)
fBmavg
Step 1: Interpolate Conditioning Points
Step 3: Interpolate fBm at Cond. Points
0 200 400 600 800 1000 1200 −1.5 −1 −0.5 0 Step 2: Generate fBm fBm 0 200 400 600 800 1000 1200 −1 −0.5 0 0.5 1
1.5 Step 4: Generate Conditioned K−field
K = Kavg + fBm − fBmavg
Step 2: Generate fBm
Step 4: Generate Conditioned K-field
fBm K avg
itioning Points
fBm (Conditioning Points) fBmavglate fBm at Cond. Points
K = Kavg + fBm − fBmavg
Figure 3.4: Illustration of creating a fBm K-field in 1-dimension. The first step is creating a smooth trend of the original dataset. The second step is generating an unconditional fBm by convolving Gaussian noise and the filter kernel. The third step is generating a smooth trend through the unconditional fBm at the same locations as the actual K data. The fourth and final step is calculating the residuals of the unconditional fBm and the smoothed trend through the unconditional fBm and adding them to the original smoothed trend in Step 1. (Monnig, 2012).
Figure 3.5: osfBm hydraulic conductivity (K) field. This K field will be denoted f Bm K. The units of K are in meters per hour.
3.3 Upscaling
To investigate the ADE and t-fADE models at a range of hydraulic conductivity reso-lutions the fBm K field was upscaled. The original K domain (Figure 3.5) was upscaled to a progressively coarser grid by arithmetic averaging. Arithmetic averaging was used because groundwater flow is parallel to the subsurface layering. In many cases the scale of measurement is smaller than the numerical discretization size and the small-scale val-ues must be homogenized to fit this larger, regional scale (Wen and Gomez-Hernandez, 1996). The dimensions of the fBm domain are 128 x 128 x 512 cells in the x-, y-, and z- directions respectively. The larger number of grid blocks in the z-direction are be-cause of the high-resolution vertical DP measurements. The vadose zone was
incorpo-rated by adding a volume of 128 x 128 x 128 cells to the top of the existing fBm K do-main. The entire volume of the vadose zone has one saturated K (Ksat) value of 0.191 m
hr−1, thatisthearithmeticaverageof theoriginalf BmKf ield.T oupscale, f irstthefBmKf ieldwasarithmetical
0
1
2
3
4
5
fBM K
Figure 3.6: Slice through each upscaled hydraulic conductivity (K) domain from f bm K to Homogeneous K (5). This illustration is to visualize the size of each grid block rather than the K value of each grid block so a color legend is not provided.
3.4 Fine K Realizations
The f Bm K field shown in Figure 3.5 is one of an infinite number of realizations. To compare the variability in f Bm K fields with similar statistics (mean and variance) multiple fBm simulations were performed to create new realizations. The reason for generating many K field realizations is to compare BTC variability from equally probable K fields in a highly heterogeneous aquifer. Table 3.2 contains the mean, variance, and standard deviation of ln(K) for all fine realizations. These were calculated in the area of interest (center of the K
domain and area of SWIW test). Only the center of the domain is important for comparison between realizations because there are few conditioning points at the boundaries of the model domain. Few conditioning points provide less restraint on predicted values and so the K values vary widely in these areas. The boundaries of the domain are far from the area of interest and are designed to have minimal impact in the approximately 10m by 10m region of the SWIW test.
Table 3.2: The mean, variance, and standard deviation of ln(K) for each fractional Brownian motion (f Bm) hydraulic conductivity (K) field realization.
Realization Mean (ln(K))(m/hr) Variance ln(K) Stand. Dev. ln(K)
Original fBm −0.6532 4.2669 2.0656 1 −0.4924 6.355 2.5209 2 −0.322 11.0809 3.3288 3 −0.7339 11.0102 3.3182 4 −1.1473 6.703 2.589 5 −0.2084 8.6635 2.9434 6 −0.4594 9.0306 3.0051 7 −1.1145 8.8416 2.9735 8 −0.2019 7.0064 2.647 9 −0.4038 7.9964 2.8278 10 −1.0124 9.7562 3.1235 11 −0.8926 9.7174 3.1173 12 0.0686 7.9397 2.8177 13 −0.2646 6.1812 2.4862 14 −0.3834 8.0724 2.8412
A forward run to simulate the SWIW BTC will be performed for all f Bm realizations. The same parameter values will be used so the variability in BTCs on realization to
realiza-tion basis can be assessed.
3.5 Model Simulations
To run a full simulation, first Parflow was used to generate groundwater velocity, pressure, and saturation files for each domain (f Bm K through HomogeneousK). These files are then needed by the transport code, Slim-FAST, for particle tracking and obtaining a final concentration BTC. Finally the simulated BTCs are compared to the measured BTC within the parameter estimation code, PEST, to obtain best-fit parameter values.
3.5.1 Parflow
Parflow was used to calculate groundwater velocities in variably-saturated conditions. Parflow is a watershed model that can simulate surface and subsurface fluid flow in steady-state saturated, variably saturated, and integrated-watershed flow conditions (Maxwell et al., 2008). To generate 3D groundwater velocities for each domain, hydraulic conductivity (K), porosity (φ), permeability (k), specific storage (SS), and van Genuchten parameters
α, n, θr, and θs were needed (Appendix B). Van Genuchten parameters were selected from
Carsel and Parrish (1998). The van Genuchten parameters chosen are for loamy sand because the Ksat value for loamy sand was closest to the arithmetic average Ksat value of 0.191 (m
hr−1)(T able 3.3).T hisK
sat value was chosen for all of the vadose zone. A porosity 0.32 has
been used in many previous studies and was first estimated by Rehfeldt et al. (1992) and a specific storage value of 10−6 (m−1) was used. The unsaturated zone porosity was later
optimized to estimate transport in multiphase flow.
Table 3.3: Parameter values needed for Parflow simulations
Parameter Value
Porosity, φ[−] 0.32 Specific Storage, Ss[L−1] 10−6
α[L−1] 12.4
n[−] 2.28
The model domain is 50 m by 50 m with a 0.390625 m grid resolution in the x and y-directions. The z-direction has a grid resolution of 0.016797 m to a depth of 10.75 m. The following boundary conditions were used: (1) top and bottom of the domain were no flow boundaries, (2) sides of the domain were set to constant head boundaries using the water table as a reference. The well is located near the center of the domain at 25.1953125 m in the x- and y-directions and from 2.4 to 8 m in the z direction. The total run time is 460 hours and pressure, saturation, and head files were output each hour.
3.5.2 Slim-FAST
Slim-FAST is a particle transport code that uses a Lagrangian approach to track particles throughout the numerical simulation (Maxwell et al., 2010). Each particle is assigned a position, velocity, mass, and diffusivity, all of which are updated through time. Slim-FAST requires velocity files that can be derived from Parflow pressure files as well as PFB files of the van Genuchten parameters. Two versions of Slim-FAST were used; transport simulation code written by Maxwell (2010) for the ADE and a modified version updated by Dean (2011) for the t-fADE. Inject, wait, and extract periods of the SWIW test are run separately for 31, 19, and 1000 hours respectively (Appendix B). All particles begin in the mobile phase and transition into the immobile phase after an assigned random amount of time based on an exponential distribution. When a particle enters the immobile phase the time it remains there is based on a power law distribution. This power law distribution is the result of the memory function described in Section 2.3.1. Particle input during the injection period is flux weighted dependent upon the radial velocity adjacent to the well. Due to water table mounding during the injection period the screened interval is set from 2.4 m to 8.5 m. The screened interval during the wait and extract periods is 2.4 m to 8 m as the water table lowers. At the end of each period the x-, y-, and z- coordinates of each particle are recorded in a text file that then becomes the particle input file for the following time period. During the extraction period concentrations were calculated per hour.
3.5.3 Parameter Estimation and Optimization: PEST
The parameter estimation code, PEST, was used to determine optimized parameter val-ues. For the ADE model porosity (φ), longitudinal dispersivity (αL), and transverse
disper-sivity (αT) were optimized. For the t-fADE model porosity (φ), longitudinal dispersivity
(αL), transverse dispersivity (αT), and the capacity coefficient (β) were optimized. The
ad-ditional parameter in the t-fADE equation, γ, was manually fit by determining the slope of the heavy-tail on a log-log plot. For all simulations a γ value of 0.64 was used. A total of 410 measured concentrations, one for each extraction hour, were calculated with a moving aver-age and smoother function to attain the final simulated BTC (Appendix B). The smoother function treats each particle as a pseudo-Gaussian distribution and the BTC as a convolu-tion of them (Pedretti et al., 2013). Residuals were derived for each measured concentraconvolu-tion. These residuals were then squared, weighted by the reciprocal of the measured concentration (1/mi), and summed to calculate the weighted sum of square residuals (W SSR) objective
function: W SSR = n X i=1 wi(mi− si)2 (3.4)
where w is the weight, m is the measured concentration, s is the simulated concentration, and n is the total number of observations. The weights were chosen because they emphasize low concentrations over larger concentrations, putting larger importance in the heavy-tailing of the BTC rather than the early-time peak. PEST continually ran the transport code until a minimum or convergence of the WSSR value is met. The objective function quantifies error by measuring the accuracy of model fit to observation and when this value reaches a minimum the current parameter values are considered the actual. PEST also calculates parameter sensitivities as well as various model comparison statistics. All PEST input files are located in Appendix B.
3.6 Error Statistics and Model Comparison
Error was quantified for both models at each scale of measurement. When considering the suitability of a model to simulate a system one must not only consider the model’s abil-ity to accurately predict the physical system, but also consider the number of parameters needed and the amount of data required. The tradeoff between these three aspects is con-sidered model fitness which is the combination of simulated to measured fit as well as model parsimony. A model is considered to be parsimonious if it has the ability to accurately pre-dict a system without the need for many parameters or highly detailed data. In this study, a parsimonious model would be one that could predict fairly well the BTC concentrations with less detailed velocity information and fewer parameters. Although the t-fADE has more parameters than the ADE if it better matches the observed data with far less K information then it would be considered more parsimonious.
CHAPTER 4
RESULTS AND DISCUSSION
The optimized parameters from the ADE and t-fADE models were obtained for seven different hydraulic conductivity domains, totaling 14 separate optimization runs. These are summarized in Table 4.1. There is one fine-scaled K domain, f Bm K, and 6 upscaled K domains labeled K0 to Homogeneous K. All of the upscaled domains are averages of the f Bm K domain. The BTCs and optimized parameters for each upscaled domain will be analyzed separately and then compared. At the finest scale, 14 additional fBm K field realizations were created and BTCs were simulated for both models using the original f Bm K’s optimized parameter values. The resulting BTCs will be evaluated to distinguish how similar results are from equally probable K fields.
Table 4.1: Summary of Upscaling model runs, their Weighted Sum of Squared Residuals (WSSR), and % Mass Recovery.
K Domains [x, y, z] Model Simulation WSSR Mass Recovered %
f Bm K ADE 1 5.59 76.27 [128 × 128 × 512] t-fADE 2 4.59 75.95 K0 ADE 3 4.40 77.08 [32 × 32 × 32] t-fADE 4 3.45 76.87 K1 ADE 5 11.44 81.84 [16 × 16 × 16] t-fADE 6 6.86 78.15 K2 ADE 7 16.46 92.95 [8 × 8 × 8] t-fADE 8 13.05 88.21 K3 ADE 9 22.60 95.60 [4 × 4 × 4] t-fADE 10 15.72 85.20 K4 ADE 11 22.09 95.94 [2 × 2 × 2] t-fADE 12 14.89 92.53 Homogeneous K ADE 13 27.58 95.97 [1 × 1 × 1] t-fADE 14 11.84 91.72 4.1 Upscaling
Specifics for each upscaled K domain are listed in ??. Although the main goals of this section are to analyze final BTCs and optimized parameter values it is important to visualize
how the complexity of the hydraulic conductivity field affects particle movement. Plotting particle location through time from a slice of the model domain displays the differences in transport that take place dependent upon K resolution. Figure 4.1 shows time-lapse particle movement during the injection period of the finest (f Bm K), middle (K2), and coarsest (Homogeneous K) domains.
Y-Direction [m] Y-Direction [m] Y-Direction [m] Z -Dir ec tion [m] Z -Dir ec tion [m] Z -Dir ec tion [m] time 31 time 20 time 10 time 6 time 4 time 1 a. 31 20 10 6 4 1 c. b.
Figure 4.1: Time-lapse of particle injection, from left to right: a) finest (f Bm) K field, b) middle (K2), and c) coarsest (Homogeneous) K field. Particle locations are recorded at hours 1, 4, 6, 10, 20, and 31.
The particles in the Homogeneous K domain (Figure 4.1c) are injected and move radially around the well as a block. Transport properties are vertically uniform which is expected in a single material. Although the domain consists of one K value, there is evidence of dispersion around the edges of the plume. If dispersion were set to zero then there would be a sharp contact at the edge of the plume rather than the fuzzy look seen in Figure 4.1c. The dotted black line on each plot represents the top of the screened interval (8.5 m). For the homogeneous domain particles just go over the water table into the vadose zone. The time-lapse image of the f Bm K domain (Figure 4.1a) shows channels of high K zones as well as areas where particles do not flow as quickly. The higher K values are primarily concentrated towards the top of the screened interval whereas the lower sections of the screened interval have lower K values. Particle movement around the well is not symmetrical
like the homogenous K domain. Significant water table mounding occurs during the injection period in the fBm K domain transporting many particles above the water table into the vadose zone. The influence of the flux-weighted particle input is also evident in Figure 4.1. The number of particles injected vertically along the screened interval is dependent upon the radial velocities adjacent to the well. The homogeneous domain has an equal number of particles injected along the screen because the velocities are equivalent. Figure 4.1c has a solid green bar where the well is located. The velocities along the screened interval of the f Bm K domain vary and in some areas are so low no particles enter the well in the corresponding zone. This is evident in Figure 4.1a by looking at the color corresponding to Time 1. The medium (K2) upscaled K domain (Figure 4.1b) has sections of higher K but is very coarse compared to the f Bm K field. Some characteristics of the f Bm K domain can be seen, but the decrease in K complexity is evident. The extent of particle movement into the aquifer and the number of particles that have entered the vadose zone varies greatly dependent upon K resolution. Of particular importance are particles that enter the vadose zone because these can become trapped and are not captured during the extraction period resulting in low mass recovery.
4.1.1 Upscaling of the ADE Model
BTC figures for all upscaled K domains for both the ADE and t-fADE models will be represented with a cold to hot color scheme. The finest K domain, f Bm K, is blue and the coarsest K domain, the HomogeneousK, is red. Legends on the upscaled BTC figures (Figure 4.2 and Figure 4.4) also depict this.
Breakthrough Curves. The simulated BTCs as the K domain coarsens vary greatly for the ADE model (Figure 4.2). K1 through the Homogeneous K over-predict the peak concentration. The peak concentrations do not decrease progressively as the K field fines as expected. An example of this is the coarsest, Homogenous K domain (color: dark red). The Homogeneous K does not have the highest peak concentration, but is still the worst
