Simulation Of Unsaturated Flow Research Articles

Improved Lumped-Parameter and Numerical Modeling of Unsaturated Water Flow and Stable Water Isotopes.

Characterizing unsaturated water flow in the subsurface is a requirement for understanding effects of droughts on agricultural production or impacts of climate change on groundwater recharge. By employing an improved lumped-parameter model (LPM) approach that mimics variable flow we have interpreted stable water isotope data (δ18 O and δ2 H), taken over 3 years at a lysimeter site located in Germany. Lysimeter soil cores were characterized by sandy gravel (Ly1) and clayey sandy silt (Ly2), and both lysimeters were vegetated with maize. Results were compared with numerical simulation of unsaturated flow and stable water isotope transport using HYDRUS-1D. In addition, both approaches were extended by the consideration of preferential flow paths. Application of the extended LPM, and thus varying flow and transport parameters, substantially improved the description of stable water isotope observations in lysimeter seepage water. In general, findings obtained from the extended LPM were in good agreement to numerical modeling results. However, observations were more difficult to describe mathematically for Ly2, where the periodicity of seasonal stable water isotope fluctuation in seepage water was not fully met by numerical modeling. Furthermore, an extra isotopic upshift improved simulations for Ly2, probably controlled by stable water isotope exchange processes between mobile soil water and quasi-immobile water within stagnant zones. Finally, although LPM requires less input data compared with numerical models, both approaches achieve comparable decision-support integrity. The extended LPM approach can thus be a powerful tool for soil and groundwater management approaches.

An improved Picard iteration scheme for simulating unsaturated flow in porous media

The Richards equation is widely used in the simulation of unsaturated flow and related fields. In the numerical solution process, the finite volume method can be used to carry out numerical discretization, and the Picard method is used for iterative solution. However, in order to obtain a more reliable numerical solution, the space step size of the conventional uniform grid is often small. Especially under some unfavorable numerical conditions, such as infiltration into dry soil and layered soil with very different hydraulic conductivity, this will be very time-consuming and even slow to converge. Thus, combined with the non-uniform grid in the form of Chebyshev, an improved Picard method based on the non-uniform two-grid correction scheme (NTG-PI) is proposed to model 1D unsaturated flow in porous media. Through three examples of unsaturated flow under unfavorable conditions, the proposed scheme was verified. The results show that, compared with the conventional Picard method, NTG-PI can obtain higher numerical accuracy with a smaller number of nodes, and also has a faster convergence rate. This method can provide a certain reference for the numerical simulation of unsaturated flow.

Scaling Relationships between van Genuchten Model Parameters and Hydraulic Conductivity

The soil-water characteristic curve (SWCC) is indispensable for predicting the behaviour of groundwater flow in vadose zones. It remains difficult to determine SWCC, particularly for fractured rocks, under field conditions when little or no test data is available. It is generally believed that the van Genuchten model (VG) parameters α and n are related to pore size distribution (PSD) of soil and can be estimated from PSD data. In this study, correlations are proposed for predicting the VG parameters α and n from hydraulic conductivity K. The VG parameters were obtained by curve-fitting the experimental data of 993 soil samples collected from UNSODA database and published papers. The hydraulic conductivity K of the soil samples varies in six orders of magnitude from 1.30 × 10 −9 m/s to 1.22 × 10 −3 m/s. The statistical distribution of the VG parameters α and n is, respectively, examined using hypothesis test and Anderson-Darling test, showing that α follows a normal distribution after a Box-Cox transformation, and n follows a three-parameter lognormal distribution. The mean and standard deviation of α can be best estimated from K by and σ α = 1.4K 0.13, and the parameter n can be similarly estimated by and σ n = 6.6K 0.18. These correlations are also validated against site-scale test data at Yucca Mountain and published data on fractured rocks and applied to estimate the VG parameters of fractured rocks in the left bank slope at the site of Changbo Hydropower Project. The compiled data and the proposed equations are useful for estimating a narrower range of VG parameters from the more easily-obtained K values at site scale or in case of unavailable test data, and these estimates can be used as input for inverse modelling and simulation of unsaturated flow in porous or fractured media.

