Time–space fractional governing equations of transient groundwater flow in confined aquifers: Numerical investigation

In order to model non‐Fickian transport behaviour in groundwater aquifers, various forms of the time–space fractional advection–dispersion equation have been developed and used by several researchers in the last decade. The solute transport in groundwater aquifers in fractional time–space takes place by means of an underlying groundwater flow field. However, the governing equations for such groundwater flow in fractional time–space are yet to be developed in a comprehensive framework. In this study, a finite difference numerical scheme based on Caputo fractional derivative is proposed to investigate the properties of a newly developed time–space fractional governing equations of transient groundwater flow in confined aquifers in terms of the time–space fractional mass conservation equation and the time–space fractional water flux equation. Here, we apply these time–space fractional governing equations numerically to transient groundwater flow in a confined aquifer for different boundary conditions to explore their behaviour in modelling groundwater flow in fractional time–space. The numerical results demonstrate that the proposed time–space fractional governing equation for groundwater flow in confined aquifers may provide a new perspective on modelling groundwater flow and on interpreting the dynamics of groundwater level fluctuations. Additionally, the numerical results may imply that the newly derived fractional groundwater governing equation may help explain the observed heavy‐tailed solute transport behaviour in groundwater flow by incorporating nonlocal or long‐range dependence of the underlying groundwater flow field.     

Adaptive volume penalization for ocean modeling

The development of various volume penalization techniques for use in modeling topographical features in the ocean is the focus of this paper. Due to the complicated geometry inherent in ocean boundaries, the stair-step representation used in the majority of current global ocean circulation models causes accuracy and numerical stability problems. Brinkman penalization is the basis for the methods developed here and is a numerical technique used to enforce no-slip boundary conditions through the addition of a term to the governing equations. The second aspect to this proposed approach is that all governing equations are solved on a nonuniform, adaptive grid through the use of the adaptive wavelet collocation method. This method solves the governing equations on temporally and spatially varying meshes, which allows higher effective resolution to be obtained with less computational cost. When penalization methods are coupled with the adaptive wavelet collocation method, the flow near the boundary can be well-resolved. It is especially useful for simulations of boundary currents and tsunamis, where flow near the boundary is important. This paper will give a thorough analysis of these methods applied to the shallow water equations, as well as some preliminary work applying these methods to volume penalization for bathymetry representation for use in either the nonhydrostatic or hydrostatic primitive equations.     

Analytical approximations for flow in compressible,saturated, one-dimensional porous media

A nonlinear model for single-phase fluid flow in slightly compressible porous media is presented and solved approximately. The model assumes state equations for density, porosity, viscosity and permeability that are exponential functions of the fluid (either gas or liquid) pressure. The governing equation is transformed into a nonlinear diffusion equation. It is solved for a semi-infinite domain for either constant pressure or constant flux boundary conditions at the surface. The solutions obtained, although approximate, are extremely accurate as demonstrated by comparisons with numerical results. Predictions for the surface pressure resulting from a constant flux into a porous medium are compared with published experimental data.     

A generalized transformation approach for simulating steady-state variably-saturated subsurface flow

A numerical method is proposed to accurately and efficiently compute a direct steady-state solution of the nonlinear Richards equation. In the proposed method, the Kirchhoff integral transformation and a complementary transformation are applied to the governing equation in order to separate the nonlinear hyperbolic characteristic from the linear parabolic part. The separation allows the transformed governing equation to be applied to partially- to fully-saturated systems with arbitrary constitutive relations between primary (pressure head) and secondary variables (relative permeability). The transformed governing equation is then discretized with control volume finite difference/finite element approximations, followed by inverse transformation. The approach is compared to analytical and other numerical approaches for variably-saturated flow in 1-D and 3-D domains. The results clearly demonstrate that the approach is not only more computationally efficient but also more accurate than traditional numerical solutions. The approach is also applied to an example flow problem involving a regional-scale variably-saturated heterogeneous system, where the vadose zone is up to 1 km thick. The performance, stability, and effectiveness of the transform approach is exemplified for this complex heterogeneous example, which is typical of many problems encountered in the field. It is shown that computational performance can be enhanced by several orders of magnitude with the described integral transformation approach.     

Time–space fractional governing equations of one‐dimensional unsteady open channel flow process: Numerical solution and exploration

