Transport in chemically and mechanically heterogeneous porous media: V. Two-equation model for solute transport with adsorption

In this paper we develop the two-equation model for solute transport and adsorption in a two-region model of a mechanically and chemically heterogeneous porous medium. The closure problem is derived and the coefficients in both the one- and two-equation models are determined on the basis of the Darcy-scale parameters. Numerical experiments are carried out for a stratified system at the aquifer scale, and the results are compared with the one-equation model presented in Part IV and the two-equation model developed in this paper. Good agreement between the two-equation model and the numerical experiments is obtained. In addition, the two-equation model is used, in conjunction with a moment analysis, to derive a one-equation, non-equilibrium model that is valid in the asymptotic regime. Numerical results are used to identify the asymptotic regime for the one-equation, non-equilibrium model.     

Transport in chemically and mechanically heterogeneous porous media IV: large-scale mass equilibrium for solute transport with adsorption

In this article we consider the transport of an adsorbing solute in a two-region model of a chemically and mechanically heterogeneous porous medium when the condition of large-scale mechanical equilibrium is valid. Under these circumstances, a one-equation model can be used to predict the large-scale averaged velocity, but a two-equation model may be required to predict the regional velocities that are needed to accurately describe the solute transport process. If the condition of large-scale mass equilibrium is valid, the solute transport process can be represented in terms of a one-equation model and the analysis is simplified greatly. The constraints associated with the condition of large-scale mass equilibrium are developed, and when these constraints are satisfied the mass transport process can be described in terms of the large-scale average velocity, an average adsorption isotherm, and a single large-scale dispersion tensor. When the condition of large-scale mass equilibrium is not valid, two equations are required to describe the mass transfer process, and these two equations contain two adsorption isotherms, two dispersion tensors, and an exchange coefficient. The extension of the analysis to multi-region models is straight forward but tedious.     

Modeling non-equilibrium mass transport in biologically reactive porous media

We develop a one-equation non-equilibrium model to describe the Darcy-scale transport of a solute undergoing biodegradation in porous media. Most of the mathematical models that describe the macroscale transport in such systems have been developed intuitively on the basis of simple conceptual schemes. There are two problems with such a heuristic analysis. First, it is unclear how much information these models are able to capture; that is, it is not clear what the model's domain of validity is. Second, there is no obvious connection between the macroscale effective parameters and the microscopic processes and parameters. As an alternative, a number of upscaling techniques have been developed to derive the appropriate macroscale equations that are used to describe mass transport and reactions in multiphase media. These approaches have been adapted to the problem of biodegradation in porous media with biofilms, but most of the work has focused on systems that are restricted to small concentration gradients at the microscale. This assumption, referred to as the local mass equilibrium approximation, generally has constraints that are overly restrictive. In this article, we devise a model that does not require the assumption of local mass equilibrium to be valid. In this approach, one instead requires only that, at sufficiently long times, anomalous behaviors of the third and higher spatial moments can be neglected; this, in turn, implies that the macroscopic model is well represented by a convection–dispersion–reaction type equation. This strategy is very much in the spirit of the developments for Taylor dispersion presented by Aris (1956). On the basis of our numerical results, we carefully describe the domain of validity of the model and show that the time-asymptotic constraint may be adhered to even for systems that are not at local mass equilibrium.     

The impact of parameter lumping and geometric simplification in modelling runoff and erosion in the shrublands of southeast Arizona

There have been many studies of hydrologic processes and scale. However, some researchers have found that predictions from hydrologic models may not be improved by attempting to incorporate the understanding of these processes into hydrologic models. This paper quantifies the effect of simplifying watershed geometry and averaging the parameter values on simulations generated using the KINEROS2 model. Furthermore, it examines how these changes in model input effect model output. The model was applied on a small semiarid rangeland watershed in which runoff is generated by the infiltration excess mechanism. The study concludes that averaging input parameter values has little effect on runoff volume and peak in simulating runoff. However, geometric simplification does have an effect on runoff peak and volume, but it is not statistically significant. In contrast, both averaging input parameter values and geometric simplification have an effect on model‐predicted sediment yield. Copyright © 2005 John Wiley & Sons, Ltd.     

