Modeling of non-reactive solute transport in fractured clayey till during variable flow rate and time

Fractures and biopores can act as preferential flow paths in clay aquitards and may rapidly transmit contaminants into underlying aquifers. Reliable numerical models for assessment of groundwater contamination from such aquitards are needed for planning, regulatory and remediation purposes. In this investigation, high resolution preferential water-saturated flow and bromide transport data were used to evaluate the suitability of equivalent porous medium (EPM), dual porosity (DP) and discrete fracture/matrix diffusion (DFMD) numerical modeling approaches for assessment of flow and non-reactive solute transport in clayey till. The experimental data were obtained from four large undisturbed soil columns (taken from 1.5 to 3.5 m depth) in which biopores and channels along fractures controlled 96-99% of water-saturated flow. Simulating the transport data with the EPM effective porosity model (FRACTRAN in EPM mode) was not successful because calibrated effective porosity for the same column had to be varied up to 1 order of magnitude in order to simulate solute breakthrough for the applied flow rates between 11 and 49 mm/day. Attempts to simulate the same data with the DP models CXTFIT and MODFLOW/MT3D were also unsuccessful because fitted values for dispersion, mobile zone porosity, and mass transfer coefficient between mobile and immobile zones varied several orders of magnitude for the different flow rates, and because dispersion values were furthermore not physically realistic. Only the DFMD modeling approach (FRACTRAN in DFMD mode) was capable to simulate the observed changes in solute transport behavior during alternating flow rate without changing values of calibrated fracture spacing and fracture aperture to represent the macropores.     

Sorption and mobility of dithiopyr in golf course greens rooting medium

Sorption and mobility of dithiopyr in golf course greens rooting medium (RM) were studied. The sorption increased from 20 to 27 degrees C at 24 h after treatment (HAT) and reached equilibrium in 48 HAT at 20 degrees C. The sorption isotherms had Freundlich values (KF) of 1122, 27.44 to 35.16, and 0.053 to 0.168 for peat moss, the RM, and quartz sand, respectively, and solid to aqueous phase partition coefficients (Kd) of 470 to 1706 L/kg, 14.61 to 84.4 L/kg, and 0.07 to 0.29 L/kg for peat moss, RM, and quartz sand, respectively. Generally, higher dithiopyr concentration in the aqueous solution and the reduced pH of the solution corresponded to the higher Kd values. The average values for dispersion (D, cm2/min), retardation coefficient (R), beta, and omega parameters for solute transport in the RM lysimeter; obtained from CXTFIT curve fitting of Br- breakthrough curves; were 0.95, 1.01, 1, and 93.89, respectively. After elutriation by 18 L of aqueous KNO3 (10 mM), greater than 90% of the added dithiopyr remained in the top 10 cm of the RM lysimeter and no detectable dithiopyr was present at depths beyond 35 cm. The lysimeter effluent contained dithiopyr at concentrations less than 3.5 microg/L. The R value obtained from CXFIT curve fitting is 38.5. Results from both sorption and mobility experiments indicated that dithiopyr is quite immobile in golf course greens RM and has minimal potential for movement into surface water drainage or ground water.     

Continuous time random walk model better describes the tailing of atrazine transport in soil

Contaminant transport in soils is complicated and involves some physical and chemical nonequilibrium processes. In this research, the soil column displacement experiments of Cl⁻ and atrazine under different flow velocities were carried out. The data sets of Cl⁻ transport in sandy loam fitted to the convection dispersion equation (CDE) and the two-region model (TRM) indicated that the effects of physical nonequilibrium process produced by immobile water on the breakthrough curves (BTCs) of Cl⁻ and atrazine transport through the repacking soil columns were negligible. The two-site model (TSM) and the continuous time random walk (CTRW) were also used to fit atrazine transport behavior at the flow rate of 19.86 cm h⁻¹. The CTRW convincingly captured the full evolution of atrazine BTC in the soil column, especially for the part of long tailing. However, the TSM failed to characterize the tailing of atrazine BTC in the soil column. The calculated fraction of equilibrium sorption sites, F, ranging from 0.78 to 0.80 for all flow rates suggested the contribution of nonequilibrium sorption sites to the asymmetry of atrazine BTCs. Furthermore, the data sets for the flow rates of 6.68 cm h⁻¹ and 32.81 cm h⁻¹ were predicted by the TSM and the CTRW. As to the flow rate of 6.68 cm h⁻¹, the CTRW predicted the entire BTC of atrazine transport better than the TSM did. For the flow rate of 32.81 cm h⁻¹, the CTRW characterized the late part of the tail better, while the TSM failed to predict the tailings of atrazine BTC.     

