Thermodynamically Constrained Averaging Theory Approach for Modeling Flow and Transport Phenomena in Porous Medium Systems: 5. Single-Fluid-Phase Transport

This work is the fifth in a series of papers on the thermodynamically constrained averaging theory (TCAT) approach for modeling flow and transport phenomena in multiscale porous medium systems. The general TCAT framework and the mathematical foundation presented in previous works are used to develop models that describe species transport and single-fluid-phase flow through a porous medium system in varying physical regimes. Classical irreversible thermodynamics formulations for species in fluids, solids, and interfaces are developed. Two different approaches are presented, one that makes use of a momentum equation for each entity along with constitutive relations for species diffusion and dispersion, and a second approach that makes use of a momentum equation for each species in an entity. The alternative models are developed by relying upon different approaches to constrain an entropy inequality using mass, momentum, and energy conservation equations. The resultant constrained entropy inequality is simplified and used to guide the development of closed models. Specific instances of dilute and non-dilute systems are examined and compared to alternative formulation approaches.     

A displacement finite element modelling for convection-diffusion problems

The development of a displacement finite element formulation and its application to convective transport problems is presented. The formulation is based on the introduction of a generalized quantity defined as transport displacement. The governing equation is expressed in terms of this quantity and by using generalized coordinates a variational form of the governing equation is obtained. This equation may be solved by any numerical method, though it is of particular interest for application of the finite element method. Two finite element models are derived for the solution of convection-diffusion boundary value problems. The performance of the two element models is discussed and numerical results are given for different cases of convection and diffusion with two types of boundary conditions. The numerical results obtained show not only the efficiency of the numerical models in handling pure convection, pure diffusion and mixed convection-diffusion problems, but also good stability and accuracy. The applications of the developed numerical models are not limited to diffusion-convection problems but can also be applied to other types of problems such as mass transfer, hydrodynamics and wave propagation.     

Regional‐scale sedimentation process models from airborne gamma ray remote sensing and digital elevation data

Airborne gamma ray survey data were used to provide information on potassium, thorium and uranium concentrations in surface soil and rock in arid central Australia. Spatial patterns in these radioelements allow tracing of paths of sediment at catchment scale. Survey elevation data are combined with contour data to produce digital elevation models for terrain analysis, tracing of sediment flow paths and modelling of extreme floods. Gamma ray data show consistent variation with slope, a limited range of drainage areas, and erosion/deposition models derived from the conservation of mass equation. Supply‐limited sediment transport models give a reasonable reproduction of observed radioelement distribution but some elements of the distribution pattern reflect the area inundated by 500–1000 year floods rather than the effects of simple downslope movement. Partial area sediment supply models are derived by downstream accumulation of erosion and deposition rates calculated using the conservation of mass equation with transport laws based on slope alone and stream power. Comparison with observed radioelement patterns suggests that both transport laws apply in different parts of the landscape. Regional‐scale sediment transport models will require a range of models depending on location in the landscape and event frequency. This approach may allow estimation of sediment delivery ratios. Copyright © 2001 John Wiley & Sons, Ltd.     

Effective macroscopic transport parameters between parallel plates with constant concentration boundaries

A macroscopic transport model is developed, following the Taylor shear dispersion analysis procedure, for a 2D laminar shear flow between parallel plates possessing a constant specified concentration. This idealized geometry models flow with contaminant dissolution at pore-scale in a contaminant source zone and flow in a rock fracture with dissolving walls. We upscale a macroscopic transient transport model with effective transport coefficients of mean velocity, macroscopic dispersion, and first-order mass transfer rate. To validate the macroscopic model the mean concentration, covariance, and wall concentration gradient are compared to the results of numerical simulations of the advection–diffusion equation and the Graetz solution. Results indicate that in the presence of local-scale variations and constant concentration boundaries, the upscaled mean velocity and macrodispersion coefficient differ from those of the Taylor–Aris dispersion, and the mass transfer flux described by the first-order mass transfer model is larger than the diffusive mass flux from the constant wall. In addition, the upscaled first-order mass transfer coefficient in the macroscopic model depends only on the plate gap and diffusion coefficient. Therefore, the upscaled first-order mass transfer coefficient is independent of the mean velocity and travel distance, leading to a constant pore-scale Sherwood number of 12. By contrast, the effective Sherwood number determined by the diffusive mass flux is a function of the Peclet number for small Peclet number, and approaches a constant of 10.3 for large Peclet number.     