High order mimetic difference simulation of unsaturated flow using Richards equation

The vadose zone is the portion of the subsurface above the water table and its pore space usually contains air and water. Due to the presence of infiltration, erosion, plant growth, microbiota, contaminant transport, aquifer recharge, and discharge to surface water, it is crucial to predict the transport rate of water and other substances within this zone. However, flow in the vadose zone has many complications as the parameters that control it are extremely sensitive to the saturation of the media, leading to a nonlinear problem. This flow is referred as unsaturated flow and is governed by Richards equation. Analytical solutions for this equation exists only for simplified cases, so most practical situations require a numerical solution. Nevertheless, the nonlinear nature of Richards equation introduces challenges that causes numerical solutions for this problem to be computationally expensive and, in some cases, unreliable. High order mimetic finite difference operators are discrete analogs of the continuous differential operators and have been extensively used in the fields of fluid and solid mechanics. In this work, we present a numerical approach involving high order mimetic operators along with a Newton root-finding algorithm for the treatment of the nonlinear component. Fully-implicit time discretization scheme is used to deal with the problem’s stiffness.

Soil-water characteristic curves of bentonites in isochoric conditions during wetting: measurement and prediction

The determination of the soil-water characteristic curve (SWCC) of compacted bentonites in the isochoric condition is a prerequisite for unsaturated flow simulations in several geoenvironmental applications. The SWCC data are, however, not readily available for many compacted bentonites over a wide suction range due to difficulties associated with the testing. In this work, wetting SWCCs of four Indian bentonites of different plasticity were established experimentally at compaction dry densities of 1.4, 1.6, and 1.8 Mg/m3 in isochoric conditions using two independent laboratory techniques in different suction ranges. The modified Kovács (MK) model was extended to simulate the SWCC data of the Indian bentonites. Empirical relationships between MK model parameters and compaction dry density for these bentonites were established. A theoretical procedure for the estimation of SWCCs of the compacted bentonites based on the basic bentonite properties and proposed correlations was proposed and validated on the available wetting SWCC data of compacted GaoMiaoZi (GMZ), MX80, and FEBEX bentonites from the literature studies. The proposed method predicts the wetting SWCCs of compacted powder bentonites, in volume-restrained condition, from the basic properties satisfactorily. The proposed model is an approximate method for predicting the SWCCs of bentonites in the absence of measured data and in understanding the influence of dry density and plasticity.

Effective flow and transport properties of heterogeneous unsaturated soils

The heterogeneity of field scale soils poses a challenge to predictive large scale flow and transport modeling. The theory of effective macroscale parameters holds good and is applicable in dealing with such problems. But the validity of the analytic stochastic solutions obtained for randomly heterogeneous soils is debatable, as the test cases under which they are validated are of limited scope due to linearization and perturbation approximations. In this study, samples of heterogeneous soils are generated using sets of spatially correlated random field parameters that are either geometrically isotropic, or else, geometrically anisotropic with either horizontal or vertical stratification (perfect or imperfect). Several combinations of ratios of correlation length and capillary dispersion lengths are considered. Numerical simulations of unsaturated flow are performed on each randomly heterogeneous soil sample. The principal components K^ii(Ψ) of the macroscale effective unsaturated conductivity are then obtained as a function of the mean suction Ψ of the sample. They are compared to stochastic spectral perturbation theory, and to a probabilistic semi-empirical Power Average Model (PAM). They are also compared with arithmetic, geometric and harmonic mean conductivity-suction curves. The numerically upscaled principal conductivity curves match quite well the PAM, better than the classical means (Arithmetic, Geometric, Harmonic), and also somewhat better than the curves obtained from stochastic spectral perturbation theory. It is observed that the upscaled principal components Kii(ψ), obtained numerically and with the PAM along directions “i” orthogonal/parallel to perfect stratification coincide with the harmonic/arithmetic mean curves at low suctions (i.e., near saturation), but deviate from it and come closer to the geometric mean at higher suctions. The PAM appears suitable for generation of approximate upscaled conductivity curves, e.g., for obtaining the mesh-scale or block-scale conductivity curves in large scale simulation codes. Transient solute transport simulations are then performed on the detailed random velocity fields obtained from the steady state simulations of unsaturated flow in the randomly heterogeneous soil samples. Snapshots of solute concentration C(x,z,t) are taken at different times. The temporal evolution of spatial moments of concentration is analyzed in order to characterize the macroscale advection and dispersion of the unsaturated concentration plume, and in particular, its macro-dispersion coefficient (D) and dispersivity length scale (A). For the synthetic soil samples considered in this study, the macro-dispersive spreading of the solute is stronger for flow parallel to vertical stratification, compared to flow perpendicular to horizontal stratification, and also, compared to flow in statistically isotropic non-stratified soil.