Continuous time random walks for non-local radial solute transport

This study formulates and analyzes continuous time random walk (CTRW) models in radial flow geometries for the quantification of non-local solute transport induced by heterogeneous flow distributions and by mobile–immobile mass transfer processes. To this end we derive a general CTRW framework in radial coordinates starting from the random walk equations for radial particle positions and times. The particle density, or solute concentration is governed by a non-local radial advection–dispersion equation (ADE). Unlike in CTRWs for uniform flow scenarios, particle transition times here depend on the radial particle position, which renders the CTRW non-stationary. As a consequence, the memory kernel characterizing the non-local ADE, is radially dependent. Based on this general formulation, we derive radial CTRW implementations that (i) emulate non-local radial transport due to heterogeneous advection, (ii) model multirate mass transfer (MRMT) between mobile and immobile continua, and (iii) quantify both heterogeneous advection in a mobile region and mass transfer between mobile and immobile regions. The expected solute breakthrough behavior is studied using numerical random walk particle tracking simulations. This behavior is analyzed by explicit analytical expressions for the asymptotic solute breakthrough curves. We observe clear power-law tails of the solute breakthrough for broad (power-law) distributions of particle transit times (heterogeneous advection) and particle trapping times (MRMT model). The combined model displays two distinct time regimes. An intermediate regime, in which the solute breakthrough is dominated by the particle transit times in the mobile zones, and a late time regime that is governed by the distribution of particle trapping times in immobile zones. These radial CTRW formulations allow for the identification of heterogeneous advection and mobile-immobile processes as drivers of anomalous transport, under conditions relevant for field tracer tests.     

Analytical solutions for the coefficient of variation of the volume-averaged solute concentration in heterogeneous aquifers

Under the assumption that local solute dispersion is negligible, a new general formula (in the form of a convolution integral) is found for the arbitrary k-point ensemble moment of the local concentration of a solute convected in arbitrary m spatial dimensions with general sure initial conditions. From this general formula new closed-form solutions in m=2 spatial dimensions are derived for 2-point ensemble moments of the local solute concentration for the impulse (Dirac delta) and Gaussian initial conditions. When integrated over an averaging window, these solutions lead to new closed-form expressions for the first two ensemble moments of thevolume-averaged solute concentration and to the corresponding concentration coefficients of variation (CV). Also, for the impulse (Dirac delta) solute concentration initial condition, the second ensemble moment of thesolute point concentration in two spatial dimensions and the corresponding CV are demonstrated to be unbound. For impulse initial conditions the CVs for volume-averaged concentrations axe compared with each other for a tracer from the Borden aquifer experiment. The point-concentration CV is unacceptably large in the whole domain, implying that the ensemble mean concentration is inappropriate for predicting the actual concentration values. The volume-averaged concentration CV decreases significantly with an increasing averaging volume. Since local dispersion is neglected, the new solutions should be interpreted as upper limits for the yet to be derived solutions that account for local dispersion; and so should the presented CVs for Borden tracers. The new analytical solutions may be used to test the accuracy of Monte Carlo simulations or other numerical algorithms that deal with the stochastic solute transport. They may also be used to determine the size of the averaging volume needed to make a quasi-sure statement about the solute mass contained in it.     

Modelling the tidal mixing fronts and seasonal stratification of the Northwest European Continental shelf