An interpretation of potential scale dependence of the effective matrix diffusion coefficient

Matrix diffusion is an important process for solute transport in fractured rock, and the matrix diffusion coefficient is a key parameter for describing this process. Previous studies have indicated that the effective matrix diffusion coefficient values, obtained from a large number of field tracer tests, are enhanced in comparison with local values and may increase with test scale. In this study, we have performed numerical experiments to investigate potential mechanisms behind possible scale-dependent behavior. The focus of the experiments is on solute transport in flow paths having geometries consistent with percolation theories and characterized by multiple local flow loops formed mainly by small-scale fractures. The water velocity distribution through a flow path was determined using discrete fracture network flow simulations, and solute transport was calculated using a previously derived impulse-response function and a particle-tracking scheme. Values for effective (or up-scaled) transport parameters were obtained by matching breakthrough curves from numerical experiments with an analytical solution for solute transport along a single fracture. Results indicate that a combination of local flow loops and the associated matrix diffusion process, together with scaling properties in flow path geometry, seems to be an important mechanism causing the observed scale dependence of the effective matrix diffusion coefficient (at a range of scales).     

Increasing the efficiency of solute leaching: impacts of flow interruption with drainage of the “preferential flow paths”

Most soils contain preferential flow paths that can impact on solute mobility. Solutes can move rapidly down the preferential flow paths with high pore-water velocities, but can be held in the less permeable region of the soil matrix with low pore-water velocities, thereby reducing the efficiency of leaching. In this study, we conducted leaching experiments with interruption of the flow and drainage of the main flow paths to assess the efficiency of this type of leaching. We compared our experimental results to a simple analytical model, which predicts the influence of the variations in concentration gradients within a single spherical aggregate (SSA) surrounded by preferential flow paths on leaching. We used large (length: 300 mm, diameter: 216 mm) undisturbed field soil cores from two contrasting soil types. To carry out intermittent leaching experiments, the field soil cores were first saturated with tracer solution (CaBr₂), and background solution (CaCl₂) was applied to mimic a leaching event. The cores were then drained at 25- to 30-cm suction to empty the main flow paths to mimic a dry period during which solutes could redistribute within the undrained region. We also conducted continuous leaching experiments to assess the impact of the dry periods on the efficiency of leaching. The flow interruptions with drainage enhanced leaching by 10–20% for our soils, which was consistent with the model's prediction, given an optimised “equivalent aggregate radius” for each soil. This parameter quantifies the time scales that characterise diffusion within the undrained region of the soil, and allows us to calculate the duration of the leaching events and interruption periods that would lead to more efficient leaching. Application of these methodologies will aid development of strategies for improving management of chemicals in soils, needed in managing salts in soils, in improving fertiliser efficiency, and in reclaiming contaminated soils.     

Evaluating non-equilibrium solute transport in small soil columns