Unsaturated flow in a packing of swelling particles; a grain-scale model

In this work, a grain-scale modelling technique is introduced for the simulation of unsaturated flow in deforming and swelling granular materials. To do so, a pore-scale model for unsaturated flow is coupled to the discrete element method (DEM). It is assumed that initially a dry packing of particles is quickly invaded by a liquid and becomes fully saturated. Particles start absorbing the liquid and this causes a rearrangement of particles, entrance of air into the packing, and a redistribution of the liquid (i.e. unsaturated flow). Flow was computed using a scheme of implicit pressure solver and explicit saturation update (IMPES), whilst particle movement was modelled using DEM. Simulations are continued until the packing is dry again. This is the first time that such a pore-scale model has been developed. We have used the model to investigate unsaturated flow during drying of a bed of swelling particles. Results indicated that the characteristic time scales of unsaturated flow and water absorption determine the swelling behaviour of the particle packing.

Investigation of Gravity‐Driven Infiltration Instabilities in Smooth and Rough Fractures Using a Pairwise‐Force Smoothed Particle Hydrodynamics Model

Core Ideas We study infiltration dynamics in smooth and rough fractures using a 3D SPH model. Roughness induces the formation of flow fingering and preferential flow paths. Fractal roughness leads to a deceleration of average fingertip velocities. The velocity variance of fingers positively correlates with the degree of roughness. Normalized finger velocity and length scale in accordance with analytical predictions. This work investigates small‐scale infiltration dynamics in smooth and rough single fractures using a three‐dimensional multiphase pairwise‐force smoothed particle hydrodynamics (PF‐SPH) model. Gravity‐driven infiltration instabilities in fractures under unsaturated conditions can significantly influence the arrival time of tracers or contaminants, and the rapid and localized recharge dynamics in fractured–porous aquifer systems. Here, we study the influence of roughness and injection rate on fluid flow modes and flow velocity. Three types of fractures are considered with different degrees of roughness, including a smooth fracture. Both the rough and smooth fractures exhibit flow instabilities, fingering, and intermittent flow regimes for low infiltration rates. In agreement with theoretical predictions, a flat fluid front is achieved when the flux q supplied to a fracture is larger than the gravitationally driven saturated flux [q > kρg/μcos(φ), where k is the intrinsic permeability of the fracture, ρ is a density, μ is the viscosity, and φ is the fracture inclination angle measured from the vertical direction]. To characterize the flow instability, we calculate standard deviations of velocity along the fracture width. For the considered infiltration rates, we find that an increase in roughness decreases the flow velocity and increases the standard deviation of velocity. This is caused by a higher likelihood of flow discontinuities in the form of fingering and/or snapping rivulets. To validate our unsaturated flow simulations in fractures, we estimate the scaling of specific discharge with normalized finger velocity, compute the relationship between fingertip length and scaled finger velocity, and find good agreement with experimental results.

A numerical formulation with unified unilateral boundary condition for unsaturated flow problems in porous media