An inverse problem approach to modelling coastal effluent plumes

Formulated as an inverse problem, the diffusion parameters associated with length-scale dependent eddy diffusivities can be viewed as the unknowns in the mass conservation equation for coastal zone transport problems. The values of the diffusion parameters can be optimized according to an error function incorporated with observed concentration data. Examples are given for the Fickian, shear diffusion and inertial subrange diffusion models. Based on a new set of dyeplume data collected in the coastal zone off Bronte, Lake Ontario, it is shown that the predictions of turbulence closure models can be evaluated for different flow conditions. The choice of computational schemes for this diagnostic approach is based on tests with analytic solutions and observed data. It is found that the optimized shear diffusion model produced a better agreement with observations for both high and low advective flows than, e.g., the unoptimized semi-empirical model, K_y=0.075 σ_y^1.2, described by Murthy and Kenney.     

Impact of historic land‐use change on sediment delivery to a Chesapeake Bay subestuarine delta

Historic land use in the Chesapeake Bay drainage basin induced large fluxes of fluvial sediment to subestuarine tributaries. Stratigraphic and palaeoecologic analyses of deltaic deposits may be used to infer changes on the landscape, but are not sufficient to quantify past sediment supply. When viewed as an inverse boundary‐value problem, reconstruction of the sediment supply function may be achieved by combining deltaic sedimentation chronologies with an equation governing delta progradation. We propose that the diffusion equation is appropriate for simulating delta progradation and obtaining the sediment supply function provided a suitable diffusion constant (D) can be determined. Three new methods for estimating D are presented for the case of estuarine deltas. When the inverse boundary‐value technique was applied to Otter Point Creek, a tidal freshwater delta at the head of Bush River in upper Chesapeake Bay, D values ranged from 3763 to 6199 m² a^?1. Delta growth simulations showed a 1740–1760 initial pulse, a 1760–1780 erosive/redistributive interval, a 1780–1920 growth period, and a 1920‐present erosive/redistributive era. Coupling of simulated delta elevations with an empirical plant habitat predictive equation allowed for comparison of predicted versus actual relative habitat areas. Also, the model yielded reconstructed watershed erosion rates and stream suspended sediment concentrations that could be useful for development of water quality regulations. Copyright © 2001 John Wiley & Sons, Ltd.     

Modelling soil solute release and transport in run‐off on a loessial slope with and without surface stones

Stone covers on loessial slopes can increase the time of infiltration by slowing the velocity of the overland flow, which reduces the transport of solutes, but few mechanistic models have been tested under water‐scouring conditions. We carried out field experiments to test a previously proposed, physically based model of water and solute transport. The area of soil infiltration was calculated from the uncovered surface area, and Richards' equation and the kinematic wave equation were used to describe water infiltration and flow along slopes with stone covers. The transport of chemicals into the run‐off from the surface soil, presumably by diffusion, and their movement in the soil profile could be described by the convection–diffusion equations of the model. The simulated and measured data correlated well. The stones on the soil surface reduced the area available for infiltration but increased the Manning coefficient, eventually leading to increased water infiltration and decreased solute loss with run‐off. Our results indicated that the traditional model of water movement and solute migration could be used to simulate water transport and solute migration for stone‐covered soil on loessial slopes.     

A three-dimensional Eulerian model for transport and deposition of volcanic ashes

Volcanic ash fallout represents a serious threat to people living near active volcanoes because it can produce several undesirable effects such as collapse of roofs by ash loading, respiratory sickness, air traffic disruption, or damage to agriculture. The assessment of such volcanic risk is therefore an issue of vital importance for public safety and its mitigation often requires to evaluate the temporal evolution of the phenomenon through reliable computational models.We develop an Eulerian model, named FALL3D, for the transport and deposition of volcanic ashes. The model is based on the advection–diffusion–sedimentation equation with a turbulent diffusion given by the gradient transport theory, a wind field obtained from a meteorological limited area model (LAM) and the source term derived from by buoyant plume theory. It can be used to forecast either ash concentration in the atmosphere or ash loading on the ground. Model inputs are topography, meteorological data given by a LAM, mass eruption rate, and a particle settling velocity distribution. A test application to the July 2001 Etna eruption is presented.     

A Boltzmann-based mesoscopic model for contaminant transport in flow systems