Displacement studies on leaching of bromide and two pesticides (atrazine and isoproturon) were conducted under unsaturated steady state flow conditions in 24 small undisturbed soil columns (5.7 cm in diameter and 10 cm long) each collected from two sites differing in soil structure and organic carbon content in North Germany. There were large and irregular variabilities in the characteristics of both soils, as well as in the shapes of breakthrough curves (BTCs) of different columns, including some with early breakthrough and increased tailing, qualitatively indicating the presence of preferential flow. It was estimated that one preferential flow column (PFC) at site A, and four at site B, contributed, respectively to 11% and 58% of the accumulated leached fraction and to more than 80% of the maximum observed standard deviation (SD) in the field-scale concentration and mass flux of pesticides at two sites. The bromide BTCs of two sites were analyzed with the equilibrium convection-dispersion equation (CDE) and a non-equilibrium two-region/mobile-immobile model. Transport parameters of these models for individual BTCs were determined using a curve fitting program, CXTFIT, and by the time moment method. For the CDE based equilibrium model, the mean values of retardation factor, R, considered separately for all columns, PFCs or non-preferential flow columns (NPFCs) were comparable for the two methods; significant differences were observed in the values of dispersion coefficients of two sites using the two estimation methods. It was inferred from the estimated parameters of non-equilibrium model that 5-12% of water at site A, and 12% at site B, was immobile during displacement in NPFCs. The corresponding values for PFCs of two sites were much larger, ranging from 25% to 51% by CXTFIT and from 24% to 72% by the moment method, suggesting the role of certain mechanisms other than immobile water in higher degrees of non-equilibrium in these columns. Peclet numbers in PFCs of both sites were consistently smaller than five, indicating the inadequacy of the non-equilibrium model to incorporate the effect of all forms of non-equilibrium in PFCs. Overall, the BTCs of individual NPFCs, PFCs and of field average concentration at the two sites were better reproduced with parameters obtained from CXTFIT than by the moment method. The moment method failed to capture the peak concentrations in PFCs, but tended to describe the desorption and tail branches of BTCs better than the curve fitting approach.     

Fate of nitrogen for subsurface drip dispersal of effluent from small wastewater systems

Subsurface drip irrigation systems apply effluent from onsite wastewater systems in a more uniform manner at a lower rate than has been possible with other effluent dispersal methods. The effluent is dispersed in a biologically active part of the soil profile for optimal treatment and where the water and nutrients can be utilized by landscape plants. Container tests were performed to determine the fate of water and nitrogen compounds applied to packed loamy sand, sandy loam, and silt loam soils. Nitrogen removal rates measured in the container tests ranged from 63 to 95% despite relatively low levels of available carbon. A Hydrus 2D vadose zone model with nitrification and denitrification rate coefficients calculated as a function of soil moisture content fit the container test results reasonably well. Model results were sensitive to the denitrification rate moisture content function. Two-phase transport parameters were needed to model the preferential flow conditions in the finer soils. Applying the model to generic soil types, the greatest nitrogen losses (30 to 70%) were predicted for medium to fine texture soils and soils with restrictive layers or capillary breaks. The slow transport with subsurface drip irrigation enhanced total nitrogen losses and plant nitrogen uptake opportunity.     

Competitive sorption and diffusion of chromate and sulphate in a flow system with goethite in gel beads

Column experiments and model simulations were employed to evaluate the processes involved in multicomponent solute transport in a system with heterogeneous flow. Column experiments were performed with goethite embedded in polyacrylamide gel beads. The gel forms an immobile water region that can be accessed by diffusion. A two-region transport model with diffusion into spheres was combined with a surface complexation model to predict reactive transport in the goethite-gel bead system. Chromate and sulphate breakthrough curves were measured in a set of transport experiments, along with corresponding changes in the pH of the effluent. Sorption and transport of sulphate and chromate in separate columns were predicted from independently measured sorption parameters. The model overestimated the pH changes in the effluent, possibly because of proton buffering by the polyacrylamide gel. The effect of competitive sorption on transport was examined in experiments with both anions present. The model predicted the effect of competition very well in a system initially equilibrated with sulphate, followed by infiltration with chromate. However, when sulphate was infiltrated after equilibration with chromate, chromate desorption and sulphate adsorption were clearly overestimated by the transport model. The exchange between the more strongly bound chromate and the sulphate added subsequently may be too slow to cause a substantial chromate peak in the effluent. This suggests that the local equilibrium assumption was not applicable in this case.     