Although fractional integration and differentiation have found many applications in various fields of science, such as physics, finance, bioengineering, continuum mechanics, and hydrology, their engineering applications, especially in the field of fluid flow processes, are rather limited. In this study, a finite difference numerical approach is proposed to solve the time–space fractional governing equations of 1‐dimensional unsteady/non‐uniform open channel flow process. By numerical simulations, results of the proposed fractional governing equations of the open channel flow process were compared with those of the standard Saint‐Venant equations. Numerical simulations showed that flow discharge and water depth can exhibit heavier tails in downstream locations as space and time fractional derivative powers decrease from 1. The fractional governing equations under consideration are generalizations of the well‐known Saint‐Venant equations, which are written in the integer differentiation framework. The new governing equations in the fractional‐order differentiation framework have the capability of modelling nonlocal flow processes both in time and in space by taking the global correlations into consideration. Furthermore, the generalized flow process may possibly shed light on understanding the theory of the anomalous transport processes and observed heavy‐tailed distributions of particle displacements in transport processes.     

New analytical solutions for groundwater flow in wedge-shaped aquifers with various topographic boundary conditions

The hydraulic head distribution in a wedge-shaped aquifer depends on the wedge angle and the topographic and hydrogeological boundary conditions. In addition, an equation in terms of the radial distance with trigonometric functions along the boundary may be suitable to describe the water level configuration for a valley flank with a gentle sloping and rolling topography. This paper develops a general mathematical model including the governing equation and a variety of boundary conditions for the groundwater flow within a wedge-shaped aquifer. Based on the model, a new closed-form solution for transient flow in the wedge-shaped aquifer is derived via the finite sine transform and Hankel transform. In addition, a numerical approach, including the roots search scheme, the Gaussian quadrature, and Shanks’ method, is proposed for efficiently evaluating the infinite series and the infinite integral presented in the solution. This solution may be used to describe the head distribution for wedges that image theory is inapplicable, and to explore the effects of the recharge from various topographic boundaries on the groundwater flow system within a wedge-shaped aquifer.     

Finite element analysis of three-dimensional groundwater flow problems

A numerical method is presented for analysing either steady state or transient three-dimensional groundwater flow problems. The governing equation is formulated in terms of the finite element process using the Galerkin approach, and cubic isoparametric elements are used to simulate the flow domain as these permit accurate modelling of curved boundaries. Particular attention is paid to the time dependent movement of the phreatic surface where an iterative technique based on the replacement of the original transient problem by a discrete number of steady state problems is used to effect a solution. Furthermore, in tracing the movement of the surface use is made of the element formulation theory in order to compute the normal to the boundary.The validity of the technique is first established by analysing a radially symmetrical problem for which an alternative analytical solution is available. Finally, a general three-dimensional flow system is studied for which there is no known analytical solution. It is shown that relatively few elements are required to yield practical solutions.     

Numerical solution of two-phase flow for the advection-dominated and non-linear case

This work presents a highly efficient numerical scheme for solving immiscible, advection-dominated two-phase flow in heterogeneous porous media. The pressure equation is decoupled from the saturation equation using an IMPES approach, while the advective terms are decoupled from the capillary diffusive terms in the saturation equation through sequential operator splitting. The parabolic and hyperbolic equations are approximated in time by implicit and explicit schemes, respectively. Damped Newton linearization is applied to the implicit non-linear diffusive step. Mixed hybrid finite elements are applied to the global pressure equation and to the regularized capillary diffusion term. For both linear systems arising from the approximation procedure, an AMG preconditioned conjugate gradient solver is used. A finite volume scheme with slope limiter is applied to the advective step. Numerical comparison with standard preconditioners demonstrates the reliability of the proposed AMG-preconditioner. Benchmark examples illustrate the robustness of the method.     

Fractional calculus in hydrologic modeling: A numerical perspective

Fractional derivatives can be viewed either as handy extensions of classical calculus or, more deeply, as mathematical operators defined by natural phenomena. This follows the view that the diffusion equation is defined as the governing equation of a Brownian motion. In this paper, we emphasize that fractional derivatives come from the governing equations of stable Lévy motion, and that fractional integration is the corresponding inverse operator. Fractional integration, and its multi-dimensional extensions derived in this way, are intimately tied to fractional Brownian (and Lévy) motions and noises. By following these general principles, we discuss the Eulerian and Lagrangian numerical solutions to fractional partial differential equations, and Eulerian methods for stochastic integrals. These numerical approximations illuminate the essential nature of the fractional calculus.     