The objective of this paper is to demonstrate the formulation of a numerical model for mass transport based on the Bhatnagar–Gross–Krook (BGK) Boltzmann equation. To this end, the classical chemical transport equation is derived as the zeroth moment of the BGK Boltzmann differential equation. The relationship between the mass transport equation and the BGK Boltzmann equation allows an alternative approach to numerical modeling of mass transport, wherein mass fluxes are formulated indirectly from the zeroth moment of a difference model for the BGK Boltzmann equation rather than directly from the transport equation. In particular, a second-order numerical solution for the transport equation based on the discrete BGK Boltzmann equation is developed. The numerical discretization of the first-order BGK Boltzmann differential equation is straightforward and leads to diffusion effects being accounted for algebraically rather than through a second-order Fickian term. The resultant model satisfies the entropy condition, thus preventing the emergence of non-physically realizable solutions including oscillations in the vicinity of the front. Integration of the BGK Boltzmann difference equation into the particle velocity space provides the mass fluxes from the control volume and thus the difference equation for mass concentration. The difference model is a local approximation and thus may be easily included in a parallel model or in accounting for complex geometry. Numerical tests for a range of advection–diffusion transport problems, including one- and two-dimensional pure advection transport and advection–diffusion transport show the accuracy of the proposed model in comparison to analytical solutions and solutions obtained by other schemes.     

Magnetic resonance microscopy of biofouling induced scale dependent transport in porous media

Non-invasive magnetic resonance microscopy (MRM) methods are applied to study biofouling of a homogeneous model porous media. MRM of the biofilm biomass using magnetic relaxation weighting shows the heterogeneous nature of the spatial distribution of the biomass as a function of growth. Spatially resolved MRM velocity maps indicate a strong variation in the pore scale velocity as a function of biofilm growth. The hydrodynamic dispersion dynamics for flow through the porous media is quantitatively characterized using a pulsed gradient spin echo technique to measure the propagator of the motion. The propagator indicates a transition in transport dynamics from a Gaussian normal diffusion process following a normal advection diffusion equation to anomalous transport as a function of biofilm growth. Continuous time random walk models resulting in a time fractional advection diffusion equation are shown to model the transition from normal to anomalous transport in the context of a conceptual model for the biofouling. The initially homogeneous porous media is transformed into a more complex heterogeneous porous media by the biofilm growth.     

Gopher bioturbation: field evidence for non‐linear hillslope diffusion

It has generally been assumed that diffusive sediment transport on soil‐mantled hillslopes is linearly dependent on hillslope gradient. Fieldwork was done near Santa Barbara, California, to develop a sediment transport equation for bioturbation by the pocket gopher (Thomomys bottae) and to determine whether it supports linear diffusion. The route taken by the sediment is divided into two parts, a subsurface path followed by a surface path. The first is the transport of soil through the burrow to the burrow opening. The second is the discharge of sediment from the burrow opening onto the hillslope surface. The total volumetric sediment flux, as a function of hillslope gradient, is found to be: q_s (cm³ cm⁻¹ a⁻¹) = 176(dz/dx)³ − 189(dz/dx)² + 68(dz/dx) + 34(dz/dx)^0·4. This result does not support the use of linear diffusion for hillslopes where gopher bioturbation is the dominant mode of sediment transport. A one‐dimensional hillslope evolution program was used to evolve hillslope profiles according to non‐linear and linear diffusion and to compare them to a typical hillslope. The non‐linear case more closely resembles the actual profile with a convex cap at the divide leading into a straight midslope section. Copyright © 2000 John Wiley & Sons, Ltd.     

Groundwater transport modeling with nonlinear sorption and intraparticle diffusion