We investigate mixing processes under stratified conditions on the Northwest European Continental shelf using a numerical model (POLCOMS). Our results indicate that convection induced by vertical shearing of horizontal density gradients (‘shear-induced convection’) is a regularly occurring feature in the bottom and surface boundary layers in this open shelf-sea situation. Two types of turbulence models are investigated to study their capability for reproducing the observed location of tidal mixing fronts, and the physical processes occurring in seasonally stratified waters. The first model is a one-equation variant of the Mellor–Yamada model, whereas the second model combines a more recent second-momentum closure with a two-equation model. It is found that generally mean frontal positions (as estimated from ICES data) are predicted more accurately by the two-equation model. The one-equation model reproduces the mean frontal locations to 18.1 km (<3 grid spacings) and the two-equation model to 17.1 km; although in the Celtic Sea the accuracy is ∼33 and ∼12 km, respectively. Comparison with historical tide gauges, current metres, CTD stations, and thermistor chain data from the North Sea Project all show an improvement in accuracy when the two-equation model is used. This is particularly apparent in the model's ability to reproduce the spring–neap variability during stratification. We find that in the presence of shear-induced convection the routinely applied clipping of the turbulent length-scale, previously thought to be a minor ingredient in a turbulence model, has a dramatic effect on the results: if the length-scale clipping is not applied, substantial over-mixing is observed to occur. The causes and possible remedies of this effect are investigated. Overall our results demonstrate a sensitivity to the details of the turbulence model that is significantly greater than previously thought.     

Transport in chemically and mechanically heterogeneous porous media. I: Theoretical development of region-averaged equations for slightly compressible single-phase flow

Transport processes in heterogeneous porous media are often treated in terms of one-equation models. Such treatment assumes that the velocity, pressure, temperature, and concentration can be represented in terms of a single large-scale averaged quantity in regions having significantly different mechanical, thermal, and chemical properties. In this paper we explore the process of single-phase flow in a two-region model of heterogeneous porous media. The region-averaged equations are developed for the case of a slightly compressible flow which is an accurate representation for a certain class of liquid-phase flows. The analysis leads to a pair of transport equations for the region averaged pressures that are coupled through a classic exchange term, in addition to being coupled by a diffusive cross effect. The domain of validity of the theory has been identified in terms of a series of length and timescale constraints.In Part II the theory is tested, in the absence of adjustable parameters, by comparison with numerical experiments for transient, slightly compressible flow in both stratified and nodular models of heterogeneous porous media. Good agreement between theory and experiment is obtained for nodular and stratified systems, and effective transport coefficients for a wide range of conditions are presented on the basis of solutions of the three closure problems that appear in the theory. Part III of this paper deals with the principle of large-scale mechanical equilibrium and the region-averaged form of Darcy's law. This form is necessary for the development and solution of the region-averaged solute transport equations that are presented in Part IV. Finally, in Part V we present results for the dispersion tensors and the exchange coefficient associated with the two-region model of solute transport with adsorption.     

An eddifying Stommel model: fast eddy effects in a two-box ocean

ABSTRACT

A system of stochastic differential equations is formulated describing the heat and salt content of a two-box ocean. Variability in the heat and salt content and in the thermohaline circulation between the boxes is driven by fast Gaussian atmospheric forcing and by ocean-intrinsic, eddy-driven variability. The eddy forcing of the slow dynamics takes the form of a colored, non-Gaussian noise. The qualitative effects of this non-Gaussianity are investigated by comparing to two approximate models: one that includes only the mean eddy effects (the “averaged model”), and one that includes an additional Gaussian white-noise approximation of the eddy effects (the “Gaussian model”). Both of these approximate models are derived using the methods of fast averaging and homogenisation. In the parameter regime where the dynamics has a single stable equilibrium the averaged model has too little variability. The Gaussian model has accurate second-order statistics, but incorrect skew and rare-event probabilities. In the parameter regime where the dynamics has two stable equilibria the eddy noise is much smaller than the atmospheric noise. The averaged, Gaussian, and non-Gaussian models all have similar stationary distributions, but the jump rates between equilibria are too small for the averaged and Gaussian models.     

Exploring geochemical controls on weathering and erosion of convex hillslopes: beyond the empirical regolith production function