Time-of-flight (TOF)-based two-phase upscaling for subsurface flow and transport

Subsurface formations are characterized by heterogeneity over multiple length scales, which can have a strong impact on flow and transport. In this paper, we present a new upscaling approach, based on time-of-flight (TOF), to generate upscaled two-phase flow functions. The method focuses on more accurate representations of local saturation boundary conditions, which are found to have a dominant impact (in comparison to the pressure boundary conditions) on the upscaled two-phase flow models. The TOF-based upscaling approach effectively incorporates single-phase flow and transport information into local upscaling calculations, accounting for the global flow effects on saturation, as well as the local variations due to subgrid heterogeneity. The method can be categorized into quasi-global upscaling techniques, as the global single-phase flow and transport information is incorporated in the local boundary conditions. The TOF-based two-phase upscaling can be readily integrated into any existing local two-phase upscaling framework, thus more flexible than local–global two-phase upscaling approaches developed recently. The method was applied to permeability fields with different correlation lengths and various fluid-mobility ratios. It was shown that the new method consistently outperforms existing local two-phase upscaling techniques, including recently developed methods with improved local boundary conditions (such as effective flux boundary conditions), and provides accurate coarse-scale models for both flow and transport.     

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.     

On the application of the boundary layer approximation for the simulation of density stratified flows in aquifers

Stratification of the density in groundwater flow stems from the contact between water which contains minerals in low concentration with water containing a high concentration of minerals. The flow in such a flow field should be simulated by solving simultaneously the equations of continuity, motion and solute transport, because solute concentration affects the dynamics of the flow. Such an approach is generally associated with complicated calculations and numerical schemes subject to problems of convergence and stability, as the basic equations are highly nonlinear.This study applies the phenomenological boundary layer approximation, and suggests a reference to three different zones in the flow field: (a) fresh water zone, (b) transition zone, and (c) mineralized water zone. In zones (a) and (c) it is assumed that the potential flow theory can be applied. In zone (b) the flow is nonpotential but the basic similarity conditions typical to boundary layers exist.The approach suggested in this study simplifies the mathematical models that should be used for the flow field simulation. This approach is especially attractive in cases where the Dupuit approximation is applicable. In such cases very often analytical solutions can be obtained for unidirectional flows. In cases that are too complicated for representation by analytical solutions, the method can be used for the creation of simplified numerical schemes.Various examples in this study demonstrate the application of the method for various field problems associated with steady state as well as unsteady state conditions.The simplicity of the method makes it useful for variety of problems. It can be used even by small institutions and small consulting firms, who have usually access to minicomputers and microprocessors.     

Perfectly matched layer boundary conditions for the second-order acoustic wave equation solved by the rapid expansion method

We derive a governing second-order acoustic wave equation in the time domain with a perfectly matched layer absorbing boundary condition for general inhomogeneous media. Besides, a new scheme to solve the perfectly matched layer equation for absorbing reflections from the model boundaries based on the rapid expansion method is proposed. The suggested scheme can be easily applied to a wide class of wave equations and numerical methods for seismic modelling. The absorbing boundary condition method is formulated based on the split perfectly matched layer method and we employ the rapid expansion method to solve the derived new perfectly matched layer equation. The use of the rapid expansion method allows us to extrapolate wavefields with a time step larger than the ones commonly used by traditional finite-difference schemes in a stable way and free of dispersion noise. Furthermore, in order to demonstrate the efficiency and applicability of the proposed perfectly matched layer scheme, numerical modelling examples are also presented. The numerical results obtained with the put forward perfectly matched layer scheme are compared with results from traditional attenuation absorbing boundary conditions and enlarged models as well. The analysis of the numerical results indicates that the proposed perfectly matched layer scheme is significantly effective and more efficient in absorbing spurious reflections from the model boundaries.     

Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures

Contrast in capillary pressure of heterogeneous permeable media can have a significant effect on the flow path in two-phase immiscible flow. Very little work has appeared on the subject of capillary heterogeneity despite the fact that in certain cases it may be as important as permeability heterogeneity. The discontinuity in saturation as a result of capillary continuity, and in some cases capillary discontinuity may arise from contrast in capillary pressure functions in heterogeneous permeable media leading to complications in numerical modeling. There are also other challenges for accurate numerical modeling due to distorted unstructured grids because of the grid orientation and numerical dispersion effects. Limited attempts have been made in the literature to assess the accuracy of fluid flow modeling in heterogeneous permeable media with capillarity heterogeneity. The basic mixed finite element (MFE) framework is a superior method for accurate flux calculation in heterogeneous media in comparison to the conventional finite difference and finite volume approaches. However, a deficiency in the MFE from the direct use of fractional flow formulation has been recognized lately in application to flow in permeable media with capillary heterogeneity. In this work, we propose a new consistent formulation in 3D in which the total velocity is expressed in terms of the wetting-phase potential gradient and the capillary potential gradient. In our formulation, the coefficient of the wetting potential gradient is in terms of the total mobility which is smoother than the wetting mobility. We combine the MFE and discontinuous Galerkin (DG) methods to solve the pressure equation and the saturation equation, respectively. Our numerical model is verified with 1D analytical solutions in homogeneous and heterogeneous media. We also present 2D examples to demonstrate the significance of capillary heterogeneity in flow, and a 3D example to demonstrate the negligible effect of distorted meshes on the numerical solution in our proposed algorithm.     

An efficient numerical model for incompressible two-phase flow in fractured media

Various numerical methods have been used in the literature to simulate single and multiphase flow in fractured media. A promising approach is the use of the discrete-fracture model where the fracture entities in the permeable media are described explicitly in the computational grid. In this work, we present a critical review of the main conventional methods for multiphase flow in fractured media including the finite difference (FD), finite volume (FV), and finite element (FE) methods, that are coupled with the discrete-fracture model. All the conventional methods have inherent limitations in accuracy and applications. The FD method, for example, is restricted to horizontal and vertical fractures. The accuracy of the vertex-centered FV method depends on the size of the matrix gridcells next to the fractures; for an acceptable accuracy the matrix gridcells next to the fractures should be small. The FE method cannot describe properly the saturation discontinuity at the matrix–fracture interface. In this work, we introduce a new approach that is free from the limitations of the conventional methods. Our proposed approach is applicable in 2D and 3D unstructured griddings with low mesh orientation effect; it captures the saturation discontinuity from the contrast in capillary pressure between the rock matrix and fractures. The matrix–fracture and fracture–fracture fluxes are calculated based on powerful features of the mixed finite element (MFE) method which provides, in addition to the gridcell pressures, the pressures at the gridcell interfaces and can readily model the pressure discontinuities at impermeable faults in a simple way. To reduce the numerical dispersion, we use the discontinuous Galerkin (DG) method to approximate the saturation equation. We take advantage of a hybrid time scheme to alleviate the restrictions on the size of the time step in the fracture network. Several numerical examples in 2D and 3D demonstrate the robustness of the proposed model. Results show the significance of capillary pressure and orders of magnitude increase in computational speed compared to previous works.     

Closed-form solutions for unsaturated flow under variable flux boundary conditions

It is shown that from any solution of the linear diffusion equation, we may construct a solution of a realistic form of the Richards equation for unsaturated flow. Compared to the usual direct linearization method, our inverse approach involves a quite different sequence of transformations. This opens the possibility of exact solutions with a wider variety of continuously varying flux boundary conditions. Closed-form solutions are presented for two examples. In these, the varying water flux boundary conditions resemble (i) the passage of a peaking storm and (ii) the continuous opening of a valve preceding a steady water supply. Unlike earlier more systematic approaches to this problem, our method does not require the numerical solution of an integral equation.     

Analysis of hydrodynamic conditions in adjacent free and heterogeneous porous flow domains