An active region model for capturing fractal flow patterns in unsaturated soils: model development

Preferential flow commonly observed in unsaturated soils allows rapid movement of solute from the soil surface or vadose zone to the groundwater, bypassing a significant volume of unsaturated soil and increasing the risk of groundwater contamination. A variety of evidence indicates that complex preferential patterns observed from fields are fractals. In this study, we developed a relatively simple active region model to incorporate the fractal flow pattern into the continuum approach. In the model, the flow domain is divided into active and inactive regions. Flow occurs preferentially in the active region (characterized by fractals), and inactive region is simply bypassed. A new constitutive relationship (the portion of the active region as a function of saturation) was derived. The validity of the proposed model is demonstrated by the consistency between field observations and the new constitutive relationship.     

Simulations of physical nonequilibrium solute transport models: Application to a large-scale field experiment

Equations expressing the spatial moments of solute concentration distributions simulated by various models, in terms of model parameters, have recently been presented. Using independently obtained parameter values, these equations are used to compare simulations of physical non-equilibrium models with spatial moment data collected in a large-scale natural gradient experiment on solute transport. The physical nonequilibrium models examined postulate the existence of layered zones of immobile water through which solute is transported by a diffusion mechanism. It is found that the qualitative aspects of the measured moment behavior are simulated by the physical nonequilibrium models if the independently obtained parameters are modified somewhat on the basis of reasonable corrective assumptions. It is further demonstrated that the physical nonequilibrium models, using parameter values obtained from spatial data, can qualitatively simulate temporal behavior at individual well points in this relatively homogeneous aquifer.     

Modelling of atrazine transport in the presence of surfactants

Laboratory experiments were conducted to examine the effect of detergents on transport of atrazine [2-chloro-4-(ethylamino)-6-(isopropylamino)-s-triazine] through loam and sandy loam soils under saturation conditions. The Convection Dispersion Equation (CDE) was used to model and quantify the effects of detergents on atrazine model parameters: the retardation factor (R), pore velocity (v) and dispersion coefficient (D). The transport parameters were estimated using moment technique and partition coefficient obtained from batch experiments and compared with best-fitted parameters, R and D, keeping pore velocity constant. Results indicated the CDE model was not successful in predicting atrazine transport in the presence of surfactants at high concentrations. In the case of anionic surfactant with Elora loam, the average predicted R and D from moment technique of 3.4 and 11.1 cm2/h, respectively were significantly different than fitted parameters (R = 39 and D = 227 cm2/h). The poor performance of CDE in the presence of surfactants results from physiochemical changes in herbicide solubility and retention to the soil matrix rather than changes in soil hydraulic properties since the predicted pore water velocities from moment technique were similar to those measured during leaching experiments. Nevertheless, BTC analysis with CDE showed that land application of anionic surfactant (sulphonic) significantly increased R and D and decrease v for both soils. Addition of sulphonic increased R of atrazine by 12 and 26 folds for loam and sandy loam soils, respectively. On the other hand non-ionic surfactants seemed to decrease R, especially in sandy loam soil, thus facilitating atrazine leaching through soil. Non-equilibrium conditions seemed to govern atrazine transport in the presence of surfactants; double peaks in breakthrough curves were observed, indicating a need for mathematical models to account for such phenomena. Atrazine dispersion and tailing seemed to be higher through Elora loam compared to Caledon sandy loam due to higher aggregation of the Elora soil.     

Transport modeling of copper and cadmium with linear and nonlinear retardation factors