Despite recent advances in the mechanistic understanding of sorption in groundwater systems, most contaminant transport models provide limited support for nonideal sorption processes such as nonlinear isotherms and/or diffusion-limited sorption. However, recent developments in the conceptualization of “dual mode” sorption for hydrophobic organic contaminants have provided more realistic and mechanistically sound alternatives to the commonly used Langmuir and Freundlich models. To support the inclusion of both nonlinear and diffusion-limited sorption processes in groundwater transport models, this paper presents two numerical algorithms based on the split operator approach. For the nonlinear equilibrium scenario, the commonly used two-step split operator algorithm has been modified to provide a more robust treatment of complex multi-parameter isotherms such as the Polanyi-partitioning model. For diffusion-limited sorption, a flexible three step split-operator procedure is presented to simulate intraparticle diffusion in multiple spherical particles with different sizes and nonlinear isotherms. Numerical experiments confirmed the accuracy of both algorithms for several candidate isotherms. However, the primary advantages of the algorithms are: (1) flexibility to accommodate any isotherm equation including “dual mode” and similar expressions, and (2) ease of adapting existing grid-based transport models of any dimensionality to include nonlinear sorption and/or intraparticle diffusion. Comparisons are developed for one-dimensional transport scenarios with different isotherms and particle configurations. Illustrative results highlight (1) the potential influence of isotherm model selection on solute transport predictions, and (2) the combined effects of intraparticle diffusion and nonlinear sorption on the plume transport and flushing for both single-particle and multi-particle scenarios.     

Importance of advective zone in longitudinal mixing experiments

One-dimensional Fickian dispersion models such as the advection diffusion equation (ADE) are commonly used to analyse and predict concentration distributions downstream of contamination events in watercourses. Such models are only valid once the tracer had entered the equilibrium zone. This paper compares previous theoretical, experimental and numerical estimates of the distance to reach the equilibrium zone with new experimental values, obtained by examining the change of skewness in a tracer profile, downstream of a cross-sectionally well mixed source. Closer agreement was found with Fischers’ theoretical estimate than prior experimental and numerical studies.     

First-principles derivation of reactive transport modeling parameters for particle tracking and PDE approaches

Both Eulerian and Lagrangian reactive transport simulations in natural media require selection of a parameter that controls the “promiscuity” of the reacting particles. In Eulerian models, measurement of this parameter may be difficult because its value will generally differ between natural (diffusion-limited) systems and batch experiments, even though both are modeled by reaction terms of the same form. And in Lagrangian models, there previously has been no a priori way to compute this parameter. In both cases, then, selection is typically done by calibration, or ad hoc. This paper addresses the parameter selection problem for Fickian transport by deriving, from first principles and D (the diffusion constant) the reaction-rate-controlling parameters for particle tracking (PT) codes and for the diffusion–reaction equation (DRE). Using continuous time random walk analysis, exact reaction probabilities are derived for pairs of potentially reactive particles based on D and their probability of reaction provided that they collocate. Simultaneously, a second PT scheme directly employing collocation probabilities is derived. One-to-one correspondence between each of D, the reaction radius specified for a PT scheme, and the DRE decay constant are then developed. These results serve to ground reactive transport simulations in their underlying thermodynamics, and are confirmed by simulations.     

Effective hydraulic conductivity of nonstationary aquifers

Assuming that the ln hydraulic conductivity in an aquifer is mathematically approximated by a spatial deterministic surface, or trend, plus a stationary random noise, we treat the problem of finding what the effective hydraulic conductivity of that aquifer is. This problem is tackled by spectral methods applied to a type of diffusion equation of groundwater flow, together with suitable coordinate transformations. Analytical (exact) solutions in terms of elementary functions are presented for one- and three-dimensional finite and infinite domains. Stability criteria are obtained for the solutions, in terms of a critical parameter, that turns out to involve the product of correlation scale and trend gradient. For the case of finite and symmetrical domains, additional provisions to insure the stability of numerical calculations of effective hydraulic conductivity are provided. Effective hydraulic conductivity is an important property, with potential applications in the calibrations of groundwater and transport numerical models.     

Effective hydraulic conductivity of nonstationary aquifers

Assuming that the ln hydraulic conductivity in an aquifer is mathematically approximated by a spatial deterministic surface, or trend, plus a stationary random noise, we treat the problem of finding what the effective hydraulic conductivity of that aquifer is. This problem is tackled by spectral methods applied to a type of diffusion equation of groundwater flow, together with suitable coordinate transformations. Analytical (exact) solutions in terms of elementary functions are presented for one- and three-dimensional finite and infinite domains. Stability criteria are obtained for the solutions, in terms of a critical parameter, that turns out to involve the product of correlation scale and trend gradient. For the case of finite and symmetrical domains, additional provisions to insure the stability of numerical calculations of effective hydraulic conductivity are provided. Effective hydraulic conductivity is an important property, with potential applications in the calibrations of groundwater and transport numerical models.     

Non-Fickian mass transport in fractured porous media