This paper proposed a numerical formulation for unsaturated flow problems with nonlinear boundaries of seepage face and soil–atmosphere interface via the concept of parabolic variational inequality (PVI) method. A unified unilateral boundary condition was first proposed to represent the conditions on the seepage face and soil–atmosphere interface boundaries within the partial differential equation (PDE) formulation. A PVI formulation mathematically equivalent to the PDE formulation was then proposed, which automatically transforms the flux part of the unified unilateral boundary condition into the natural boundary condition and eliminates the singularity at seepage points. By discretizing the PVI formulation, a finite element procedure together with an iterative algorithm was suggested. An existing experiment of unsaturated flow in a layered hillside and a laboratory test of unsaturated flow through sand flume performed in this study were used to validate the proposed method, with a good agreement between the measured and computed results and a satisfactory balance of mass being maintained during the simulations. The numerical results also indicated that the problem of mesh dependence associated with unsaturated flow simulations is well addressed with the proposed numerical method. Finally, the process of unsaturated flow in a soil slope with layers of horizontal drains subjected to rainfall/evaporation was further examined. The numerical results reveal that the deployment of drains in a soil slope can significantly lower the pore water pressure around the drains, with the bottom layer drains being most effective in controlling the seepage flow.

Volume Change Consideration in Determining Unsaturated Soil Properties for Geotechnical Applications

AbstractA common underlying assumption in the development of soil-water characteristic curves (SWCCs) is that the soil is sufficiently stiff such that no changes occur in the volume of the soil associated with soil-suction change. A number of researchers recently recognized the importance of considering soil-volume changes in SWCCs; however, few researchers have evaluated the effect of volume change on the computed rate and degree of wetting/drying of unsaturated soils and the associated impact on computed soil-volume change response. This paper focuses on the impacts of considering soil-volume change in the determination of SWCCs for high-volume-change soils, with emphasis on the estimated unsaturated hydraulic conductivity function and response to wetting. Numerical simulations of unsaturated flow and suction-induced volume change were used to study the effect of consideration of volume change in the determination of SWCCs. Numerical analyses used laboratory-determined SWCCs with and without the inclusi...

Numerical simulation of water flow in a large waste rock pile, Haizhou coal mine, China

Waste rock is the large volume of broken rock and low grade ore materials tend to create potential problems, such as acid mine drainage (AMD), leaching of heavy metals, and slope stability problems. These problems are related to the flow of water through the waste rock. In this study, long term numerical simulations of unsaturated flow in a large waste rock pile are conducted to investigate the effect of internal variability in water flow and saturation. The approach first consists of a basic geotechnical properties investigation to define the main features of the pile. This is followed by numerical simulations of the unsaturated flow through the pile based on the results from the field and laboratory characterization. The simulation results showed that pile is unable to fully saturate and drain during the simulated period. It was also found that high evaporation and the high bottom flux rate and also coarse grain nature of waste rock materials results holding few amounts of water in the system.

Storage and mobilization of natural and septic nitrate in thick unsaturated zones, California

Summary Mobilization of natural and septic nitrate from the unsaturated zone as a result of managed aquifer recharge has degraded water quality from public-supply wells near Yucca Valley in the western Mojave Desert, California. The effect of nitrate storage and potential for denitrification in the unsaturated zone to mitigate increasing nitrate concentrations were investigated. Storage of water extractable nitrate in unsaturated alluvium up to 160 meters (m) thick, ranged from 420 to 6600 kilograms per hectare (kg/ha) as nitrogen (N) beneath undeveloped sites, from 6100 to 9200 kg/ha as N beneath unsewered sites. Nitrate reducing and denitrifying bacteria were less abundant under undeveloped sites and more abundant under unsewered sites; however, δ 15 N–NO 3 , and δ 18 O–NO 3 data show only about 5–10% denitrification of septic nitrate in most samples—although as much as 40% denitrification occurred in some parts the unsaturated zone and near the top of the water table. Storage of nitrate in thick unsaturated zones and dilution with low-nitrate groundwater are the primary attenuation mechanisms for nitrate from septic discharges in the study area. Numerical simulations of unsaturated flow, using the computer program TOUGH2, showed septic effluent movement through the unsaturated zone increased as the number and density of the septic tanks increased, and decreased with increased layering, and increased slope of layers, within the unsaturated zone. Managing housing density can delay arrival of septic discharges at the water table, especially in layered unsaturated alluvium, allowing time for development of strategies to address future water-quality issues.