The predictive accuracy of using the one-dimensional advection–dispersion equation to evaluate the fate and transport of solute in a soil column is usually dependent on the proper determination of chemical retardation factors. Typically, the distribution coefficient (K_d) obtained by fitting the linear sorption isotherm has been extensively used to consider general geochemical reactions on solute transport in a low-concentration range. However, the linear distribution coefficient cannot be adequately utilized to describe the solute fate at a higher concentration level. This study employed the nonlinear equilibrium-controlled sorption parameters to determine the retardation factor used in column leaching experiments. Copper and cadmium transportation in a lateritic silty-clay soil column was examined. Through the explicit finite-difference calculations with a third-order total-variation-diminishing (TVD) numerical solution scheme, all results of the theoretical copper and cadmium breakthrough curves (BTCs) simulated by using the Freundlich nonlinear retardation factors revealed good agreement with the experimental observations.     

Investigation of hydrogeologic processes in a dipping layer structure: 2. Transport and biodegradation of organics

Numerical simulation tools have been used to study the dominating processes during transport of aromatic hydrocarbons in the unsaturated soil zone. Simulations were based on field observations at an experimental site located on a glacial delta plain with pronounced layered sedimentary structures. A numerical model for transport in the unsaturated zone, SWMS-3D, has been extended to incorporate coupled multispecies transport, microbial degradation following Monod kinetics and gas diffusive transport of oxygen and hydrocarbons. The flow field parameters were derived from previous work using nonreactive tracers. Breakthrough curves (BTC) from the hydrocarbon field experiment were used to determine sorption parameters and Monod kinetic parameters using a fitting procedure. The numerical simulations revealed that the assumption of homogeneous layers resulted in deviations from the field observations. The deviations were more pronounced with incorporation of reactive transport, compared with earlier work on nonreactive transport. To be able to model reasonable BTC, sorption had to be reduced compared to laboratory experiments. The initial biomass and the maximum utilisation rate could be adjusted to capture both the initial lag phase and the overall degradation rate. Nevertheless, local oxygen limitation is predicted by the model, which was not observed in the field experiment. Incorporation of evaporation and diffusive gas transport of the hydrocarbons did not significantly change the local oxygen demand. The main cause of the observed discrepancies between model and field are attributed to channelling as a result of small-scale heterogeneities such as biopores.     

Influence of reactor configuration on reliability of activated sludge process

The effect of uncertainty in system parameters and the accuracy of the mathematical model used on the design and reliability of the activated sludge process is investigated by Monte Carlo simulations. Simulations indicate that the coefficient of variation for a reactor volume varies from 0.56 for a single mixed tank reactor to 0.44 for a plug flow reactor. The coefficient of variation for effluent for reactors designed on nominal values was found to be 0.56 for a mixed-tank reactor, 1.28 for two reactors in series, 1.56 for three reactors in series and 1.6 for a plug flow reactor. Significant contributing parameters to the reliability of the process are established. Reactor volumes for desired reliability levels are also calculated.     

An independent prediction of the effect of dissolved organic matter on the transport of polycyclic aromatic hydrocarbons

The transport of polycyclic aromatic hydrocarbons (PAH) in porous media in the presence of dissolved organic matter (DOM) was predicted with a transport bicontinuum model using independently obtained relationships to derive transport parameters for describing the effect of PAH binding to the DOM. The sorption constants of PAHs to soil and their binding constants to DOM were derived from basic correlations with K(ow) (indicator of hydrophobicity). The kinetic (rate) constants were derived from previously published correlations with K(p) (sorption constant). The independently obtained sorption and rate constants were corrected for binding to DOM and were used to predict the breakthrough curves (BTC) of contaminants in the presence and the absence of DOM. Column results confirmed the independently predicted BTC of PAHs in the presence of DOM that did not sorb to the solid phase, as well as the effect of DOM on the rate of the sorption and desorption processes. These findings confirm the ability to quantitatively describe how DOM facilitates transport of contaminants in the subsurface using independently derived parameters.     

Modelling the solute transport under nonequilibrium conditions on the basis of mass transfer equations