Landscape curvature evolves in response to physical, chemical, and biological influences that cannot yet be quantified in models. Nonetheless, the simplest models predict the existence of equilibrium hillslope profiles. Here, we develop a model describing steady‐state regolith production caused by mineral dissolution on hillslopes which have attained an equilibrium parabolic profile. When the hillslope lowers at a constant rate, the rate of chemical weathering is highest at the ridgetop where curvature is highest and the ridge develops the thickest regolith. This result derives from inclusion of all the terms in the mathematical definition of curvature. Including these terms shows that the curvature of a parabolic hillslope profile varies with distance from the ridge. The hillslope model (meter‐scale) is similar to models of weathering rind formation (centimeter‐scale) where curvature‐driven solute transport causes development of the thickest rinds at highly curved clast corners. At the clast scale, models fit observations. Here, we similarly explore model predictions of the effect of curvature at the hillslope scale. The hillslope model shows that when erosion rates are small and vertical porefluid infiltration is moderate, the hill weathers at both ridge and valley in the erosive transport‐limited regime. For this regime, the reacting mineral is weathered away before it reaches the land surface: in other words, the model predicts completely developed element‐depth profiles at both ridge and valley. In contrast, when the erosion rate increases or porefluid velocity decreases, denudation occurs in the weathering‐limited regime. In this regime, the reacting mineral does not weather away before it reaches the land surface and simulations predict incompletely developed profiles at both ridge and valley. These predictions are broadly consistent with observations of completely developed element‐depth profiles along hillslopes denuding under erosive transport‐limitation but incompletely developed profiles along hillslopes denuding under weathering limitation in some field settings. Copyright © 2013 John Wiley & Sons, Ltd.     

Modelling of solute transport in rivers under different flow rates: A case study without transient storage

A methodology to derive solute transport models at any flow rate is presented. The novelty of the proposed approach lies in the assessment of uncertainty of predictions that incorporate parameterisation based on flow rate. A simple treatment of uncertainty takes into account heteroscedastic modelling errors related to tracer experiments performed over a range of flow rates, as well as the uncertainty of the observed flow rates themselves. The proposed approach is illustrated using two models for the transport of a conservative solute: a physically based, deterministic, advection-dispersion model (ADE), and a stochastic, transfer function based, active mixing volume model (AMV). For both models the uncertainty of any parameter increases with increasing flow rate (reflecting the heteroscedastic treatment of modelling errors at different observed flow rates), but in contrast the uncertainty of travel time, computed from the predicted model parameters, was found to decrease with increasing flow rate.     

Non-equilibrium interphase heat and mass transfer during two-phase flow in porous media—Theoretical considerations and modeling

Mass and heat transfer occurring across phase-interfaces in multi-phase flow in porous media are mostly approximated using equilibrium relationships or empirical kinetic models. However, when the characteristic time of flow is smaller than that of mass or heat transfer, non-equilibrium situations may arise. Commonly, empirical approaches are used in such cases. There are only few works in the literature that use physically-based models for these transfer terms. In fact, one would expect physical approaches to modeling kinetic interphase mass and heat transfer to contain the interfacial area between the phases as a variable. Recently, a two-phase flow and solute transport model was developed that included interfacial area as a state variable [36]. In that model, interphase mass transfer was modeled as a kinetic process.     

A coupled model for simulating surface water and groundwater interactions in coastal wetlands