Effect of Dense Material Layers on Unsaturated Water Flow Inside a Large Waste Rock Pile: A Numerical Investigation

The construction method used to build waste rock piles influences their internal structure. Commonly used methods typically lead to the creation of compacted material layers within otherwise loose, coarse-grained waste rock. These dense layers, which typically have a finer grain size, affect the movement and distribution of water inside the pile. Long-term numerical simulations of unsaturated flow in a large pile were conducted to investigate the effect of such layers. The simulations led to various observations that provide a better understanding of the hydrogeological behaviour of the modeled pile (based on an actual case). The results show how water distribution and seepage within the pile are influenced by the presence of these layers. Other factors, including the magnitude of precipitation (or recharge) and pile size, were also investigated. This article presents the main results of the simulations, with some comments on their practical implications for pile design.

Performance assessment of exponential Rosenbrock methods for large systems of ODEs

This article studies time integration methods for stiff systems of ordinary differential equations of large dimension. For such problems, implicit methods generally outperform explicit methods because the step size is usually less restricted by stability constraints. Recently, however, a family of explicit methods, called exponential integrators, have become popular for large stiff problems due to their favourable stability properties and the rapid convergence of non-preconditioned Krylov subspace methods for computing matrix-vector products involving exponential-like functions of the Jacobian matrix. In this article, we implement the so-called exponential Rosenbrock methods using Krylov subspaces. Numerical experiments on a challenging real-world test problem reveal that these methods are a promising preconditioner-free alternative to well-established approaches based on preconditioned Newton--Krylov implementations of the backward differentiation formulas. References E. J. Carr, T. J. Moroney and I. W. Turner. Efficient simulation of unsaturated flow using exponential time integration. Appl. Math. Comput., 217(14): 6587--6596, 2011. doi:10.1016/j.amc.2011.01.041 E. J. Carr, I. W. Turner and M. Ilic. Krylov subspace approximations for the exponential Euler method: error estimates and the harmonic Ritz approximant. ANZIAM Journal, 52: C612--C627, 2011. http://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/view/3938 E. J. Carr, I. W. Turner and P. Perre. A variable step size Jacobian-free exponential integrator for simulating transport in heterogeneous porous media: application to wood drying. J. Comput. Phys., 233: 66-82, 2013. doi:10.1016/j.jcp.2012.07.024 Sundials user documentation. https://computation.llnl.gov/casc/sundials/documentation/documentation.html M. Hochbruck, C. Lubich and H. Selhofer. Exponential integrators for large systems of differential equations. SIAM J. Sci. Comput., 19(5): 1552--1574, 1998. doi:10.1137/S1064827595295337 M. Hochbruck, A. Ostermann and J. Schweitzer. Exponential Rosenbrock-type methods. SIAM J. Numer. Anal., 47(1): 786--803, 2009. doi:10.1137/080717717 C. T. Kelley. Solving nonlinear equations with Newton's method. SIAM, USA, 2003. D. A. Knoll and D. E. Keyes. Jacobian-free Newton--Krylov methods: a survey of approaches and applications. J. Comput. Phys., 193(2): 357--397, 2004. doi:10.1016/j.jcp.2003.08.010 M. Tokman. Efficient integration of large stiff systems of ODEs with exponential propagation iteration (EPI) methods. J. Comput. Phys., 213: 748--776, 2006. doi:10.1016/j.jcp.2005.08.032

Modeling Sedimentation-Filtration Basins for Urban Watersheds Using Soil and Water Assessment Tool