A solute transport model that describes nonequilibrium adsorption in soil/groundwater systems by mass transfer equations for film and intraparticle diffusion is presented. The model is useful in cases where breakthrough curve spreading cannot be explained by dispersion only. To evaluate its validity, the model was applied to several data sets from column experiments. The validity was also proved by a comparison with an analytical solution for the limiting case of predominating dispersion. Furthermore, a sensitivity analysis was performed to illustrate the influence of different process and sorption parameters (pore water velocity, intraparticle mass transfer coefficient, isotherm nonlinearity) on the shape of the calculated breakthrough curves. The application of the proposed model is discussed in comparison to the widely used dispersed flow/local equilibrium model, and a relationship between both models, which is based on a lumped parameter approach, is shown.     

Interpretation of injection-withdrawal tracer experiments conducted between two wells in a large single fracture

Tracer experiments conducted using a flow field established by injecting water into one borehole and withdrawing water from another are often used to establish connections and investigate dispersion in fractured rock. As a result of uncertainty in the uniqueness of existing models used for interpretation, this method has not been widely used to investigate more general transport processes including matrix diffusion or advective solute exchange between mobile and immobile zones of fluid. To explore the utility of the injection-withdrawal method as a general investigative tool and with the intent to resolve the transport processes in a discrete fracture, two tracer experiments were conducted using the injection-withdrawal configuration. The experiments were conducted in a fracture which has a large aperture (>500 microm) and horizontally pervades a dolostone formation. One experiment was conducted in the direction of the hydraulic gradient and the other in the direction opposite to the natural gradient. Two tracers having significantly different values of the free-water diffusion coefficient were used. To interpret the experiments, a hybrid numerical-analytical model was developed which accounts for the arcuate shape of the flow field, advection-dispersion in the fracture, diffusion into the matrix adjacent to the fracture, and the presence of natural flow in the fracture. The model was verified by comparison to a fully analytical solution and to a well-known finite-element model. Interpretation of the tracer experiments showed that when only one tracer, advection-dispersion, and matrix diffusion are considered, non-unique results were obtained. However, by using multiple tracers and by accounting for the presence of natural flow in the fracture, unique interpretations were obtained in which a single value of matrix porosity was estimated from the results of both experiments. The estimate of porosity agrees well with independent measurements of porosity obtained from core samples. This suggests that: (i) the injection-withdrawal method is a viable tool for the investigation of general transport processes provided all relevant experimental conditions are considered and multiple conservative tracers are used; and (ii) for the conditions of the experiments conducted in this study, the dominant mechanism for exchange of solute between the fracture and surrounding medium is matrix diffusion.     

Modeling solute transport in one-dimensional homogeneous and heterogeneous soil columns with continuous time random walk

In this paper, we used the continuous time random walk (CTRW) framework to characterize the transport process in 1250-cm long one-dimensional homogenous and heterogeneous soil columns at the experiments conducted by Huang et al. [Huang, K., Toride, N., van Genuchten, M.Th., 1995. Experimental investigation of solute transport in large, homogeneous and heterogeneous, saturated soil columns. Trans. Porous Media. 18, 283-302]. The transport process was also simulated by using the advection-dispersion equation (ADE) and the spatial fractional advection-dispersion equation (FADE) for comparison. In the homogeneous soil column, the non-Fickian behavior is found at the distances less than 1000cm with beta values larger than 1.60, but less than 2, and Fickian form transport is obtained at distances larger than 1000cm with beta values larger than 2. In the heterogeneous soil column, we found the most anomalous behavior at distances from 200cm to 700cm with beta values ranging from 0.894 to 0.958, and non-Fickian transport process is observed at distances larger than 800cm with beta values in the range between 1 and 1.3. More significant non-Fickian behavior is found for transport in the heterogeneous soil column than that in the homogeneous soil column. The CTRW fits to the breakthrough curves (BTCs) have lower values of root mean square error (RMSE) and higher values of determination coefficient (r(2)), with respect to the fits of ADE and FADE. The CTRW model also is better captures the full evolution of BTCs, and especially their tails.