The paper provides an introduction to fundamental concepts of mathematical modeling of mass transport in fractured porous heterogeneous rocks. Keeping aside many important factors that can affect mass transport in subsurface, our main concern is the multi-scale character of the rock formation, which is constituted by porous domains dissected by the network of fractures. Taking into account the well-documented fact that porous rocks can be considered as a fractal medium and assuming that sizes of pores vary significantly (i.e. have different characteristic scales), the fractional-order differential equations that model the anomalous diffusive mass transport in such type of domains are derived and justified analytically. Analytical solutions of some particular problems of anomalous diffusion in the fractal media of various geometries are obtained. Extending this approach to more complex situation when diffusion is accompanied by advection, solute transport in a fractured porous medium is modeled by the advection-dispersion equation with fractional time derivative. In the case of confined fractured porous aquifer, accounting for anomalous non-Fickian diffusion in the surrounding rock mass, the adopted approach leads to introduction of an additional fractional time derivative in the equation for solute transport. The closed-form solutions for concentrations in the aquifer and surrounding rocks are obtained for the arbitrary time-dependent source of contamination located in the inlet of the aquifer. Based on these solutions, different regimes of contamination of the aquifers with different physical properties can be readily modeled and analyzed.     

Ammonium Transport Simulation in Horizontal Soil Columns from a Natural Inland Alkaline Wetland

Ammonium transport was simulated in horizontal soil columns from an inland alkaline wetland (Fulaowenpao wetland) of Northeast China. The primary objectives of this work are to investigate the changes in ammonium transport rate with increasing distances along horizontal soil column and to determine the effects of water diffusion rate and volumetric water content on ammonium transport rate. Our results showed that water diffusion coefficient was the lowest at the soil layer of 10–20 cm, followed by the 0–10 cm soil layer, and the highest value occured at the soil layer of 20–60 cm. The highest ammonium transport rate also appeared at the soil layer of 20–60 cm, while the lowest value was observed at the soil layer of 10–20 cm. Ammonium transport rates decreased with increasing distances along horizontal soil columns. The ammonium transport rates showed higher values at the distance from 0 to 6 cm and then decreased rapidly from 6 to 18 cm. However, they nearly kept constant and approached to zero after exceeding the distance of 18 cm. Ammonium transport rates increased exponentially with increasing volumetric water contents and water diffusion rates. The alkaline wetland soils prevented ammonium from horizontal diffusion at all soil layers under drying conditions.     

Two‐dimensional nonlinear diffusive numerical simulation of geomorphic modifications to cinder cones

The temporal evolution of simple landforms such as cinder cones by nonlinear diffusive processes is studied through the use of a new 2D numerical model using well‐established and accurate numerical mathematics and high‐resolution digital elevation models (DEMs). Extending 1D (profile) nonlinear diffusion analyses used in cinder cone, hillslope and fault scarp evolution studies, we have implemented a 2D numerical model with a spatially and temporally varying sediment transport rate coefficient scaled nonlinearly by the ratio of local slope to critical slope. The high accuracy and efficient numerical implementation are documented in the paper and the MATLAB toolkit developed is used to solve for the developmentof an initial 2D cone form. First, we examine the nonlinear transport rule and suggest a refinement that accounts explicitly for flux at threshold slopes. We find that the maximum diffusion (necessarily introduced in the numerical model to avoid infinite rates) at the critical slope controls the final morphology, especially approaching steady state. Secondly, solving the landscape evolution problem in 2D enables a natural accounting for sediment flux convergence or divergence in the profile. Thirdly, the boundary behavior of a given landscape element controls much of what happens in that domain and so we allow for arbitrary flux magnitude or elevation boundary conditions. Fourthly, landscapes are heterogeneous in their surface cover and so we allow for spatially and temporally varying transport rate k and we permit an arbitrary vertical displacement field within the model domain. To test the new formulation for the nonlinear term, the effect of variable diffusivity k and the numerical schemes implemented, we apply the model to cinder cones built on the flanks of Mount Etna in 2001 and 2002–2003. We explore the effects of DEM resolution with data from the 2001 cone and the utility of spatially variable diffusivity to explain the variation in erosion measured by differencing repeat light detection and ranging (LIDAR) surveys gathered in 2004 and 2007 over the 2002–2003 cone complex. Copyright © 2013 John Wiley & Sons, Ltd.