AbstractSedimentation-filtration (SedFil) basins are one of the storm-water best management practices (BMPs) that are intended to mitigate water quality problems in urban creeks and rivers. A new physically based model of variably saturated flows was developed for simulating flow and sediment in SedFils within the Soil and Water Assessment Tool (SWAT). The integrated SWAT-SedFil model allows for simulation of unsaturated flow in the filtration basin during small storms and fully saturated flow. Unsaturated flow was modeled using a modified Green and Ampt equation, and saturated flow was simulated with Darcy’s Law. Unsaturated flow comprises only a small fraction of large storm events; however, many regular storms are small and may not generate sufficient runoff to create a saturated flow in the filtration basin. Therefore, the combined unsaturated/saturated flow approach for modeling SedFils improved the accuracy of the model, especially in long-term evaluations. The model performs well with respect to es...

Deep subsurface drip irrigation using coal-bed sodic water: Part I. Water and solute movement

Water co-produced with coal-bed methane (CBM) in the semi-arid Powder River Basin of Wyoming and Montana commonly has relatively low salinity and high sodium adsorption ratios that can degrade soil permeability where used for irrigation. Nevertheless, a desire to derive beneficial use from the water and a need to dispose of large volumes of it have motivated the design of a deep subsurface drip irrigation (SDI) system capable of utilizing that water. Drip tubing is buried 92cm deep and irrigates at a relatively constant rate year-round, while evapotranspiration by the alfalfa and grass crops grown is seasonal. We use field data from two sites and computer simulations of unsaturated flow to understand water and solute movements in the SDI fields. Combined irrigation and precipitation exceed potential evapotranspiration by 300–480mm annually. Initially, excess water contributes to increased storage in the unsaturated zone, and then drainage causes cyclical rises in the water table beneath the fields. Native chloride and nitrate below 200cm depth are leached by the drainage. Some CBM water moves upward from the drip tubing, drawn by drier conditions above. Chloride from CBM water accumulates there as root uptake removes the water. Year over year accumulations indicated by computer simulations illustrate that infiltration of precipitation water from the surface only partially leaches such accumulations away. Field data show that 7% and 27% of added chloride has accumulated above the drip tubing in an alfalfa and grass field, respectively, following 6 years of irrigation. Maximum chloride concentrations in the alfalfa field are around 45cm depth but reach the surface in parts of the grass field, illustrating differences driven by crop physiology. Deep SDI offers a means of utilizing marginal quality irrigation waters and managing the accumulation of their associated solutes in the crop rooting zone.

Simulation of saturated and unsaturated flow in karst systems at catchment scale using a double continuum approach

Abstract. The objective of this work is the simulation of saturated and unsaturated flow in a karstified aquifer using a double continuum approach. The HydroGeoSphere code (Therrien et al., 2006) is employed to simulate spring discharge with the Richards equations and van Genuchten parameters to represent flow in the (1) fractured matrix and (2) conduit continuum coupled by a linear exchange term. Rapid vertical small-scale flow processes in the unsaturated conduit continuum are accounted for by applying recharge boundary conditions at the bottom of the saturated model domain. An extensive sensitivity analysis is performed on single parameters as well as parameter combinations. The transient hydraulic response of the karst spring is strongly controlled by the matrix porosity as well as the van Genuchten parameters of the unsaturated matrix, which determine the head dependent inter-continuum water transfer when the conduits are draining the matrix. Sensitivities of parameter combinations partially reveal a non-linear dependence over the parameter space. This can be observed for parameters not belonging to the same continuum as well as combinations, which involve the exchange parameter, showing that results of the double continuum model may depict a certain degree of ambiguity. The application of van Genuchten parameters for simulation of unsaturated flow in karst systems is critically discussed.

On the importance of hysteresis and heterogeneity in the numerical simulation of unsaturated flow