Coastal wetlands are characterized by strong, dynamic interactions between surface water and groundwater. This paper presents a coupled model that simulates interacting surface water and groundwater flow and solute transport processes in these wetlands. The coupled model is based on two existing (sub) models for surface water and groundwater, respectively: ELCIRC (a three‐dimensional (3‐D) finite‐volume/finite‐difference model for simulating shallow water flow and solute transport in rivers, estuaries and coastal seas) and SUTRA (a 3‐D finite‐element/finite‐difference model for simulating variably saturated, variable‐density fluid flow and solute transport in porous media). Both submodels, using compatible unstructured meshes, are coupled spatially at the common interface between the surface water and groundwater bodies. The surface water level and solute concentrations computed by the ELCIRC model are used to determine the boundary conditions of the SUTRA‐based groundwater model at the interface. In turn, the groundwater model provides water and solute fluxes as inputs for the continuity equations of surface water flow and solute transport to account for the mass exchange across the interface. Additionally, flux from the seepage face was routed instantaneously to the nearest surface water cell according to the local sediment surface slope. With an external coupling approach, these two submodels run in parallel using time steps of different sizes. The time step (Δt_g) for the groundwater model is set to be larger than that (Δt_s) used by the surface water model for computational efficiency: Δt_g = M × Δt_s where M is an integer greater than 1. Data exchange takes place between the two submodels through a common database at synchronized times (e.g. end of each Δt_g). The coupled model was validated against two previously reported experiments on surface water and groundwater interactions in coastal lagoons. The results suggest that the model represents well the interacting surface water and groundwater flow and solute transport processes in the lagoons. Copyright © 2011 John Wiley & Sons, Ltd.     

Assessment of alternative adsorption models and global sensitivity analysis to characterize hexavalent chromium loss from soil to surface runoff

We investigate our ability to assess transfer of hexavalent chromium, Cr(VI), from the soil to surface runoff by considering the effect of coupling diverse adsorption models with a two‐layer solute transfer model. Our analyses are grounded on a set of two experiments associated with soils characterized by diverse particle size distributions. Our study is motivated by the observation that Cr(VI) is receiving much attention for the assessment of environmental risks due to its high solubility, mobility, and toxicological significance. Adsorption of Cr(VI) is considered to be at equilibrium in the mixing layer under our experimental conditions. Four adsorption models, that is, the Langmuir, Freundlich, Temkin, and linear models, constitute our set of alternative (competing) mathematical formulations. Experimental results reveal that the soil samples characterized by the finest grain sizes are associated with the highest release of Cr(VI) to runoff. We compare the relative abilities of the four models to interpret experimental results through maximum likelihood model calibration and four model identification criteria (i.e., the Akaike information criteria [AIC and AIC_C] and the Bayesian and Kashyap information criteria). Our study results enable us to rank the tested models on the basis of a set of posterior weights assigned to each of them. A classical variance‐based global sensitivity analysis is then performed to assess the relative importance of the uncertain parameters associated with each of the models considered, within subregions of the parameter space. In this context, the modelling strategy resulting from coupling the Langmuir isotherm with a two‐layer solute transfer model is then evaluated as the most skilful for the overall interpretation of both sets of experiments. Our results document that (a) the depth of the mixing layer is the most influential factor for all models tested, with the exception of the Freundlich isotherm, and (b) the total sensitivity of the adsorption parameters varies in time, with a trend to increase as time progresses for all of the models. These results suggest that adsorption has a significant effect on the uncertainty associated with the release of Cr(VI) from the soil to the surface runoff component.     

Modeling preasymptotic transport in flows with significant inertial and trapping effects – The importance of velocity correlations and a spatial Markov model

We study solute transport in a periodic channel with a sinusoidal wavy boundary when inertial flow effects are sufficiently large to be important, but do not give rise to turbulence. This configuration and setup are known to result in large recirculation zones that can act as traps for solutes; these traps can significantly affect dispersion of the solute as it moves through the domain. Previous studies have considered the effect of inertia on asymptotic dispersion in such geometries. Here we develop an effective spatial Markov model that aims to describe transport all the way from preasymptotic to asymptotic times. In particular we demonstrate that correlation effects must be included in such an effective model when Péclet numbers are larger than O(100) in order to reliably predict observed breakthrough curves and the temporal evolution of second centered moments. For smaller Péclet numbers correlation effects, while present, are weak and do not appear to play a significant role. For many systems of practical interest, if Reynolds numbers are large, it may be typical that Péclet numbers are large also given that Schmidt numbers for typical fluids and solutes can vary between 1 and 500. This suggests that when Reynolds numbers are large, any effective theories of transport should incorporate correlation as part of the upscaling procedure, which many conventional approaches currently do not do. We define a novel parameter to quantify the importance of this correlation. Next, using the theory of CTRWs we explain a to date unexplained phenomenon of why dispersion coefficients for a fixed Péclet number increase with increasing Reynolds number, but saturate above a certain value. Finally we also demonstrate that effective preasymptotic models that do not adequately account for velocity correlations will also not predict asymptotic dispersion coefficients correctly.     