The existence of a free‐flow domain (e.g. a liquid layer) adjacent to a porous medium is a common occurrence in many environmental and petroleum engineering problems. The porous media may often contain various forms of heterogeneity, e.g. layers, fractures, micro‐scale lenses, etc. These heterogeneities affect the pressure distribution within the porous domain. This may influence the hydrodynamic conditions at the free–porous domain interface and, hence, the combined flow behaviour. Under steady‐state conditions, the heterogeneities are known to have negligible effects on the coupled flow behaviour. However, the significance of the heterogeneity effects on coupled free and porous flow under transient conditions is not certain. In this study, numerical simulations have been carried out to investigate the effects of heterogeneous (layered) porous media on the hydrodynamics conditions in determining the behaviour of combined free and porous regimes. Heterogeneity in the porous media is introduced by defining a domain composed of two layers of porous media with different values of intrinsic permeability. The coupling of the governing equations of motion in free and porous domains has been achieved through the well‐known Beavers and Joseph interfacial condition. Of special interest in this work are porous domains with flow‐through ends. They represent the general class of problems where large physical domains are truncated to smaller sections for ease of mathematical analysis. However, this causes a practical difficulty in modelling such systems. This is because the information on flow behaviour, i.e. boundary conditions at the truncated sections, is usually not available. Use of artificial boundary conditions to solve these problems effectively implies the imposition of conditions that do not necessarily match with the solutions required for the interior of the domain. This difficulty is resolved in this study by employing ‘stress‐free boundary conditions’ at the open ends of the domains, which have been shown to provide accurate results by a number of previous workers. Copyright © 2005 John Wiley & Sons, Ltd.     

Numerical implementation of a backward probabilistic model of ground water contamination

Backward location and travel time probabilities can be used to characterize known and unknown sources or prior positions of ground water contamination. Backward location probability describes the position of the observed contamination at some time in the past; backward travel time probability describes the amount of time prior to observation that the contamination was released from its source or was at a particular upgradient location. The governing equation for backward probabilities is the adjoint of the governing equation for contaminant transport, but with new load terms. Numerical codes that have been written to solve the forward equations of contaminant transport, e.g., the advection-dispersion equation, can also be used to solve the adjoint equation for location and travel time probabilities; however, the interpretation of the results is different and some new approximations must be made for the load terms. We present the governing equations for backward location and travel time probabilities, and provide appropriate numerical approximations for these load terms using the cell-centered finite difference method, one of the most popular numerical methods in ground water hydrology. We discuss some additional numerical considerations for the backward model including boundary conditions, reversal of the flow field, and interpretation of the results. We illustrate the implementation of the backward probability model using hypothetical examples in one- and two-dimensional domains. We also present a three-dimensional application of a pump-and-treat remediation capture zone delineation at the Massachusetts Military Reservation. The illustrations are performed using MODFLOW-96 for flow simulations and MT3DMS for transport simulations.     

1-D, 2-D, and 3-D analytical solutions of unsaturated flow in groundwater

An excellent tool for checking numerical models of unsaturated flow in groundwater is analytical solutions. However, because of the highly nonlinear nature of the governing partial differential equation, only a limited number of analytical solutions are available. This paper first gives some simple 1-D solutions. Next, by use of a transformation, the nonlinear partial differential equation is converted to a linear one for a specific form of the moisture content vs. pressure head and relative hydraulic conductivity vs. pressure head curves. This allows both 2-D and 3-D solutions to be derived, which is done in this paper. Finally, computations from a finite element computer program are compared with results from one of the analytical solutions to illustrate the use of the derived equations.     

Curvilinear immersed boundary method for simulating coupled flow and bed morphodynamic interactions due to sediment transport phenomena

The fluid-structure interaction curvilinear immersed boundary (FSI-CURVIB) numerical method of Borazjani et al. [3] is extended to simulate coupled flow and sediment transport phenomena in turbulent open-channel flows. The mobile channel bed is discretized with an unstructured triangular mesh and is treated as a sharp-interface immersed boundary embedded in a background curvilinear mesh used to discretize the general channel outline. The unsteady Reynolds-averaged Navier-Stokes (URANS) equations closed with the k − ω turbulence model are solved numerically on a hybrid staggered/non-staggered grid using a second-order accurate fractional step method. The bed deformation is calculated by solving the sediment continuity equation in the bed-load layer using an unstructured, finite-volume formulation that is consistent with the CURVIB framework. Both the first-order upwind and the higher-order hybrid GAMMA schemes [12] are implemented to discretize the bed-load flux gradients and their relative accuracy is evaluated through a systematic grid refinement study. The GAMMA scheme is employed in conjunction with a sand-slide algorithm for limiting the bed slope at locations where the material angle of repose condition is violated. The flow and bed deformation equations are coupled using the partitioned loose-coupling FSI-CURVIB approach [3]. The hydrodynamic module of the method is validated by applying it to simulate the flow in an 180° open-channel bend with fixed bed. To demonstrate the ability of the model to simulate bed morphodynamics and evaluate its accuracy, we apply it to calculate turbulent flow through two mobile-bed open channels, with 90° and 135° bends, respectively, for which experimental measurements are available.