Soil hydraulic properties in field scale porous media are distributed with large spatial heterogeneity. Thus, detailed understanding of micro-physical properties leading to macro-scale heterogeneity is important for generalizing purposes in modeling. Hysteretic water retention properties greatly affect solute transport under unsaturated condition. In this study, laboratory infiltration experiments were carried out to investigate the importance of hysteresis in numerical modeling of water flow under heterogeneous distribution of hydraulic conductivity. Soil water content was observed by time domain reflectometry using small printed circuit board probes (PCBP). The PCBP observations were compared to results from numerical simulations involving saturated-unsaturated seepage analysis and two-phase flow with/without hysteresis effects. Observations were compared to both uniform and heterogeneous parameter modeling. Differences between unsaturated-saturated and two-phase flow were quite small in the numerical results for soil water content. The uniform parameter model, however, could not properly reproduce experimental water content change. The numerical results of water content change including hysteresis matched better observations for the unsaturated-saturated and two-phase flow cases. Results displayed that it is important to take heterogeneity and hysteresis into consideration in the numerical modeling in both two-phase and saturated-unsaturated flow analysis. (Less)

Krylov subspace approximations for the exponential Euler method: error estimates and the harmonic Ritz approximant

We study Krylov subspace methods for approximating the matrix-function vector product $\varphi(tA)b$ where $\varphi(z) = [\exp(z)-1]/z$. This product arises in the numerical integration of large stiff systems of differential equations by the Exponential Euler Method, where $A$~is the Jacobian matrix of the system. Recently, this method has found application in the simulation of transport phenomena in porous media within mathematical models of wood drying and groundwater flow. We develop an a posteriori upper bound on the Krylov subspace approximation error and provide a new interpretation of a previously published error estimate. This leads to an alternative Krylov approximation to $\varphi(tA)b$, the so-called Harmonic Ritz approximant, which we find does not exhibit oscillatory behaviour of the residual error. References E. J. Carr, T. J. Moroney and I. W. Turner. Efficient simulation of unsaturated flow using exponential time integration. Appl. Math. Comput., 217(14):6587--6596, 2011. doi:10.1016/j.amc.2011.01.041 E. J. Carr, I. W. Turner and P. Perre. A Jacobian-free exponential integrator for simulating transport in heterogeneous porous media: application to the drying of softwood, submitted for publication. E. J. Carr, I. W. Turner and P. Perre. A new control-volume finite-element scheme for heterogeneous porous media:application to the drying of softwood. Chem. Eng. Technol., 34(7):1143--1150, 2011. doi:10.1002/ceat.201100060 E. Celledoni and I. Moret. A Krylov projection method for systems of ODEs. Appl. Numer. Math., 24(2--3):365--378, 1997. doi:10.1016/S0168-9274(97)00033-0 N. J. Higham. Functions of matrices: theory and computation. SIAM, Philadelphia, PA, USA, 2008. M. Hochbruck, M. E. Hochstenbach. Subspace extraction for matrix functions, submitted for publication. http://www.win.tue.nl/ hochsten/pdf/funext.pdf M. Hochbruck, C. Lubich and H. Selhofer. Exponential integrators for large systems of differential equations. SIAM J. Sci. Comput., 19(5):1552--1574, 1998. doi:10.1137/S1064827595295337 B. V. Minchev and W. M. Wright. A review of exponential integrators for first order semi-linear problems. Numerics No. 2/05, Norwegian University of Science and Technology, Trondheim, Norway, 2005. http://www.math.ntnu.no/preprint/numerics/2005/N2-2005.ps P. Perre and I. Turner. A heterogeneous wood drying computational model that accounts for material property variation across growth rings. Chem. Eng. J., 86(1--2):117--131, 2002. doi:10.1016/S1385-8947(01)00270-4

Efficient simulation of unsaturated flow using exponential time integration

We assess the performance of an exponential integrator for advancing stiff, semidiscrete formulations of the unsaturated Richards equation in time. The scheme is of second order and explicit in nature but requires the action of the matrix function φ( A) = A −1( e A − I) on a suitability defined vector v at each time step. When the matrix A is large and sparse, φ( A) v can be approximated by Krylov subspace methods that require only matrix–vector products with A. We prove that despite the use of this approximation the scheme remains second order. Furthermore, we provide a practical variable-stepsize implementation of the integrator by deriving an estimate of the local error that requires only a single additional function evaluation. Numerical experiments performed on two-dimensional test problems demonstrate that this implementation outperforms second-order, variable-stepsize implementations of the backward differentiation formulae.