Lattice Boltzmann Models for Flow and Transport in Saturated Karst

Flow and transport simulation in karst aquifers remains a significant challenge for the ground water modeling community. Darcy's law–based models cannot simulate the inertial flows characteristic of many karst aquifers. Eddies in these flows can strongly affect solute transport. The simple two-region conduit/matrix paradigm is inadequate for many purposes because it considers only a capacitance rather than a physical domain. Relatively new lattice Boltzmann methods (LBMs) are capable of solving inertial flows and associated solute transport in geometrically complex domains involving karst conduits and heterogeneous matrix rock. LBMs for flow and transport in heterogeneous porous media, which are needed to make the models applicable to large-scale problems, are still under development. Here we explore aspects of these future LBMs, present simple examples illustrating some of the processes that can be simulated, and compare the results with available analytical solutions. Simulations are contrived to mimic simple capacitance-based two-region models involving conduit (mobile) and matrix (immobile) regions and are compared against the analytical solution. There is a high correlation between LBM simulations and the analytical solution for two different mobile region fractions. In more realistic conduit/matrix simulation, the breakthrough curve showed classic features and the two-region model fit slightly better than the advection-dispersion equation (ADE). An LBM-based anisotropic dispersion solver is applied to simulate breakthrough curves from a heterogeneous porous medium, which fit the ADE solution. Finally, breakthrough from a karst-like system consisting of a conduit with inertial regime flow in a heterogeneous aquifer is compared with the advection-dispersion and two-region analytical solutions.     

The impact of boundary on the fractional advection–dispersion equation for solute transport in soil: Defining the fractional dispersive flux with the Caputo derivatives

The inherent heterogeneity of geological media often results in anomalous dispersion for solute transport through them, and how to model it has been an interest over the past few decades. One promising approach that has been increasingly used to simulate the anomalous transport in surface and subsurface water is the fractional advection–dispersion equation (FADE), derived as a special case of the more general continuous time random walk or the stochastic continuum model. In FADE, the dispersion is not local and the solutes have appreciable probability to move long distances, and thus reach the boundary faster than predicted by the classical advection–dispersion equation (ADE). How to deal with different boundaries associated with FADE and their consequent impact is an issue that has not been thoroughly explored. In this paper we address this by taking one-dimensional solute movement in soil columns as an example. We show that the commonly used FADE with its fractional derivatives defined by the Riemann–Liouville definition is problematic and could result in unphysical results for solute transport in bounded domains; a modified method with the fractional dispersive flux defined by the Caputo derivatives is presented to overcome this problem. A finite volume approach is given to numerically solve the modified FADE and its associated boundaries. With the numerical model, we analyse the inlet-boundary treatment in displacement experiments in soil columns, and find that, as in ADE, treating the inlet as a prescribed concentration boundary gives rise to mass-balance errors and such errors could be more significant in FADE because of its non-local dispersion. We also discuss a less-documented but important issue in hydrology: how to treat the upstream boundary in analysing the lateral movement of tracer in an aquifer when the tracer is injected as a pulse. It is shown that the use of an infinite domain, as commonly assumed in literature, leads to unphysical backward dispersion, which has a significant impact on data interpretation. To avoid this, the upstream boundary should be flux-prescribed and located at the upstream edge of the injecting point. We apply the model to simulate the movement of Cl⁻ in a tracer experiment conducted in a saturated hillslope, and analyse in details the significance of upstream-boundary treatments in parameter estimation.