Groundwater and Surface Water Interactions: Controlling Water Levels and Aquifer Recharge, Apuntes de Ciencias de la Educación

¡Descarga Groundwater and Surface Water Interactions: Controlling Water Levels and Aquifer Recharge y más Apuntes en PDF de Ciencias de la Educación solo en Docsity! Modeling Groundwater Flow and Contaminant Transport Prof. Dr. Jacob Bear Technion-Israel Institute of Technology Dept. Civil & Environmental Engineering 32000 Haifa Israel cvrbear@tx.technion.ac.il Dr. Alexander H.-D. Cheng University of Mississippi Dept. Civil Engineering P.O.Box 1848 University MS 38677-1848 USA acheng@olemiss.edu ISBN 978-1-4020-6681-8 e-ISBN 978-1-4020-6682-5 DOI 10.1007/978-1-4020-6682-5 Springer Dordrecht Heidelberg London New York Library of Congress Control Number: 2009938711 c© Springer Science+Business Media B.V. 2010 No part of this work may be reproduced, stored in a retrieval system, or transmitted in any form or by any means, electronic, mechanical, photocopying, microfilming, recording or otherwise, without written permission from the Publisher, with the exception of any material supplied specifically for the purpose of being entered and executed on a computer system, for exclusive use by the purchaser of the work. Printed on acid-free paper Springer is part of Springer Science+Business Media (www.springer.com) Contents Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi List of Main Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Groundwater in Water Resources Systems . . . . . . . . . . . . . . . . . 2 1.1.1 Hydrological cycle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Surface water versus groundwater . . . . . . . . . . . . . . . . . . 3 1.1.3 Characteristics of groundwater . . . . . . . . . . . . . . . . . . . . . 5 1.1.4 Functions of aquifers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.1.5 Subsurface contamination . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.1.6 Sustainable yield . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 1.2 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 1.2.1 Modeling concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 1.2.2 Modeling process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 1.2.3 Model use . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 1.3 Continuum Approach to Transport in Porous Media . . . . . . . . 42 1.3.1 Phases, chemical species and components . . . . . . . . . . . . 42 1.3.2 Need for continuum approach . . . . . . . . . . . . . . . . . . . . . . 43 1.3.3 Representative elementary volume and averages . . . . . . 45 1.3.4 Scale of heterogeneity in continuum models . . . . . . . . . . 48 1.3.5 Homogenization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 1.4 Scope and Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 2 GROUNDWATER AND AQUIFERS . . . . . . . . . . . . . . . . . . . . . 65 2.1 Definitions of Aquifers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 2.2 Moisture Distribution in Vertical Soil Profile . . . . . . . . . . . . . . . 66 2.3 Classification of Aquifers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 2.4 Solid Matrix Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 2.4.1 Soil classification based on grain size distribution . . . . . 71 2.4.2 Porosity and void ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 2.4.3 Specific surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 2.5 Inhomogeneity and Anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 v vi 2.6 Hydraulic Approach to Flow in Aquifers . . . . . . . . . . . . . . . . . . . 78 3 REGIONAL GROUNDWATER BALANCE . . . . . . . . . . . . . . 81 3.1 Groundwater Flow and Leakage . . . . . . . . . . . . . . . . . . . . . . . . . . 82 3.1.1 Inflow and outflow through aquifer boundaries . . . . . . . 82 3.1.2 Leakage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.2 Natural Replenishment from Precipitation . . . . . . . . . . . . . . . . . 84 3.3 Return Flow from Irrigation and Sewage . . . . . . . . . . . . . . . . . . 88 3.4 Artificial Recharge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 3.4.1 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 3.4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 3.5 River-Aquifer Interrelationships . . . . . . . . . . . . . . . . . . . . . . . . . . 97 3.6 Springs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.7 Evapotranspiration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 3.8 Pumping and Drainage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 3.9 Change in Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 3.10 Regional Groundwater Balance . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4 GROUNDWATER MOTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.1 Darcy’s Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.1.1 Empirical law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.1.2 Extension to three-dimensional space . . . . . . . . . . . . . . . 116 4.1.3 Hydraulic conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 4.1.4 Extension to anisotropic porous media . . . . . . . . . . . . . . 120 4.2 Darcy’s Law as Momentum Balance Equation . . . . . . . . . . . . . . 124 4.2.1 Darcy’s law by volume averaging . . . . . . . . . . . . . . . . . . . 125 4.2.2 Darcy’s law by homogenization . . . . . . . . . . . . . . . . . . . . . 128 4.2.3 Effective hydraulic conductivity by homogenization . . . 140 4.3 Non-Darcy Laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 4.3.1 Range of validity of Darcy’s law . . . . . . . . . . . . . . . . . . . . 145 4.3.2 Non-Darcian motion equations . . . . . . . . . . . . . . . . . . . . . 147 4.4 Aquifer Transmissivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 4.5 Dupuit Assumption for Phreatic Aquifer . . . . . . . . . . . . . . . . . . 152 5 WATER BALANCE AND COMPLETE FLOW MODEL . 161 5.1 Mass Balance Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162 5.1.1 Fundamental mass balance equations . . . . . . . . . . . . . . . 162 5.1.2 Deformable porous medium . . . . . . . . . . . . . . . . . . . . . . . . 167 5.1.3 Specific storativity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 5.1.4 Flow equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 5.2 Initial and Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 181 5.2.1 Boundary surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 5.2.2 Initial and general boundary conditions . . . . . . . . . . . . . 185 5.2.3 Particular boundary conditions . . . . . . . . . . . . . . . . . . . . . 187 5.3 Complete 3-D Mathematical Flow Model . . . . . . . . . . . . . . . . . . 203 Contents ix 7.10.2 Caps and cutoff walls . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 515 7.10.3 Pump-and-treat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 516 7.10.4 Soil vapor extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 519 7.10.5 Air sparging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 519 7.10.6 Permeable reactive barrier . . . . . . . . . . . . . . . . . . . . . . . . . 522 8 NUMERICAL MODELS AND COMPUTER CODES . . . . 525 8.1 Finite Difference Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 527 8.1.1 Laplace equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 528 8.1.2 Diffusion equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 531 8.1.3 Cell-centered approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533 8.1.4 Boundary and boundary conditions . . . . . . . . . . . . . . . . . 535 8.2 Finite Volume Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 537 8.3 Finite Element Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 541 8.3.1 Weighted residual methods . . . . . . . . . . . . . . . . . . . . . . . . 542 8.3.2 Galerkin finite element methods . . . . . . . . . . . . . . . . . . . . 550 8.3.3 Meshless finite element methods . . . . . . . . . . . . . . . . . . . . 558 8.3.4 Control volume finite element methods . . . . . . . . . . . . . . 559 8.4 Boundary Element Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 560 8.5 Radial Basis Function Collocation Methods . . . . . . . . . . . . . . . . 564 8.6 Eulerian-Lagrangian Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 570 8.6.1 Lagrangian method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 570 8.6.2 Method of characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . 574 8.6.3 Random walk method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 576 8.6.4 Modified Eulerian-Lagrangian method . . . . . . . . . . . . . . 578 8.7 Matrix Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 579 8.7.1 Conjugate gradient method . . . . . . . . . . . . . . . . . . . . . . . . 580 8.7.2 Preconditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 581 8.8 Computer Codes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 583 9 SEAWATER INTRUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593 9.1 Occurrence and Exploration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593 9.1.1 Occurrence of seawater intrusion . . . . . . . . . . . . . . . . . . . 593 9.1.2 Exploration of saltwater intrusion . . . . . . . . . . . . . . . . . . 596 9.2 Sharp Interface Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 601 9.2.1 Sharp interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 601 9.2.2 Ghyben-Herzberg approximation . . . . . . . . . . . . . . . . . . . 605 9.2.3 Upconing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 607 9.2.4 Essentially horizontal flow model . . . . . . . . . . . . . . . . . . . 610 9.2.5 Some analytical solutions for stationary interface . . . . . 613 9.2.6 Multilayered aquifers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 620 9.3 Transition Zone Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 620 9.3.1 Variable density model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 621 9.3.2 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 628 9.4 Management of Coastal Aquifer . . . . . . . . . . . . . . . . . . . . . . . . . . 633 Contents x 10 MODELING UNDER UNCERTAINTY . . . . . . . . . . . . . . . . . . 637 10.1 Stochastic Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 639 10.1.1 Random process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 639 10.1.2 Quantifying uncertainty as stochastic process . . . . . . . . 640 10.1.3 Ensemble statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 643 10.1.4 Spatial (or temporal) statistics . . . . . . . . . . . . . . . . . . . . . 648 10.1.5 Ergodicity hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 650 10.2 Tools for Uncertainty Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 652 10.2.1 Kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653 10.2.2 Sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 661 10.2.3 Monte Carlo simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 666 10.2.4 Generation of random field . . . . . . . . . . . . . . . . . . . . . . . . 670 10.3 Examples of Uncertainty Problems . . . . . . . . . . . . . . . . . . . . . . . 676 10.3.1 Random boundary conditions . . . . . . . . . . . . . . . . . . . . . . 678 10.3.2 Uncertain parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 686 11 OPTIMIZATION, INVERSE, AND MANAGEMENT TOOLS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 695 11.1 Groundwater Management . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 696 11.2 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 698 11.2.1 Optimization problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 698 11.2.2 Linear programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 700 11.2.3 Nonlinear problems and unconstrained optimization . . 712 11.2.4 Gradient search method . . . . . . . . . . . . . . . . . . . . . . . . . . . 715 11.2.5 Genetic algorithm and simulated annealing . . . . . . . . . . 720 11.2.6 Chance constrained optimization . . . . . . . . . . . . . . . . . . . 726 11.2.7 Multiobjective optimization . . . . . . . . . . . . . . . . . . . . . . . . 728 11.3 Inverse Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 742 11.3.1 Pumping test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 745 11.3.2 Regional scale parameter estimation . . . . . . . . . . . . . . . . 755 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 759 Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 815 Contents Preface Groundwater, extracted from deep geological formations (called aquifers) through pumping wells, constitutes an important component of many wa- ter resource systems. A spring constitutes an outlet for groundwater from an underlying aquifer to ground surface; its discharge rate may be strongly af- fected by pumping from the same aquifer in the vicinity of the spring. Water can be injected through specially designed wells into an aquifer, say, for stor- age purposes. A water table aquifer can also be artificially recharged through infiltration ponds. These are just a few examples of factor that may affect the management of a groundwater system. Decisions associated with such management include, for example, • The volume that can be safely withdrawn annually from the aquifer. • The location of pumping and artificial recharge wells, and their rates. • The quality of the water to be maintained in the aquifer, and/or to be pumped from it. In fact, in the management of water resources, the quantity and quality problems cannot be separated from each other. In many parts of the world, as a result of increased withdrawal of groundwater, often beyond permissi- ble limits, the quality of groundwater has been continuously deteriorating, causing much concern to both suppliers and users. The quality deterioration may manifest itself in the form of an increase in the total salinity, or as in- creased concentrations of nitrates and other undesirable chemical species, or as increased concentrations of harmful pathogens. Traditionally, hydrogeologist dealt with flow in aquifers, and with certain water quality aspects, e.g., salinization. Soil physicist and agronomists, in connection with agricultural activities, have modeled the movement of water and chemicals (e.g., fertilizers) in the unsaturated zone. The hydrogeologist, whose primary interest has been water in aquifers, regarded the unsaturated zone only as the domain through which water from precipitation passes on its way to replenish an underlying water table aquifer. The details of the actual movement of water through the unsaturated zone have been of little or no interest. The situation has completely changed with the rising interest in subsurface contamination. Clearly, interdisciplinary efforts, straddling nu- xi xiv Preface Three additional topics, strongly related to use of models for predicting flow and transport regimes in aquifers, within the framework of management, are discussed: • The use of numerical models and computer codes as practical tools for solving the mathematical models. • The issues of uncertainty associated with modeling. • Certain mathematical tools for groundwater management. With these objectives in mind, the book is aimed at practitioners, model- ers, water resources managers, scientists, and researchers, who face the need to build and solve models of flow and contaminant transport in the subsur- face. It is also suitable for graduate and upper level undergraduate students who are interested in such topics as groundwater, water resources, and envi- ronmental engineering. The basic scientific background needed is the concepts and terminologies of hydrology and hydraulics. Finally, we wish to acknowledge the following colleagues who have provided useful advice for various parts of the book: Yunwei Sun and Walt McNab of Lawrence Livermore National Laboratory; Quanlin Zhou of Lawrence Berke- ley National Laboratory; Peter Lichtner of Los Alamos National Labora- tory; Vicky Freeman of Pacific Northwest National Laboratory; Ajit Sadana of the University of Mississippi; Randy Gentry of the University of Ten- nessee; Shlomo Neuman of the University of Arizona; Prabhakar Clement of Auburn University; T.N. Narasimhan of the University of California, Berke- ley; Shlomo Orr of MRDS, Inc.; Shaul Sorek of Ben-Gurion University, Israel; Uri Shavit and Leonid Fel of Technion, Israel; Jacob Bensabat of Environmen- tal and Water Resources Engineering Inc., Israel; Dalila Loudyi of Hassan II University, Morocco; and Don Nield of University of Auckland, New Zealand. Jacob Bear Haifa, Israel Alexander H.-D. Cheng Oxford, Mississippi, USA 2009 List of Main Symbols a Dispersivity of a porous medium (fourth rank tensor). aijkl Component of a. aL, aT Longitudinal and transversal dispersivities (isotropic medium), respectively. B Thickness of a confined aquifer. B′ Thickness of a semipervious layer (aquitard). B Mass balance operator. cab, caa Covariance between a and b, and autocovariance of a. cr Hydraulic resistance of a semipervious layer. cγα Mass concentration of a γ-species in an α-phase. [cγα] Molar concentration of a γ-species in an α-phase. C Normalized salt mass fraction. Cv Consolidation coefficient. Cw Moisture (water) capacity. d Effective grain diameter. Microscopic length characterizing the void space. d∗ Length characterizing macroscopic heterogeneity. Da Darcy number. Dm Damköhler number. Dγα Coefficient of dispersion of the mass of a γ-species in an α-phase. Dγαh Coefficient of hydrodynamic dispersion of a γ-species in an α-phase (= D∗γα + Dγα). DL, DT Longitudinal and transversal dispersion coefficients. Dγα Coefficient of molecular diffusion of a γ-species in an α-phase. D∗γα Coefficient of molecular diffusion of a γ-species in an α-phase in a porous medium (= DγαT∗). Dw Moisture diffusivity. e An intensive quantity, density of an extensive quantity, E. Void ratio (= Uov/Uos). E An extensive quantity. eγα Density of Eγα (= E of γ in α per unit volume of α-phase). Eγα An extensive quantity, E, of a γ-species in an α-phase. fEα→β Rate of transfer of E from an α-phase to a β-phase, across their xv xvi List of Main Symbols common microscopic interface, per unit volume of porous medium. fα Fractional flow rate of an α-phase in two-phase flow. F Concentration of a species adsorbed on a solid (= mγs/ms). Function representing a surface. F Body force. F Faraday’s constant. Fo Fourier number. Fr Froude number. g Gravity acceleration. G Gibbs free energy. Green’s function. h Piezometric head. h∗ Hubbert’s potential. hc Height of capillary fringe. hr Relative humidity. H Height of water table in phreatic aquifer. Heaviside unit step function. H Henry’s coefficient. H Hessian matrix. I Rate of infiltration. Ionic strength of solution. Integral scale. I Unit tensor. jγ Microscopic diffusive mass flux of a γ-species, relative to the mass weighted velocity (= ργ(Vγ −V)). jtE γ α Total microscopic flux of Eγα (= e γ αV Eγα). Jγα Macroscopic diffusive flux of a γ-species in an α-phase (≡ jγ α ). JE Dispersive flux of E. JtE γ α Total macroscopic flux of Eγα. JEαh Sum of diffusive and dispersive fluxes of E in an α-phase. J Hydraulic gradient. k Permeability (scalar). k Permeability (second rank tensor). kα Effective permeability of an α-phase. krα Relative permeability of an α-phase, in an isotropic porous medium. kγα Degradation rate constant of a γ-species in an α-phase. K Hydraulic conductivity (scalar). K Hydraulic conductivity (second rank tensor). K′ Hydraulic conductivity of semi-pervious layer. Kd Partitioning coefficient. Keq Equilibrium coefficient. Ksp Solubility product.  Characteristic size (diameter) of REV. L Length of a column. List of Main Symbols xix homogenization process. z Vertical coordinate (positive upward). zγ Electrical charge of ion of γ species. Z Objective function for optimization. Greek letters α Symbol for an α-phase. Coefficient of soil compressibility. α′ Coefficient of rock compressibility. αE Transfer coefficient of E. βp, βT Fluid compressibilities at constant p and constant T , respectively. γ Symbol denoting a γ-species. Euler constant (= 0.5772157 . . .). Specific weight. Semivariogram. γ(x) Characteristic function of the void space. γα Unit weight of an α-phase (= ραg). γαβ Surface tension between α- and β-phases. ΓE γ α Rate of production of Eγα, per unit mass of an α-phase. δ Dirac delta function. δij Components of Kronecker delta. Δ Characteristic distance from solid surface to the fluid in an REV. Δα Hydraulic radius of an α-phase. ε Volumetric strain. Small parameter in perturbation and homogenization. ε Strain tensor (Components εij). ε̇ Rate of strain. εsk Dilatations of the solid matrix (skeleton). ζ Saltwater-freshwater interface location. θ Weighting factor for Crank-Nicolson scheme. θLG Contact angle (between liquid and gas). θα Volumetric fraction of an α-phase (≡ Uα/Uo = φSα). θαr Irreducible, or residual volumetric fraction of an α-phase. λ Coefficient of radioactive decay (= 1/t1/2, where t1/2 = mean half life). Viscosity associated with fluid compressibility. Leakage factor. λ′′s Lamé constant of an elastic solid matrix. μ′s Lamé constant of an elastic solid matrix. μα Dynamic viscosity of an α-phase. μγα Chemical potential of a γ-species of an α-phase. να Kinematic viscosity of an α-phase (≡ μα/ρα). νγ Stoichiometric coefficients of γ species. ξ Location of particle in a Lagrangian system. πα Equivalent pressure head in an α-phase. xx List of Main Symbols ρb Bulk mass density of soil. ρα Mass density of an α-phase. ργα Mass density of a γ-species in an α-phase. σ Standard deviation. σ2 Variance. σ Stress tensor. σ′s Effective stress (≡ σ′s). Σαβ Specific area of Sαβ . τ Shear stress. φ, φeff Porosity and effective porosity, respectively. ϕα Piezometric head in an α-phase, in terms of a reference α-fluid. ϕ Surface potential. Strack’s potential. χ(θ) Function of moisture content. χγ Stoichiometric coefficients of γ species. ψ Suction, or matric suction. ψab, ψaa Spatial (temporal) covariance of a and b, and autocovariance of a. Ψab, Ψaa Spatial (temporal) correlation of a and b, and autocorrelation of a. Ψ Stream function. Ψg, Ψm Gravity and matric potentials, respectively. Ψp, Ψsol Pressure and solute (osmotic) potentials, respectively. Ψsw, ΨT Soil-water and thermal potentials, respectively. Ψtotal Total potential. ωγα Mass fraction of a γ-species in an α-phase. Ω Porous medium domain. ∂Ω Boundary of porous medium domain Ω. Subscripts a Air. c Characteristic value. g Gas. f Fluid. i Intermediate wetting phase. im Immobile phase.  Liquid. m Mobile phase. n Nonwetting phase. N Nonaqueous liquid phase (NAPL). o Organic liquid. Oil. pm Porous medium. s Solid. v Void space. w Water. Wetting phase. List of Main Symbols xxi α α-phase. β β-phase. Also, a symbol for all other phases, except α. Superscripts a Air, as a species. H Heat. T Transpose of a matrix. v Vapor. w Water, as a component. γ γ-component. ε Quantities dependent on the small parameter ε in homogenization. Special symbols (..) Average, volume average, or phase average of (..) ( = 1Uo ∫ Uo(..) dU). Ensemble average. (..) α Intrinsic phase average of (..) (= 1Uoα ∫ Uo(..) dU). G̊ Deviation of G from its intrinsic phase average, G α , over an REV. ︷ ︸ (..) αβ Average of (..) over the Sαβ-surface. (..) Average of (..) over the thickness of thin domain. ∇ ·A Divergence of a vector A ( ≡ div A). ∇A Gradient of a scalar A (≡ gradA). DE(..) Dt Material derivative of (..), as observed by the E-continuum.〈..〉 Temporal or spatial average. [[ A ]]α,β Jump in A across an α-β-interface (≡ A|α −A|β). {A} Activity of a species A. Water Resources Systems 3 happens in the unsaturated zone is essential, especially in the connection with subsurface contamination. The role of the unsaturated zone is clearly depicted in Fig. 1.1.1. A por- tion of the precipitation infiltrates through land surface, and then percolates through the unsaturated zone to an underlying groundwater reservoir. On its way, part of it is taken up by plants and transpires to the atmosphere. An- other part may evaporate back to the atmosphere through the ground surface. Water in the unsaturated zone, often referred to as soil moisture, amounts to only about 1.6% of the generally accessible subsurface water (at less than half a mile deep), on a global basis. This amount is of little significance as a source of water. However, the moisture in the unsaturated zone plays a very important role in providing nutrients for plants, as well as transport- ing, transforming, and temporarily storing dissolved substances originating at ground surface and carried downward by the water. 1.1.2 Surface water versus groundwater In order to discuss the role that groundwater may play in the management of regional water resources, let us assume that both surface water and ground- water are present in relatively significant quantities in a region. Actually, surface water (in lakes and streams) and groundwater (in aqui- fers) are not necessarily independent water resources. Consider, for example, the interrelations between a river (or a lake) and an adjacent aquifer, or a river passing through a region in which a phreatic aquifer exists. If the river (or lake) bed is not completely clogged, water will flow through it from the river into the aquifer when water levels in the former are higher than in the latter, and vice versa. Base flow in streams is nothing more than groundwater emerging at ground surface. In this way, rivers and lakes in direct, continuous hydraulic contact with adjacent, or underlying aquifers serve as boundaries to the flow domain in the latter. By controlling water levels on such boundaries, we can control the flow of water through them into or out of an aquifer. Discharge from a spring is another example of groundwater emerging, un- der certain conditions, at ground surface and becoming surface runoff. By controlling groundwater levels in the vicinity of a spring, its discharge is controlled, or even stopped completely. The above considerations apply not only to water quantity, but also to wa- ter quality, defined, for example, by some chemical species or microorganisms carried with the water. Polluted surface water may easily reach and pollute groundwater and vice versa. Thus, it is obvious that the management of regional water resources should always include both resources, incorporating each of them in the overall sys- tem according to its individual features. In one way or another, any control of one resource will eventually, if not immediately, affect the other. The pos- sible time lag may be due to storage and/or the relatively slow movement of groundwater and of pollutants carried by it. One should note, however, 4 INTRODUCTION that the water divide delineating a groundwater basin and that delineating a surface one are not necessarily geographically identical; in fact, generally, they are not. Depending on the geographical boundaries of any such region, management may include transfer of water from one basin to another within the framework of regional conjunctive use. Although it seems obvious that groundwater, when present in a region, should be used conjunctively with surface water within the framework of any development and management scheme, in many parts of the world, one finds a certain degree of reluctance to include groundwater in such schemes. Perhaps, in part, this attitude stems from the fact that unlike water in streams and lakes, one cannot actually see groundwater in aquifers. In trying to rationalize this attitude, the following reasons were given (Wiener, 1972): • Exploitation of groundwater is energy consuming and expensive, especially when the water table, or the piezometric surface, is deep. • Planning the development of groundwater resources requires long-term data, which, usually, are not available. • Evaluation and planning groundwater resources requires highly trained personnel that are not available. • It is difficult to predict the response of an aquifer (in terms of both water quantity and quality) to proposed management schemes. • Groundwater projects are usually single purpose ones, namely, to supply pumped water, whereas most surface water projects have multiple pur- poses, e.g., to supply water, produce energy, and provide recreation. Obviously, in order to examine these arguments and to compare surface water with groundwater, one must know the local conditions. In general, however, it seems that at least in part, these arguments are based on lack of knowledge. For example: • It is true that when pumping heads are large, energy costs may be sig- nificant (whereas energy may be produced from surface water). However, if one includes in the annual expenditures also the relatively high initial investments required for hydraulic structures, such as storage dams, diver- sion structures, canals, and pipes, the overall economic picture may sway in the direction favoring groundwater. • Because of the large storage and slow motion involved, groundwater levels at any instant reflect the accumulated effect of a rather long period of time; changes are relatively small and slow, in comparison with those of surface water. Hence, in general, shorter groundwater records give suffi- cient information for planning purposes, whereas much longer records are required in order to obtain a complete picture of the more frequent and rapid fluctuating behavior of surface water. • It is true that a certain amount of knowledge is required, but this knowl- edge has been developed to the extent that nowadays this knowledge is included in most training programs of hydrologists and civil and environ- mental engineers, or in special courses of continuing education. Most of Water Resources Systems 5 the necessary theory is also included in the present text, as a contribu- tion to the dissemination of knowledge on groundwater. Consequently, the lack of skilled personnel can easily be overcome, even in regions where this subject has been neglected in the past. • With modern hydrological tools, there is no difficulty in modeling the behavior of a groundwater system and forecasting its response (both quantity- and quality-wise). In general, the forecasts are reliable. In most cases, because of the system’s complexity, digital computers and computer programs have to be used. • Certainly, unlike surface water, groundwater cannot be used for recreation. Nevertheless, and perhaps to a more limited extent, groundwater projects may also serve multiple purposes. For example, in addition to water supply, drainage and reclamation of land may be achieved. We have already men- tioned above that by controlling groundwater levels in adjacent aquifers, we can control base flow in streams, and this, in turn, has water quality as well as ecological consequences. Using the purifying and mixing properties of an aquifer, artificial recharge can be used for the disposal of reclaimed sewage water, thus augmenting groundwater quantity. 1.1.3 Characteristics of groundwater Our main purpose in bringing these arguments is not to show that groundwa- ter utilization is always superior and more advantageous, but to emphasize, again, that whenever the two resources are present, they should be used con- junctively, according to their individual features. Following are some of the main characteristics of groundwater (Wiener, 1972). Location. Surface water flows along fixed curved paths. Their utilization usually requires the construction of regulative facilities that will make the water available only along certain portions of their path. On the other hand, groundwater, when present, underlies extended areas. If these coincide also with demand areas, there is almost no need for a surface distribution sys- tem, as the aquifer acts also as a conduit, and consumers can pump their share directly from the aquifer. This feature is of special interest in regions where development is gradual. More wells are sunk whenever an increase in pumping is required. Often control structures for surface water (e.g., dams, or diversions) cannot be built in stages, and the one-time investment is large. Flow and availability. Fluctuations in surface flow may be significant. Minimum flows, including zero flow, occur often during the season of highest demand for water. On the other hand, climatic fluctuations in groundwater levels are usually small relative to the aquifer’s thickness, so that the large volume of water stored in the aquifer may serve as a buffer that can supply water in periods of drought. Whereas the regulation of surface flow requires hydraulic structures, which are often rather costly, the regulation of ground- water flow is incorporated in the implementation of management schemes, 8 INTRODUCTION management of surface water. Because of the interrelations between surface water and groundwater, and the recommended conjunctive use wherever and whenever called for, these two kinds of water should be incorporated within a single, unified legal and administrative framework. In establishing the legal and administrative framework for the exploitation of groundwater, it is important to consider some of its basic features. The entire aquifer may be regarded as a single basin from the point of view of its water balance. In the long run, all consumers together cannot withdraw more than what is made available by the water balance, which takes into account all inputs (natural and artificial) and all outputs. Temporarily, excess of outflow over inflow can be allowed by reducing storage, i.e., decline of water levels. By pumping, each well produces a drawdown also in its vicinity and may affect the pumping of neighboring wells. The aquifer is also a basin with a certain internal flow pattern established by the pattern of pumping and recharge, both natural and artificial. Pollutants reaching the aquifer will be transported according to this flow pattern, and may reach areas at large distances from where they were originally introduced. In this way, many wells downstream of a source of pollution may be affected. All these factors call for basin-wide regulation of groundwater withdrawal, if aquifer sustainability is set as a goal. Finally, the problems of groundwater quantity and quality cannot be han- dled separately, as they are interrelated. For example, seawater intrusion depends on the rate of freshwater flow from the aquifer to the sea, and the movement of contaminants depends on the flow pattern. All these considerations lead to the conclusion that the management of an aquifer should be centralized and that it requires an appropriate legal and institutional framework. One cannot leave individual land owners to pump according to their needs, or to allow them to dump pollutants on their land. 1.1.4 Functions of aquifers From the discussion presented above, it becomes obvious that beyond serving as just a source for water, an aquifer is a system that should be managed and operated as a unit to achieve various objectives. For example: Source for water. This is the more obvious function. When an aquifer contains only water stored in it from the far past, usually under different cli- matic conditions, this water should be considered a nonrenewable resource. However, in general, an aquifer is replenished annually from precipitation over the region overlying it, or over its intake region (if it is confined). Thus, in general, groundwater is a renewable resource. Obviously only a certain part of the precipitation, depending on the distribution of storms, land to- pography and cover, permeability of soil, etc., infiltrates through the ground surface and replenishes the underlying phreatic aquifer. Aquifers can also be replenished from streams (with permeable beds) and floods. In many arid regions, aquifers in the low lands are replenished during a very short period Water Resources Systems 9 once in several years, from flash floods originating in the mountains. Under natural conditions, a quasi-equilibrium situation is maintained with inflow equal to outflow. Part of the replenishment can be intercepted by pumping, thus reducing the outflow and establishing a new quasi-equilibrium. In this way the aquifer serves as a source for water. Storage reservoir. Every water supply system requires storage, especially when replenishment is intermittent and is subject to random fluctuations. Demand may also fluctuate, e.g., seasonally. A large volume of storage is available in the void space of a phreatic (water table) aquifer. Just to give a rough idea, we can store 15 × 106 m3 of water in a portion of an aquifer of 10 km × 10 km, with a storativity of 15%, raising the water table by 1 m. Using the technique of artificial recharge, large quantities of water can, thus, be stored in a phreatic aquifer. By doing so, water levels rise and outward (from the recharge area) gradients are established. These cause the stored wa- ter to spread over ever increasing areas and/or to leave the aquifer through its boundaries. In this way, if not used, the stored water is gradually lost. Nevertheless, due to the slow movement of groundwater, and with appro- priate management, these losses can be minimized. A possible management procedure is first to lower the water table by pumping in excess of natural replenishment, withdrawing the volume of water as a one time reserve stored in the aquifer between the initial and planned water levels, and then use the dewatered void space for storage. In this way, storage is provided without raising the water table to excessive elevations. Sometimes, we even start by producing a crater in the groundwater table and then filling it up for storage. Again, losses from storage are minimized (recalling always that losses exist also for surface storage by evaporation and seepage). An aquifer can be used for long term storage, e.g., from a wet sequence of years to a dry one, for seasonal storage, from the rainy season to the dry one, or even for shorter periods. The selection of the type of storage, the right combination of surface and underground storage, etc., depends on local conditions, economics, etc. Conduit. Using the techniques of artificial recharge, water can be intro- duced into an aquifer at one point and be withdrawn by pumping at another (or at several other points). The injected water will flow through the aquifer from regions of high water levels, produced by the recharge, to region of low water levels, produced by the pumping. The rate of flow will depend also on the aquifer’s transmissivity. In this way, large distribution systems, say to individual consumers spread out over large areas, may be avoided. Obviously, there is a limit to the permissible rise in water levels, as well as to permissible drawdown; these impose a limit on the use of an aquifer as a conduit. Filter plant. Using the techniques of artificial recharge, an aquifer may serve as a filter and purifier for water of inferior quality injected into it. This may take several forms: 10 INTRODUCTION • By recharging an aquifer (through infiltration ponds) with surface water containing fine suspended load, we remove the fine suspension by the time the water reaches the water table. The bottom of the pond and the soil column act as a filter. • The subsurface in the unsaturated zone and the aquifer material act as huge chemical reactors. Various chemical species may be removed from the flowing water by chemical reactions, by adsorption and by ion exchange phenomena on the solid surface of the porous matrix, especially when clay colloids are present. Obviously, by these reactions, the chemical species are removed from the water, but they remain on the solid. If these chemical species are considered as contaminants, the fact that they remain in the aquifer means that clean-up (remediation) activities will still be required in order to remove them from the aquifer. Certain chemical species will undergo natural attenuation and decay as they move through the aquifer. Often, such attenuation is associated with and enhanced by biological ac- tivities. Of special interest is the reduction in organic matter content as well as the removal of taste, or bacteria and viruses, especially if the flow is of sufficient length and duration. Sometimes, minerals are added to the water by dissolution. • Mixing of injected water with the indigenous water of an aquifer is achieved by their simultaneous movement in the aquifer due to both the mechanism of hydrodynamic dispersion and the geometry of the flow pattern. • Pumping near a river induces recharge from the latter into the aquifer. Fil- tering of the river water and purification are achieved by the flow through the aquifer material from the river to the wells. In each case, the ability of an aquifer to improve the quality of water depends on the chemical and physical properties of the aquifer material and on the type of mineral and organic impurities contained in the water, taking into account that a significant part of the removal of impurities takes place at the phase of entry into the soil material (that is, bottom of an infiltration pond, or vicinity of a recharge well). Control of base flow. This can be achieved in springs and streams by controlling water levels in the aquifer supplying water to them. Water mine. We have already mentioned above the possibility of mining a one time reserve stored in an aquifer between some initial phreatic sur- face and a planned ultimate one. The same is true (using some coefficient of efficiency of mining, due to hydrodynamic dispersion) in the case of the advancement of saltwater in a coastal aquifer toward its planned position. In general, the yield of an aquifer is a long-term average of part of its replenishment (renewable resource). However, under certain conditions, albeit very rarely, we may plan to completely mine an aquifer (like any other non- renewable resource), not worrying about what will happen once this source has been depleted. In this case mining is, usually, based only on economic considerations. Water Resources Systems 13 wastewater by surface irrigation, or dumping of wastewater sludge). Solution mining and in-situ mining activities are also included in this category. Category II: Sources designed to store, treat, and/or dispose of substances, as well as discharge resulting from unplanned release. Here the OTA report lists, primarily, landfills for municipal, and industrial (hazardous and non-hazardous) waste. However, included in this category are also open (legal and illegal) waste dumps, surface impoundments of haz- ardous and non-hazardous liquid wastes, waste tailings, (non-waste) material stockpiles, graveyards and animal burial sites, above-ground storage tanks for hazardous and non-hazardous waste materials, containers for all kinds of non-waste materials, and disposal sites for radioactive materials. Category III: Sources that retain substances during transporta- tion or transmission. Pipelines and material transportation and transfer operations of hazardous and non-hazardous materials, as well as non-waste materials are included in this category. Category IV: Sources that discharge substances to the environ- ment as part of various planned activities. Examples are: irrigation practices, pesticide and fertilizer applications, animal feeding operations, ap- plication of de-icing salts, urban runoff, structures for infiltrating storm water precipitation (carrying atmospheric pollutants), surface and subsurface min- ing, and mine-drainage operations. Category V: Wells and construction excavation. Included here are (oil and gas) production wells, geothermal and heat recovery wells, water supply wells, monitoring wells, and exploration wells. Category VI: Naturally occurring sources whose discharge is cre- ated or exacerbated by human activities. This category includes wa- ter infiltrating from precipitation and carrying atmospheric pollutants, natu- ral leaching, saltwater intrusion, and encroachment of poor quality water as a result of man-made changes in the flow regime in an aquifer. Another way of classifying sources is according to their geometry: Point sources. Here, ‘point’ has to be interpreted as ‘of small areal ex- tent’, relative to the subsurface domain under consideration. For example, a large landfill may be considered as a ‘point source’ relative to an underlying aquifer contaminated by the leachate from the landfill, once the plume of contaminants has reached a distance which is much larger than the dimen- sions of the landfill itself. An oil spill from a ruptured pipeline, leaks from an above-ground, or buried tank, a radioactive waste repository, and a septic tank, may serve as additional examples of point sources. However, if we have a large number of septic tanks within a certain area, we may average (or in- tegrate) their effect, and regard the source created by them as a distributed, rather than point, source. Distributed sources, also called non-point or diffuse sources. Here a source extendsover a large horizontal area relative to that of the contam- 14 INTRODUCTION inated subsurface. The application of pesticides, herbicides, and fertilizers in agriculture are examples. Another example is the infiltration of rainwater carrying atmospheric pollutants (e.g., acid rain). As already emphasized in Subs. 1.1.1, liquid-borne contaminants originat- ing from sources at or near ground surface travel first through the unsaturated zone (≡ vadose zone), (primarily) downward to an underlying aquifer. The rate at which contaminants are transported to the aquifer depends on the quantity of the contaminant, and on the quantity and flux of the percolating water passing through the contaminated soil volume, dissolving and carrying the contaminants. Of special interest is the case of a contaminant that constitutes a non- aqueous liquid phase, i.e., a liquid phase different from water. This case is discussed in detail in Subs. 1.1.5E. In some cases, the moving contaminated liquid, once in the subsurface, may spread out horizontally, beyond the relatively small spreading produced by capillarity, because of anisotropy in soil permeability and/or because of the presence of horizontal lenses and strata of relatively low permeability. Altogether, the issue of whether to regard a source as a ‘point’ source or a ‘non-point’ one depends on the scale of the problem. This scale may evolve in time, such that after a certain period, a plume originating from a diffuse source may be approximated as originating from a point one. Sources may also be classified according to the temporal variation in their rate of release, e.g., a one-time, short duration spill, a slug over some time pe- riod, occurring once or repetitively, or a continuous release over an extended period of time. Some specific sources are described below. C. Point sources Following are examples of some of the more commonly encountered pollution point sources. Septic tanks. These are used as a means of disposal of domestic sewage in many (especially rural) areas. Waste water enters first the septic tank, where scum, grease and heavier than water solids are removed by gravity segregation. The clarified liquid proceeds to the subsurface infiltration sys- tem, where it is discharged to the soil. In general, a properly designed and maintained septic tank should be regarded as an efficient and economical means of domestic sewage disposal. However, even when each unit in itself is well designed and maintained, a high density of septic tanks may exceed the natural ability of the subsurface environment to absorb and purify effluents, thus causing a degradation of groundwater quality due to release of bacteria, organic contaminants and nitrates. Also, because non-domestic sewage is di- verted to them in many cases, and because of unfavorable soil and climatic conditions, septic tanks are considered a potential source for groundwater contamination. Water Resources Systems 15 Raw sewage contains biological contaminants (bacteria and viruses), in- organic contaminants (phosphorus, nitrogen, and metals), and organic ones (synthetic organics, pesticides, and hydrocarbons, which are compounds of hydrogen and carbon). Bacteria tend to be removed from percolating water through the physical process of straining and/or the chemical process of ad- sorption onto soil particles. Phosphorus movement is retarded by chemical transformations and adsorption. Ammoniacal-nitrogen is retarded primarily by adsorption, but can also be subject to cation exchange or incorporation into microbial biomass. Nitrate-nitrogen tends to be highly mobile, moving essentially with the water. Nitrates may also be removed through plant uptake and microbiological denitrification. Metals in soils can be rendered immobile, primarily by adsorption, but also through ion exchange, chemical precipi- tation, or complexation with organic substances. The transport and fate of organic contaminants is affected by volatilization, adsorption, incorporation into plants and microbial biomass, and bacterial degradation. Storage tanks. In terms of the number of incidents (e.g., gasoline tanks and service stations), this is probably the major source of subsurface con- tamination in the USA and other industrialized countries. The tanks may contain hydrocarbons, organic compounds, and inorganic liquid chemicals. The main cause of leakage from steel tanks is corrosion. Gasoline, an example of a hydrocarbon, actually contains dozens of dif- ferent hydrocarbon compounds. It contains aliphatic compounds, such as pentane and butane, which are characterized by carbon atoms linked to an open chain, and aromatic compounds such as benzene and toluene, which are characterized by a ring structure of carbon atoms. Being lighter than water, gasoline, spilled in a sufficiently large quantity, will tend to create a lens which will float and spread out on the water table. Depending on the pressure and temperature prevailing in the subsurface, certain chemical species in gasoline may volatilize, and then diffuse within the gaseous phase. Because diffusion in a gaseous phase is much faster than in a liquid, the vapor of a volatile species may spread extensively. As such vapor spreads in the gaseous phase, it may reenter the water phase, in an attempt to achieve equilibrium be- tween the two phases. Chemical species, such as benzene, toluene, xylene (or BTX), dissolve in water. These species may also be adsorb onto soil particles. Although the solubility level may be low, concentrations above permissible levels are often reached. In many cases, the leakage is not from the tank itself, but from the inlet and outlet pipes and valves at loading and unloading facilities, and from trucks or tank cars while they are being cleaned. Landfills (or sanitary landfills). If properly designed, maintained, and managed (including enforcement of their designation as solely for domestic refuse), landfills should not pose any threat to the environment. In practice, however, they do pose a serious threat to groundwater resources, due to the actual way in which they are operated and the sheer number of such facilities. 18 INTRODUCTION D. Distributed sources Agriculture. Many agricultural activities produce potential sources of groundwater contamination. Among such sources, we may mention pesti- cides and herbicides, fertilizers, animal feed and waste, irrigation, and plant residues. The first three are associated with nitrate (NO−3 ), which is the most common environmental form of nitrogen, because it is the end product of the aerobic biological process called nitrification. Fertilizers (chemical and manure) constitute a serious danger to the sub- surface, both at the handling stage (transportation and storage) and when applied in the field. Irrigation, especially excessive irrigation, may leach and transport significant quantities of nitrogen fertilizers, in the form of nitrate, to underlying groundwater. Animal raising activities produce contaminants that include nitrogen com- pounds, phosphates, chlorides, bacteria, and sometimes heavy metals (Ba, Cd, Cr, Cu, Pb, Ag, and Zn). Irrigation water often contains high levels of total dissolved solids (TDS), especially salts. By taking up water, plants raise the concentration of TDS in the water remaining in the soil within the root zone. As this concentration rises, soil productivity is reduced. When the TDS concentration becomes too high to be flushed by infiltration from precipitation, the salinity can be flushed by excess irrigation. Obviously, the flushed salts end up in the underlying groundwater. Also, as part of the usual irrigation practice, the amount of water applied is often in excess of field capacity (Subs. 6.1.9), so that part of it (called ‘return flow from irrigation’) continues to migrate downward, carrying whatever salinity is contained in the water. Acid precipitation. Acid rain and other airborne contaminants, consti- tute a source for groundwater contamination. The danger is due not only to the mere fact that certain polluting constituents, e.g., sulfates contained in acid rain, are introduced into groundwater; as the acid water travels through the vadose zone, it will increase leaching and mobilization of other chemi- cal species. It is often found that sulfate and aluminum content increases in acidic groundwater. The mobility of heavy metals may also be increased. E. Non-Aqueous Liquid Phase (NAPL) Too often, a contaminant in the form of a third fluid phase—a liquid that is practically immiscible in water—is introduced into the air-water system of the unsaturated zone through ground surface. Organic liquids, e.g., hydro- carbons and various organic solvents, such as TCE (Trichloroethylene), used in industry, may serve as an examples. Another example is pesticides, such as DBCP, which are highly toxic. They are introduced into the unsaturated zone as a result of spills, leaks from faulty storage tanks or pipes, leakage from corroded drums, and burst of pipelines. Referring to water as the aqueous phase, this third fluid is often referred to as a nonaqueous phase liquid, abbreviated, NAPL. It is sometimes referred to Water Resources Systems 19 as free product. In this book, we shall sometimes use ‘oil’ as a typical example of NAPL. The term ‘oil’ will then be used generically, without any implied meaning regarding its actual composition. Thus, when a NAPL is introduced into the unsaturated zone, the latter becomes a three-phase flow domain. Some brief comments on three fluid phases are presented in Sec. 6.5. When its quantity is small, the NAPL will move and spread out in the unsaturated zone, occupying part of the void space within a certain domain, jointly with water and air. If the quantity of NAPL is sufficiently large, this soil domain will continue to expand, primarily downward, until reaching an underlying water table. In the subsurface, certain species of an organic liquid may dissolve in the water moving in the void space (e.g., from infiltration). In this way, these species will reach and contaminate an underlying aquifer. Volatile species may evaporate and spread out by diffusion in the gaseous phase occupying part of the void space. Sharp (microscopic) interfaces are maintained between a NAPL and both the gaseous and the aqueous phases within each pore. Some NAPLs, like gasoline, are less dense than water, while others, such as chlorinated solvents, are denser than water. When the NAPL is less dense than water, it is called ‘light NAPL’, or LNAPL, and, conversely, when it is denser than water, it is called ‘dense NAPL’, or DNAPL. NAPLs may be pure organic compounds, or, like gasoline, complex mixtures of a large number of compounds. The various NAPL components may dissolve in the aqueous phase in small quantities, each according to its own water solubility. Thus, the term ‘immiscible’ is used here in the sense that the water and the NAPL are separated (inside the void space) by a sharp physical interface, despite the small amount of transfer of components between the phases. Unfortunately, some chemical species in drinking water may have deleterious effects on human beings, even at concentrations as low as a few parts per billion (ppb). Because the water solubilities of commonly encountered NAPL components are high, relative to those concentrations (e.g., TCE is soluble to about 1.1 × 106 ppb at 20◦C), the solubility of NAPLs in water is an important factor in such subsurface contamination problems. If a small quantity of NAPL is spilled on ground surface, it will percolate through the vadose zone, primarily downward. As the NAPL travels down- ward, it leaves behind isolated blobs and ganglia of NAPL. If the volume of the spill is not too large, and the water table is at a large depth, after some distance the movement of all the spilled NAPL will be arrested, creat- ing a zone containing immobile NAPL at residual NAPL saturation. When percolating water (from precipitation) passes through this zone, a certain quantity of the NAPL dissolves in the water, and is carried away. In this way, a secondary source of NAPL pollution is created. When the quantity of spilled LNAPL is larger, the downward moving LNAPL may reach the capillary fringe and form a lens floating and mov- ing according to the prevailing hydraulic gradient. The moving lens leaves 20 INTRODUCTION Zone containing free product Zone with NAPL components dissolved in water Vadose zone Groundwater flow Impervious Capillary fringe Infiltration Vadose zone (b) Groundwater flow Impervious Capillary fringe Spill Zone with NAPL components as vapor Zone with NAPL components dissolved in water and as vapor Infiltration Ground surface Ground surface (a) Spill Figure 1.1.2: NAPL spills in large volume: (a) LNAPL, (b) DNAPL. behind a trail of ganglia at residual LNAPL saturation. Eventually the float- ing LNAPL lens may become immobile at residual LNAPL saturation. Water passing through this domain of immobile water will be contaminated. Again, the remaining LNAPL acts as an immobile source for a long time, until all LNAPL has been dissolved. Figure 1.1.2a illustrates an example of a contamination event by LNAPL, e.g., gasoline. From its point of infiltration into the subsurface, the LNAPL flows essentially downward through the unsaturated zone simultaneously with water. Due to heterogeneities, such as layering and the occurrence of low permeability lenses, lateral movement of the NAPL often occurs. Some of the LNAPL will be retained along its flow path by adsorptive and capillary forces at the residual saturation (Subs. 6.2.3). As long as part of the LNAPL is at a saturation that exceeds this critical value, it will continue to move Water Resources Systems 23 Δ Δ Δ Δ I II IIIQ P N R River Aquifer Figure 1.1.3: Yield of an aquifer. or should we adjust it annually, say in response to changes in water table elevations?’. Obviously, to answer this question, we have also to consider (1) the nature of the consumers (e.g., can they tolerate fluctuations in the supply, or is their demand rigid), and (2) is the aquifer (and the consumers) part of a larger system, such that reduction of supply from one source, say the aquifer, can be compensated by increased supply from another source. A starting point for this discussion is the understanding that the yield of an aquifer is, actually, a decision variable, i.e., its value may vary (at least, up to a certain upper limit) according to some management scheme that involves various objectives and constraints. As a simple example, consider steady state flow in an aquifer, isolated on the right side and connected to a river on its left side (Fig. 1.1.3). The aquifer is naturally replenished at an annual rate N , and artificially at a rate R. With no pumping and no artificial recharge of the aquifer, i.e., P = 0 and R = 0, so called ‘virgin conditions’, the discharge, Q, leaving the aquifer to an outlet, say a river, is Q = N , and the water table is in position I (Fig. 1.1.3). As pumping is increased, say to P1, Q decreases to Q = N − P1 and the water table drops to position II. Recall that a water balance must always be satisfied ; in steady state, there is no change in storage in the aquifer. If we further increase P , the water table continues to drop, until, at some value of P = P2 = N , we will have Q = 0. Beyond that P -value, say, P = P3 > N , the water table will drop to below river level, as in position III, and Q will reverse direction. This shows that the extraction rate can exceed natural replenishment, still with a steady state maintained. However, if the water in the river is polluted, or there is a requirement to maintain a certain base flow to the river, these will serve as constraints on the permissible P . Pumping yield can also be constrained by restriction on water table elevations, e.g., due to legal considerations. So, what is the permissible yield? Altogether the allowed maximum abstraction will depend on the constraints imposed on the system, based on various hydrological and socio-economic considerations. Another example is the case 24 INTRODUCTION of the coastal aquifer discussed in Chap. 9, and especially in Subs. 9.2.5A, where the rate of pumping is shown to affect the extent of seawater intrusion, thus making the latter a decision variable. However, constraints need not be related only to water levels. For exam- ple, all aquifers, especially water table aquifers, are continuously threatened by pollution originating from human activities at ground surface, reaching the aquifer after passing through the unsaturated zone. The passage through the unsaturated zone may take some time, often a long time, depending on the depth of the water table and on soil properties. To mitigate groundwater contamination from such sources, we should allow a certain volume of water, Qmin, to continuously leave the aquifer, carrying with it (= flushing) pollu- tants (assuming that there is enough mixing in the aquifer). This means that the yield is also constrained by the need to maintain a certain minimum value of outflow, Qmin, to ensure water quality. Obviously, another option that will achieve the same goal is to undertake a comprehensive approach that protects the land above the aquifer against pollution, as a part of a comprehensive approach to regional management, incorporating the management of land use with that of the underlying water resources. One should note that if we aim at a constant annual yield, the latter is not simply equal to the long term average of the natural replenishment minus the long term average drainage required for flushing purposes. The reason is that it is not a priori obvious that the aquifer can store all excess water of rainy years for use in dry years; this has to be tested and proven, say, by appropriate modeling. It may be interesting to ask: ‘if a decision to increase annual pumping will (eventually) lead to a new lower steady state water table, where did the water that occupied the pore space between the initial water table and the new lower water table go?’ Obviously, it was drained to the outlet (e.g., river, or sea) during the transition from the initial water table to the new one. If there exists no prior obligation/constraint to discharge this water volume to the river, why not pump it (in excess of a yield that is based on natural replenishment alone)? We refer to this scheme as ‘pumping a one-time reserve’. Let us now allow the import of water from some external source (and/or surface water in the same area) and use it for artificially recharging the aquifer at an annual rate R (obviously, if a management decision is made to recharge and then pump, rather than supplying the imported water directly to con- sumers). To maintain a water balance, we can, in principle, increase pumping (above the yield) by the amount of recharge, without lowering water table elevations. On the other hand, if we recharge the aquifer, but we do not in- crease the pumping rate by the same amount, water levels will rise and so will Q. Thus, artificial recharge can be a tool for increasing the yield of an aquifer, only if it is properly incorporated in management considerations. Obviously, the pumping and the recharge (e.g., timing, spatial distribution, and water quality) must be coordinated. Just as an example, if we recharge an aquifer, Water Resources Systems 25 but do not pump in excess of the yield that corresponds to no recharge, then the recharged water will gradually be drained out of the aquifer, and we shall gain nothing by the fact that we have artificially recharged the aquifer. Simi- larly, water quality of pumped water is generally some average of the quality of the recharged water and that of the indigenous water in the aquifer. Obvi- ously, in this preliminary discussion, no attention has been given to economic considerations and constraints, although these may play a significant role. In addition to induced recharge (Sec. 3.5) of polluted water from a river, there may be other sources of inferior water quality (e.g., high salinity) within an aquifer, or along its boundaries with adjacent aquifers, such that this infe- rior quality water may be mobilized by lowering water levels, and eventually pollute the aquifer and the pumped water. Obviously, this may affect the aquifer’s safe yield. So far, we have discussed the relatively simple case of a steady state flow. However, how can we determine the yield in view of the fact that natural replenishment varies from year to year? One way is to take some long term average replenishment as a base for estimating aquifer yield; this is allowed, provided we are sure that the aquifer itself will store the water from rainy years for use in drier years. To examine this point, as part of a management procedure, we run models (Sec. 1.2) in order to discover the highest value possible for a constant aquifer yield, say, over a design period of 30 years. Or, we can consider a yield that varies, say, from year to year, as a function of the volume of water stored in the aquifer (as manifested by water table elevations), and run the model for the same period. Up to this point, the discussion has been based on the assumption that the aquifer should be preserved as a source of water for any foreseeable future. This underlying guideline should be compared, using socio-economic criteria, with (the opposite) approach that regards the aquifer as a mine of water, such that water, similar to oil, can be extracted until the aquifer is ‘empty’, or destroyed, say, by the invasion by water of inferior quality, e.g., seawater in a coastal aquifer. At this stage, let us define ‘sustainable yield’ of an aquifer (again, in fact, of any water resources system) as the maximum volume of water that can be annually withdrawn from the aquifer, for an indefinite period, without causing any undesired effects. In the past, hydrologists used the term ‘safe yield’ to indicate that the yield should be determined such that the aquifer as a source of water will be sustained forever. Its determination is based on regional water balances for the considered aquifer, using as input from precipitation a synthetic sequence (or synthetic sequences of various probability of occurrence) based on the statistics of precipitation. A design period of 30 years is usually used. To emphasize the underlying (and indisputable) requirement that an aquifer should serve as a source of water that must be sustained for an indefi- 28 INTRODUCTION • Return flow is a term used for that part of the water used for irrigation that is not consumed by the vegetation, but infiltrates and reaches the water table. Again, we can relate return flow to the (measurable) volume used for irrigation through certain parameters and obtain estimates of these parameters in the process of model calibration. • Artificial recharge, inflow from lakes and streams and spring discharge re- quire no further explanation. Artificial recharge is a measured quantity. Inflow from lakes and streams has to be estimated as function of the dif- ference in water levels, and coefficients (= calibration parameters) that express soil properties and geometry. • Evapotranspiration is the term used for water lost to the atmosphere by evaporation from the water table and transpiration through the vegetation. This loss, which is of significance only in the case of a shallow water table, has to be estimated. • As a consequence of excess of all inflows over outflows, the water table (or the piezometric head) will rise, thus storing this excess. Here, the relevant coefficient is aquifer storativity (S or Sy), again an unknown coefficient that has to be estimated in the process of model calibration. Altogether, on the basis of water level observations during a sufficiently large number of years, we obtain a set of annual balance equations and use them to derive the unknown parameters and coefficients for the considered aquifer. Although there exist a number of software packages that help in per- forming this process of parameter estimation, often, in practice, calibration is performed by a process of trial and error, in which, basically, we assume a set of coefficients and examine whether, by using them in a model, we can reconstruct (‘predict’ the observed past) the measured water levels. The main outcome of such process is information on the various parameters that enable us to estimate the replenishment and its relationship to precipi- tation (on which much more data are available), and aquifer coefficients that are required for modeling the flow regime in the aquifer. Once we have a calibrated model of the considered aquifer, we can run models for various values of aquifer yield, predict water levels (which enable us to estimate groundwater flow and storage changes) and examine whether or not that yield causes the violation of imposed constraints. In this way, we can derive (or estimate) the sustainable yield of the aquifer. If quality of water is a constraint, we also have to run calibrated solute transport models. Altogether, by using appropriate flow and solute transport models, we can also derive the sustainable yield as a variable that depends on the varying elevations of the water table in the aquifer. We can also investigate artificial recharge (using imported water) as a tool for increasing the sustainable yield. Modeling 29 1.2 Modeling The management of a groundwater resources system, which can be an aquifer, or a system of aquifers, alone or conjunctively with surface water sources, aims at achieving certain goals through a set of decisions concerning the de- velopment and/or operation of the system. A more detailed discussion on the management of groundwater resources, particularly the quantitative analysis involved, is presented in Chap. 11. The remarks in the present section are aimed at providing the background needed for the discussion on modeling of the physical and chemical behavior of water in aquifers, presented through- out this book. Particularly, we wish to emphasize the need for employing the various models introduced in this book as essential management tools. The modeling concept and process are elaborated in the following sections. 1.2.1 Modeling concepts A model may be defined as a selected simplified version of a real system and phenomena that take place within it, which approximately simulates the system’s excitation-response relationships that are of interest. For example, a groundwater system may be ‘excited’ by pumping, by the introduction of a contaminant, by artificial recharge, or by changing a boundary condition. Its ‘response’ takes the form of spatial and temporal changes in water levels and in contaminant concentrations. With the above definition, we emphasize that all that a model can do is to predict the future behavior of an investigated, say aquifer, system. However, this information may be used in different contexts and for different purposes. Modeling activities may be conducted to achieve any of the following objec- tives: (a) To predict the behavior of a system, say, an aquifer, in response to exci- tations that stem from the implementation of management decisions. (b) To obtain a better understanding of a system from the geological, hydro- logical, and chemical points of view. (c) To provide information required in order to comply with regulations. (d) To provide information for the design of a monitoring network, by pre- dicting the system’s future behavior. (e) To provide information for the design of field experiments. In principle, for both (a) and (b), we need to write a flow and, sometimes, also a solute transport model, but the requirements for details and accuracy may be significantly different. This difference may strongly affect the amount of (geological, hydrological and geochemical) data needed and its accuracy requirements. Figure 1.2.1 shows how modeling activities are incorporated in the man- agement process of a given management problem. We always start from the given management problem, identifying what we wish to achieve in order to 30 INTRODUCTION Management problem  Objectives  Propose alternative solutions  Use models to predict system’s response to alternatives  Evaluate objective function and constraints   Select preferred alternative  Implement alternative  Monitor  Figure 1.2.1: The role of modeling in the decision-making process. solve it. We express what we wish to achieve in terms of an objective function (Sec. 11.1). For example, suppose we wish to cleanup a contaminated domain within a certain time interval at minimum costs. The objective function ex- presses the cost of each alternative solution, and we wish to minimize this cost. In this example, the prescribed time interval is a constraint that the solution must satisfy. Obviously, the cost is associated with the behavior of the system. We then identify a number of alternative solutions, say, various well locations and pumping rates, with the intention of selecting the one that will be the ‘best’ (or the ‘optimal’) from the point of view of the selected objective function. We determine the information that is required in order to enable this selection, e.g., information on how will the system respond to the implementation of the various management alternatives. The ‘response’ is then expressed in terms of ‘cost’. The model is the tool that will provide the information on the system’s response. We then construct the model (and that includes model validation, calibration, etc., to be discussed below), and use it to obtain the required information. The latter is then used to select the most desired alternative. The subject of aquifer management is discussed in detail in Chap. 11. As a solution is being implemented, additional data is continuously being assembled and used to reevaluate the problem, to update Modeling 33 • The domain’s hydrogeology, stratigraphy, etc. • The dimensionality of the model (one, two, or three dimensions), and the geometry of the boundary of the domain of interest. • The behavior of the system: steady state or time-dependent. • The kind of soil and rock materials comprising the domain, as well as inhomogeneity, anisotropy, and deformability of these materials. • The number and kinds of fluid phases (e.g., water, air, NAPL), and the relevant chemical species. • The extensive quantities transported within the domain. • The relevant material properties of the fluid phases (density, viscosity, compressibility, presence of solutes). • The relevant transport mechanisms within the domain. • The possibility of phase change and exchange of chemical species between adjacent phases. • The relevant chemical, physical, and biological processes that take place in the domain. • The flow regimes of the involved fluids (e.g., laminar or non-laminar). • The existence of nonisothermal conditions and their influence on fluid and solid properties and on chemical–biological processes. • The presence of assumed sharp macroscopic fluid-fluid boundaries, such as a phreatic surface. • The relevant state variables, and the areas or volumes over which averages of such variables should be taken. • The presence of sources and sinks of fluids and contaminants within the domain, and their nature (spatial distribution and temporal variation). • The initial conditions within the domain, and conditions on its boundaries. Obviously, more items may be included in the conceptual model of any spe- cific case. The set of assumptions comprising the conceptual model is expressed in words. It is recommended that these assumptions be numbered, say [A1], [A2], etc., as is done for equations, so that they can be referred to. The selection of the appropriate model for a particular case depends on three main factors: • The objective(s) of the investigations, i.e., what kind of information is the model required to provide for the purpose of making management deci- sions. Here we may include rough preliminary estimates vs. more detailed predictions. • The available resources required to construct and solve the model. Included here are the availability of expertise, skilled personnel, and computers. Also included is the ability to describe processes that take place, and the availability of field data required to validate the model and determine the numerical values of its coefficients. • The legal and regulatory framework which pertains to the considered case. 34 INTRODUCTION The objectives of an investigation dictate which features of the domain and its behavior should be represented in the model, and to what degree of accuracy and detail. The constraints imposed by limited resources are very real and, although we shall not dwell on them in this book, cannot be overlooked. For example, a more detailed model is, obviously, more costly and requires more skilled modelers, more sophisticated computer codes, and larger computers. As we shall emphasize below, it is important to understand that a more detailed model requires more field measurements in order to calibrate it and to de- termine its coefficients. Data acquisition is usually more expensive than code and computer costs. Selecting the appropriate conceptual model for a given problem is the most important step in the modeling process. If we oversimplify, we may not pro- duce the required information. If we undersimplify, we may have neither the information required for model calibration (see below) and coefficient deter- mination, nor the resources to solve it. If we select inappropriate or wrong assumptions, our model may not represent those features of the system be- havior that are relevant to the management problem on hand. The set of assumptions serves as a ‘prescription label’ of the model. One should not use a model developed for a different problem, unless one examines its ‘label’ to ensure that it fits one’s problem. Step 3: Development of a mathematical model. In this step, the conceptual model is expressed in the form of a mathematical model. The continuum type of model, discussed in Sec. 1.3, is usually (but not always) employed. The mathematical model consists of: • A definition of the geometry of the surfaces that bound the domain. • Equations that express the balances of the relevant extensive quantities (e.g., mass of fluids, mass of chemical species, energy). • Flux equations that relate the fluxes of the extensive quantities to the relevant state variables of the problem (e.g., Fick’s law for the diffusive mass flux of a chemical species in a fluid phase). • Constitutive equations that define the behavior of the particular phases and chemical species involved (e.g., dependence of density and viscosity on pressure, temperature, and solute concentration). • Sources and sinks, often referred to as forcing functions, of the relevant extensive quantities. • Initial conditions that describe the known state of the system at some initial time. • Boundary conditions that describe the interaction of the domain with its environment (i.e., outside the delineated domain) across their common boundaries. Mathematical models of flow nd solute transport are discussed in Chaps. 5 and 7, respectively. Modeling 35 In a continuum model (Sec. 1.3), a partial differential balance equation describes the behavior at every point within the domain. However, sometimes we are not interested in what happens at every point. Instead, we need in- formation on the lumped, or averaged, behavior of an entire domain, or of parts of it. Another kind of model is called for in such cases, referred to as a lumped parameter model, a compartmental model, a multi-cell model, or an input-output one (Subs. 7.2.1). In such a model, balances of the relevant extensive quantities are stated for ‘cells’ of different shapes and sizes of the domain; state variables are averages over these cells. Sometimes, the hetero- geneity of a domain is such that an appropriate representative elementary volume (Sec. 1.3) cannot be found for it and the continuum approach cannot be applied. A lumped parameter model may be required. In this book we focus our attention mainly on continuum models. The core of any model of transport on an extensive quantity is the balance equation of that quantity. In this book, we deal with models of transport of two extensive quantities: mass of a fluid phase, and mass of a component of a fluid phase. The extensive quantity of momentum of a phase is not referred to directly, as we replace it by an approximation that takes the form of Darcy’s law. We do, however, refer to the momentum of a solid phase when we deal with porous medium deformation. In passing from a model at the microscopic scale to the macroscopic scale model (these terms are defined in Subs. 1.3.1), by some process of averaging, various coefficients of transport, transformation, and storage of the extensive quantities are introduced. The permeability of a porous medium (Subs. 4.1.3), moisture diffusivity (Subs. 6.2.2), and dispersivity (Subs. 7.1.6) are examples of such coefficients. Permeability and dispersivity of a phase are examples of coefficients that express the macroscopic effects of the microscopic config- uration of the boundaries between a considered phase and all other phases present in a representative elementary volume (REV) of the medium. It is important to realize that the coefficients derived by field experiments, in which (in principle) we compare measurements of certain state variables with the corresponding values predicted by a model, actually correspond to that particular model. If possible, one should not employ coefficients derived for one model in another one. When the coefficients developed for one model are used in another model, errors may result. The magnitude of the error will depend on the differences between the two models. In principle, to employ a particular model for a particular domain, the values of the coefficients that appear in the model should be determined from field experiments conducted in the domain. Typically, such an experiment involves a comparison between actual field observations and predictions made by the model, employing some parameter identification technique (Sec. 11.3). Step 4: Development of a numerical model and code. Having con- structed a mathematical model, in terms of relevant state variables, it has to be solved for cases of interest. The preferred method of solution is the 38 INTRODUCTION ill-posed, and there is no reason to expect a unique solution. The inverse problem is discussed in Sec. 11.3. When conditions for determining coefficients, as described above, do not exist, they can be created as a laboratory experiment, usually on a soil core taken from the considered formation, or as a field test. In such experiments, we create a situation (in the laboratory, or in the field), in which state variables can be observed and for which a solution (preferably, analytical) of the model can be determined. It is important to emphasize again that coefficients for any particular site must be determined by making use of data assembled from that particular location. In addition to model calibration, or solving the inverse problem for the considered domain, model coefficients and parameters can also be obtained (better: estimated) from: • Literature survey. In many cases, we can find in the literature values of coefficients and parameters that have been derived and used in modeling similar situations, similar soils, fluid phases and chemical species, etc. One should be careful in using this information, at least as far as soils are concerned, as soils at different sites, even when belonging to the same class in some standard classification, seldom behave identically. Nevertheless, the values found in the literature may be employed as first estimates in a model calibration process, or in sensitivity analysis runs. Such runs are conducted in the process of planning field experiments aimed at obtaining the site specific values for such coefficients. • Laboratory experiments. Laboratory experiments provide good in- formation on parameters that are independent of soil characteristics. As for soil dependent coefficients, we should recall that our modeling is at the field scale, which always involves heterogeneities, etc. Laboratory experi- ments are usually carried out on relatively small samples (‘cores’), and not always under undisturbed conditions. Hence, at best, the values obtained in such experiments should be considered as first estimates to be used in model calibration runs. • Small scale field experiments. These include both standard tests, like pumping tests, slug tests, etc., as well as specially designed ones. If properly designed, they could serve for model validation as well as provide site-specific information on values of coefficients. Step 8: Model application. Once a model has been calibrated for a considered problem (and this includes all the required site-specific coeffi- cients and parameters), the model is ready for use. Computer runs are then conducted to provide the required forecasts. Step 9: Sensitivity analysis. The term sensitivity analysis is used here to describe tools that help the modeler evaluate the impact of uncertainty, say, in the values of model coefficients, on the results predicted by the model. Briefly, we want to know how sensitive are the predicted values to changes in the values of model coefficients. If these effects are not significant (from the Modeling 39 point of view of the decision maker, who makes use of these predictions in the management process), we can accept the predicted values and make decisions. If, however, the predicted values are sensitive to changes in parameter values, we must reduce the range of uncertainty in the values of these parameters. In most cases, this means that we must invest more resources in order to acquire more and more accurate data. Sensitivity analysis can also be used to assess the reliability of parameters determined in the calibration, or parameter estimation procedure described in Step 7 above. Typically, a residual error is usually expressed in the cali- bration process, e.g., as the sum of the squared differences between the mea- sured and the predicted water levels. The optimal set of model parameters is the one that minimizes this error. This subject is discussed in Chap. 10. There, we show also how to enable a more accurate determination of model parameters and how to asses the reliability of the obtained parameter val- ues. In Subs. 11.3.1, we also present the sensitivity analysis in connection with the analysis of pumping tests used for determining aquifer parameters. Altogether, if a small change in a parameter causes a large change in the residual error, we may say that the residual error is sensitive to that param- eter. Hence, our information on that parameter is of high quality, or reliable. Conversely, if a large change in the parameter causes only a small change in the residual error, we say that the calibration process is insensitive to that parameter; hence, our information on that parameter is of low quality, or less reliable. Poor quality information means that more data is required in order to improve calibration results, and vice versa. When the quantity of available data (e.g, water levels at observation wells) is insufficient for performing regional model calibration, we use information on the domains stratigraphy, which should be available, to estimate values of model coefficients. Then, it is certainly of interest to study the sensitivity of the model’s predictions to the uncertainty in these input parameters. More information on how such analysis can be performed is presented in Chap. 10. Again, a high sensitivity indicates that it may be worthwhile spending extra efforts and resources, such as commissioning a geological or geophysical sur- vey, in order to improve our information concerning that feature. Sensitivity analysis is discussed in more details in Subs. 10.2.2. Altogether, the result of a sensitivity analysis can either increase the mod- eler’s confidence in the model, or, conversely, reveal the deficiency of the model. In the latter case, the sensitivity analysis can be used to help the manager to efficiently allocate resources to expand data acquisition, thus improving model prediction. Step 10: Stochastic analysis. The sensitivity analysis discussed in Step 9 provides a qualitative description of uncertainty. It addresses questions like ‘what if there exists 20% uncertainty in this parameter or that condition’; it does not express the range of uncertainty in either the input or the output in terms of statistical measures, such as mean and standard deviation. 40 INTRODUCTION A stochastic analysis not only takes into account the simple statistical measures of mean and standard deviation of the input data, such as the hydraulic conductivity and the natural replenishment, but it also examines the temporal and spatial correlations of these data. Following are some typical questions that can be asked : • What is the predicted mean piezometric head at a given time and location? • What is the standard deviation of that prediction? • If a 90% reliability is needed in the manager’s recommendation, what is the allowable pumping rate? • If a head is observed at a certain location at a certain time, what is the probability that the value will stay correlated a certain distance away or after a certain time period? In Chap. 10 we shall present and discuss the stochastic analysis technique that can provide answers to the above typical questions. Step 11: Summary, conclusions, and reporting. The summary and conclusions should include the information that the model was expected to provide, including additional information concerning the accuracy of the in- formation, the uncertainty involved, and suggested follow-up work. The re- port on the modeling activities may be part of the report on solving the management problem, say as an appendix, or as a separate report. 1.2.3 Model use Admittedly, the subsurface is highly heterogeneous, and we seldom have enough data to calibrate the model of an aquifer domain in a way that will accurately describe its heterogeneity. We are also uncertain about model boundaries and conditions occurring on them, and even less certain about what will occur in the future. In contamination problems, we never have sufficient knowledge and data concerning all the chemicals and biological transformations that may take place. Usually, we have insufficient knowl- edge of the thermodynamic information required for modeling the complex case of multiple multicomponent phases, which may be under nonisothermal conditions. What, then, is the use of a model? Can the model be expected to predict future states to any desired, or acceptable, degree of accuracy? In some cases, such as in the case of a repository for radioactive waste located at depth in the unsaturated zone of a highly heterogeneous fractured rock, the required extrapolation in time is for hundreds and thousands of years, and the domain of interest extends for many kilometers. Can we validate a model under lab- oratory, or small scale field conditions and then apply it to such a large scale and long term problem? The answer is, generally, not. An argument for the use of models is that, in most cases, decisions will be made anyway, with or without sufficient information. The use of a model can, at least, provide some information, even if it lacks the desirable accuracy. Continuum Approach 43 croscopic representative elementary volume (microREV), denoted as μREV, around every point within the domain occupied by the phase. Thus, a ‘fluid phase’ and a ‘solid phase’ are continuum concepts, obtained from the molec- ular level by volume averaging over appropriate representative elementary volumes. This idea will be further explained below. A phase may be composed of a large number of different chemical species. A chemical species is an atom, a molecule or an ion distinguishable from the rest of the phase due to its chemical composition. Often, the number of species in a phase is very large and in a liquid they may also interact chemically. However, under conditions of chemical equilibrium, the minimum number of independent chemical species necessary to completely describe the composition of a given phase may be much smaller. We use the term com- ponent to denote a chemical species that belongs to the smallest set of such species that is required in order to completely define the chemical composi- tion of a phase under equilibrium conditions. As a simple example, consider liquid water in contact with its vapor. The water is composed of oxygen and hydrogen, with these two elements being always present in fixed and defi- nite proportions. Therefore, the system is composed of two chemical species, but only a single component. When chemical equilibrium is not assumed, all species are defined to be components. In this book, we use the terms chemical species and component interchangeably, unless we wish to emphasize that we are referring specifically only to one of them. Often, only a subset of the species within a phase is in equilibrium. Then, the set of components consists of the components of the system that are in equilibrium, together with the species that are not in equilibrium. For the sake of simplicity, we shall often use the term ‘component’ also to denote the mixture of a number of independent chemical species in a liquid or a gas. The selection of components is not unique, in the sense that different chemical species, may be selected as components of a given phase. 1.3.2 Need for continuum approach A spatial domain is said to behave as a continuum if state variables and properties that describe the behavior of the material occupying it can be assigned to every point within it. With this definition in mind, if we consider a porous medium domain, containing a solid matrix and at least one fluid phase, we can identify within this domain the subdomains occupied, respectively, by the solid and by each fluid phase. The domain occupied by water behaves as a continuum with respect to water, and so is the domain occupied by solid. As defined earlier, we refer to the level of description of phenomena within each phase as a description at the microscopic level. As an example, let us consider a porous medium domain whose void space is entirely occupied by water. The domain occupied by solid behaves as a continuum with respect to the solid. For example, we can define solid density at every point within the solid subdomain. Similarly, we can describe certain 44 INTRODUCTION properties or behavior variables at every point in the domain occupied by wa- ter, e.g., water density, pressure, and velocity. As described in the preceding section, at this microscopic level of description, the water phase is regarded as a continuum within the domain it occupies, and so is the solid. For example, to obtain the velocity distribution within the water, regarded as a Newtonian fluid, we could make use of the Navier-Stokes equation (which is a PDE that describes the linear momentum balance of a Newtonian fluid), and solve it within the water occupied domain, subject to boundary conditions on the solid-liquid interface that bounds this domain. Obviously, for flow through porous medium problems, this approach is, usually, impractical, due to our inability to describe the complex configuration of the solid-liquid boundary over large porous medium domains. Moreover, even if we could solve for val- ues of state variables, e.g., pressure, it would be impractical to verify the solution by measurements at this level. To circumvent these difficulties, associated with trying to solve problems at the microscopic level, another level of description is introduced, referred to as the macroscopic level. At this level, properties are defined at every point in the porous medium domain, such that the knowledge of the complex interphase geometry is not needed. The continuum approach, described in the preceding subsection for the transport of a species within a phase, is obtained by averaging the behavior of the phase at the molecular level (Bear, 1972, p. 17). We referred to this description as one at the microscopic level. This approach of averaging can be extended to a multiphase system such as a porous medium domain, where the various phases are separated from each other by abrupt interfaces. To achieve this goal, the real porous medium domain, containing two, or more phases (each of which is already regarded as a continuum at the microscopic level that occupies a certain portion of space that together completely occupy the porous medium domain), is replaced by a model, in which each of the phases is assumed to behave as a continuum that fills up the entire porous medium domain. We then speak of ‘overlapping continua’, each continuum corresponding to one of the phases present in the domain. If the individual phases interact with each other, so will these continua. The space occupied by these overlapping continua will be referred to as the macroscopic space. For each point within this macroscopic space, average values of the variables that describe the behavior of a phase are taken over elementary volumes, centered at the point, regardless of whether, in the real domain, this point falls within that phase, or outside it. The averaged values are referred to as macroscopic values of the considered variables. By traversing the entire porous medium domain with a moving elementary volumes, thus assigning averaged values to every point, we obtain fields of macroscopic variables, which are differentiable functions of the spatial coordinates. In this way, we have turned the porous medium domain into a continuum; the behavior of each of the phases within this continuum is described at every point by the averaged values of variables and material properties. Continuum Approach 45 The advantages of the continuum model of a porous medium are: • It circumvents the need to specify the exact configuration of the interphase boundaries, the knowledge of which is, anyway, not available. • It describes processes occurring in porous media in terms of differentiable quantities, thus enabling the solution of problems by employing methods of mathematical analysis. • The (macroscopic) quantities mentioned above are measurable, and can, therefore, be useful in solving field problems of practical interest. These advantages are at the expense of the loss of information concerning (1) the microscopic distributions of variables within each phase, and (2) the microscopic configuration of the interphase boundaries. However, the macro- scopic effects of this configuration are retained in the form of coefficients that are created in the process of averaging. The structure and relationship of these coefficients to the statistical properties of the void space (or phase) configuration within the elementary volumes can be determined. Examples of resulting coefficients are the porosity, permeability, and dispersivity. For a specific porous medium, the numerical values of these coefficients must be determined experimentally, in the laboratory, or in the field. In the following subsection, we shall discuss the procedure for passing from the microscopic level to the macroscopic, or continuum, one, by averaging over an REV. In Subs. 1.3.5, we shall introduce another approach for smoothing out heterogeneities. 1.3.3 Representative elementary volume and averages A basic feature, which is common to the porous materials that occupy the subsurface, in fact, common to all porous media, is that both the solid matrix and the void space are distributed throughout the porous medium domain. This implies that samples of a sufficiently large volume, taken at different locations within the domain, will always contain both a solid phase and a void space. How large should such volume be? A porosity (= ratio between the volume of void space to the sample’s volume) can be defined for the sample. Obviously, there is no meaning to porosity at the microscopic level of description. For the time being, let us refer to the volume of such a sample as a representative elementary volume (abbreviated REV). We shall discuss the size of an REV in Subs. 1.3.4A. A. Definition of representative elementary volume Based on this brief discussion, we may now define a porous medium as a portion of space (1) that is occupied by a number of phases, at least one of which is a solid, and (2) for which an REV can be found. This means that a porous medium is not merely ‘a domain containing a void space and a solid phase’. Our definition implies that if a common REV cannot be found for all points of a domain, that domain cannot qualify as a porous medium. 48 INTRODUCTION where θα = Uoα Uo (1.3.6) is the volumetric fraction of the α-phase within Uo. The kind of average to be used in each case depends on the way the averaged quantity is actually measured in the field. For example, if, at a point, we take a liquid sample out of a porous medium domain, in order to determine the concentration of a solute in it, the latter is an intrinsic phase average, as it is taken only over the liquid phase. The measuring device (e.g., screened portion of a piezometer, or porous cup of a tensiometer (Fig. 6.1.9)) of an instrument designed to measure an averaged, macroscopic quantity, must also be of the size of an REV, in order to yield observations compatible with the averaged values calculated by a macroscopic model for that point. Throughout this book, it is assumed that the porous medium comprising the subsurface may be considered as a continuum in the sense explained above. Accordingly, the phenomena of fluid mass transport, solute transport, and heat transport (not considered in his book) are described (modeled) at the macroscopic level. This is the level at which engineers and hydrologists make predictions and measure state variables in the field. However, certain phenomena (e.g., capillary pressure and dispersion of a solute) cannot be understood unless we consider and understand them first at the microscopic level. Once these phenomena are understood and described at the microscopic level, they are averaged to yield their macroscopic description. 1.3.4 Scale of heterogeneity in continuum models A. Size of representative elementary volume We must still select the appropriate size and shape of the averaging volume, earlier referred to as REV. Usually, this size is selected such that (1) the average value of any geometrical characteristic of the microstructure of the void space, at any point in a porous medium, will be a unique function (or almost so, within an acceptable error) of the location of that point only, and (2) the measured averaged value should be independent of small perturbations in the size of the REV. This means that the average value should remain more or less constant over a range of REV volumes that correspond to the range of variation in the size of the sample or instrument that measures that average. Denoting the characteristic dimension of an REV by , e.g., diameter of a spherical REV, and the length characterizing the microscopic structure (or, heterogeneity) of the void space by d (say, the typical size of grain or pore, or the hydraulic radius, which is proportional to the reciprocal of the specific surface area of the solids within an REV), a necessary condition for obtaining non-random estimates of the geometrical characteristics of the void space at any point within a porous medium domain is  d. (1.3.7) Continuum Approach 49 Another condition that sets an upper limit to the size of the REV is  max, (1.3.8) where max is the distance beyond which the spatial distribution of the rele- vant macroscopic coefficients that characterize the configuration of the void space (e.g., porosity, permeability) deviates from the linear one by not more than some acceptable value (Bear and Bachmat, 1990, p. 22). The selection of the size of the REV is also constrained by the requirement that  L, (1.3.9) where L is a characteristic length of the porous medium domain. We shall dis- cuss these limits in greater details in Subs. 1.3.4B, together with the general question of scales and their corresponding elementary volumes. Although we have shown here how the macroscopic level of description is obtained from the microscopic one by volume averaging, other techniques that lead to macroscopic flow and solute transport models are also presented in the literature. Among them, we may mention another volume averaging approach proposed by Whitaker and coworkers (e.g., Hassanizadeh, 1986; Whitaker, 1967, 1986a, 1986b; Plumb and Whitaker, 1990) and the homoge- nization technique (e.g., Sanchez-Palencia, 1980; Hornung, 1997), which aims at smoothing out the heterogeneity at the pore scale, as well as at other scales. The homogenization technique is discussed in Subs. 1.3.5. In the previous section, we have discussed two levels of description of phe- nomena: the microscopic level, obtained by averaging over a μREV (read: microREV), to smooth out spatial variations at the molecular level, and the macroscopic level, obtained by averaging over an REV, to smooth out varia- tions at the microscopic level, resulting from the presence of a solid matrix and a void space within a porous medium domain. With the above in mind, we can extend this smoothing out approach to other scales of heterogeneity. Altogether, we shall have: • molecular scale, when we consider the behavior of individual molecules (obviously, impractical). • microscopic scale, at which we describe what happens at points within every phase present in the void space. • macroscopic scale, which is the usual scale of describing phenomena in porous media. • megascopic scale, which is the scale at which we describe phenomena of transport in the field. The latter is, usually, very large and very heteroge- neous in a random fashion. By smoothing out heterogeneity at one level we obtain a model at the higher level. The smoothing operation introduces ‘coefficients’ that represent the effect of the smoothed out properties. Note that the spatial distribution of the coefficients at the higher level may still be heterogeneous. Heterogeneity 50 INTRODUCTION of a continuum domain with respect to a given property at the macroscopic level will be defined in Sec. 2.5. B. Averaging over microscopic heterogeneity For the purpose of the discussion here, let us refer to a given domain (re- garded as a continuum) as homogeneous (= opposite of heterogeneous, or inhomogeneous) with respect to a given property, if that property has the same value at all points of the domain. For example, if the property of inter- est is ‘the presence of solid’ at a point, then, by the very definition of a porous medium (Subs. 1.3.3A), this domain is always heterogeneous with respect to that property. As an example, let us compare two porous medium domains, one filled with large gravel, say, 1 cm in mean radius, and the other with sand, say, 1 mm in radius. Both contain a solid matrix and void space, and both are heterogeneous at the microscopic level. If the granules are spherical and in cubical arrangement, both will have the same porosity. What is the difference in heterogeneity between these two porous media? If we pick a point inside a solid grain, we have a good chance of finding solid at points within a distance of 1 cm in the first case, and within a distance of 1 mm in the second. We can say that their ‘scales of heterogeneity’ are 1 cm and 1 mm, respectively. Extending the above example to a porous medium with random grain arrangement, we can define a distance along which a selected property (here the ‘presence of solid’) at a point is strongly correlated to that at another point. We continue to discuss what happens when we consider flow and transport in large field domains, say, in an aquifer or in the unsaturated zone, which are always highly heterogeneous with respect to (macroscopic) flow and so- lute transport coefficients. For example, we consider the permeability (in a macroscopic level continuum), and ask: ‘given the permeability at a point A, how far from A is the permeability still correlated to that at A?’ This ques- tion makes sense, because the heterogeneity in permeability in the field is a consequence of geological processes that produced it in the first place. Re- peating this question for a large number of points at various distances around A, and for various points A within the considered domain, we shall find a certain distance that we refer to as the scale of heterogeneity of the given domain, with respect to the considered property, here, permeability; within that distance the permeability is still strongly (or sufficiently) correlated to that at A. It is natural to assume that the scale of heterogeneity introduced above is related to the size of the considered domain, with larger domains exhibiting higher heterogeneity; so when we consider domains with characteristic lengths of tens, hundreds and thousands of meters, the scale of heterogeneity will also gradually grow. However, we usually assume that at some field size, this scale will level off at some value. Continuum Approach 53 in which the behavior of each phase and component can be described as a continuum at the macroscopic level. As we shall see throughout this book, the core of any model that describes the transport of an extensive property (e.g., mass, mass of a component, en- ergy and momentum) in a considered domain, is the balance of that quantity at a point, meaning within a small volume centered at the point. When writ- ten per unit volume of the domain, this balance takes the form of a partial differential equation (PDE). At the microscopic level, i.e., within a phase, this volume is the μREV, and the PDE describes the (instantaneous) bal- ance of that quantity as it accumulates within the μREV. This accumulation is a consequence of the net inflow of the considered quantity (by advection and diffusion) and the rate of its production, per unit volume of the μREV centered at the point. Since averaging is actually an integration operation (dividing the result by the volume over which the integration has been per- formed), and averaging a sum is equal to the sum of averages, the averaged equation is obtained by averaging each term in the equation over the μREV centered at the point. Once we have the balance equation at the microscopic level, we can average it to obtain its macroscopic counterpart. Bear and Bachmat (1990) describe the methodology of averaging a partial differential equation. An important example of deriving an averaged equation is the derivation of the averaged momentum balance equation for a fluid that moves in the void space of a porous medium. The movement of a Newtonian fluid is described by the Navier-Stokes equation (see most textbook on Fluid Mechanics). This description of motion is at the microscopic level. As stated earlier, solving problems of flow in porous medium domains at this level is impractical, as we cannot describe the fluid-solid interface that bounds this flow domain. The construction of a continuum model at the macroscopic level is called for. Accordingly, the averaging procedure is applied to the Navier-Stokes equa- tion (e.g., Bear and Bachmat, 1990), in order to find its macroscopic scale counterpart. Often, the averaged expressions become rather complicated, and various simplifying assumptions have to be introduced in order to eventually reach a relatively simple form of the averaged equation. For example, subject to certain simplifying assumptions, the averaged Navier-Stokes equation is re- duced to Darcy’s law. The averaged equations always contain coefficients that incorporate the information on the detailed configuration of the solid-fluid microscopic boundaries that we have been trying to avoid. For example, in Darcy’s law the permeability is such a coefficient. Although these coefficients have rigorous definitions and structure, their numerical values for particu- lar porous media can be obtained only experimentally. In Subs. 1.3.5, the homogenization technique is utilized to derive some of the above-mentioned results. 54 INTRODUCTION C. Macroscopic heterogeneity We have already introduced the option of smoothing out heterogeneity at the macroscopic level, in order to obtain the description of phenomena at the megascopic level (which may still be heterogeneous). In this case, the av- eraging is performed over a representative macroscopic volume (abbreviated RMV). The characteristic size, ∗, of this volume, is constrained by d∗  ∗  L, (1.3.14) where d∗ is a length characterizing the macroscopic heterogeneity that we wish to smooth out, and L is a length characterizing the porous medium domain. In fact, the features of the REV listed in Subs. 1.3.4B, as well as the constraints imposed on its size, may, at least in principle, be repeated also here, replacing the terms ‘microscopic’ and ‘macroscopic’ by the terms ‘macroscopic’ and ‘megascopic’, respectively. Obviously, the length scale of heterogeneity at the megascopic level will be much larger than that corre- sponding to the macroscopic one. Similar to what happens at the microscopic-to-macroscopic smoothing, here also, the information about the heterogeneity at the macroscopic level appears at the megascopic one in the form of various coefficients that reflect the effect of the actual spatial distribution of the (geometrical) parameters at the macroscopic level on various phenomena of transport. In practice, when we consider field-scale problems, we use the same models as we use for the description of phenomena of transport (e.g., flow of water) at the macroscopic level, but we interpret the coefficients and the variables as average values over an RMV. The only measurable quantities is piezometric head and concentra- tions. Eventually, we use this data with an inverse approach to determine the values of field coefficients and their spatial variability. Consider the case of lenses of one material embedded in a domain com- posed of another material (Fig. 1.3.3); silt lenses in a sandy domain may serve as an example. Let the lenses be of dimensions that are much smaller than the characteristic length of the domain, L, i.e., L1, L2  L, and be randomly distributed in space. Since we do not know the detailed spatial distribution of the lenses, we wish to replace the real heterogeneous domain by a homogeneous one without the lenses. We achieve this goal by defining a representative macroscopic volume (abbreviated as RMV), which is much larger than the scale of heterogeneity at the macroscopic level (i.e., that of variations in permeability, say, spacing between lenses) and averaging over it. By doing so, the permeability of the (macroscopic level) lenticular structure of heterogeneity is replaced by that of an equivalent fictitious homogeneous material having some average, or effective, permeability that is anisotropic at the megascopic scale, with higher permeability in the direction parallel to the lenses. The effective permeability has to be determined by an appropri- ate field experiment. In Sec. 7.5, in connection with contaminant transport, we shall discuss how such a domain can be modeled as a ‘dual continuum’ Continuum Approach 55 L1 L2Clay lens (K2 < K1) K1 L Porous medium domain (enlarged) A representative macroscopic volume (RMV) Figure 1.3.3: An inhomogeneous aquifer treated as an equivalent homoge- neous one using Representative Macroscopic Volume averaging. (or ‘double porosity’) model. If the lenses are non-uniformly distributed, av- eraging will lead to an equivalent heterogeneous domain at the megascopic scale. 1.3.5 Homogenization Although the volume averaging technique discussed above has been widely used (mostly conceptually) for the passage from the microscopic level to the macroscopic one, and from the macroscopic to the megascopic one, let us in- troduce another technique that is generally acknowledged to be more appro- priate for handling multiple-scale heterogeneity. This technique is known as the mathematical theory of homogenization, a term coined by Babus̆ka (1975, 1976-1977). Since the 1970s, this technique has been applied to a range of physical problems that involve composite materials, heterogeneous geological media, and porous media (Bensoussan et al., 1978; Sanchez-Palencia, 1974, 1980; Lions, 1981; Bakhvalov and Panasenko, 1989; Jikov et al., 1994; Mikelic, 2000). Briefly, homogenization is a mathematical technique applied to differen- tial equations that describe physical phenomena associated with a domain exhibiting heterogeneities and/or geometrical features at two scales or more. Figure 1.3.4 shows an illustration of such case. The microscopic (pore) scale, and the macroscopic scale may serve as an example. By homogenization, we obtain a domain which is more homogeneous, at least locally. The co- efficients, which characterize this ‘homogenized’ medium, are referred to as ‘homogenized’, ‘equivalent’, or ‘effective’. In the process of homogenization, 58 INTRODUCTION 0.2 0.4 0.6 0.8 1 x 0.2 0.4 0.6 0.8 1 aΕ a 0.2 0.4 0.6 0.8 1 x 0.2 0.4 0.6 0.8 1 uΕ b Figure 1.3.6: Solution of (1.3.15) with a rapidly fluctuating coefficient: (a) a plot of aε as (1.3.18), with ε = 0.1, in which the dashed line indicates the ‘effective coefficient’ ao in (1.3.39); and (b) the solution (1.3.19) for uε. the asymptotic solution as ε→ 0, u(x) = lim ε→0 uε(x). (1.3.17) To demonstrate this feature, and to validate the result of homogenization, we shall use an example for which we can obtain an exact solution. In the present example, without loss of generality, we shall assume L = 1, and hence  = ε. Let us assume that a is a periodic function, with period ε, aε(x) = 1 1 + 2 sin2 πxε . (1.3.18) For a small ε, say, ε = 0.1, this coefficient is plotted as Fig. 1.3.6a, where we observe ‘rapid’ fluctuations (small period) with a large amplitude. Before applying the homogenization technique, it would be instructional to examine the exact solution of the problem represented by (1.3.15) and (1.3.16). For 1/ε an integer, this solution is uε(x) = ∫ x 0 dx aε(x) ∫ 1 0 dx aε(x) = 4πx− ε sin 2πxε 4π − ε sin 2πε = x− ε 4π sin 2πx ε , (1.3.19) plotted as Fig. 1.3.6b for ε = 0.1. As observed in the figure, and also in (1.3.19), the solution consists of two parts: a slowly varying part (linear), with rapidly fluctuating, small amplitude ‘ripples’ superposed on top. Indeed, the magnitude of the ‘disturbance’ in the solution, caused by the fluctuat- ing coefficient aε, is controlled not by the coefficient’s amplitude, but by its period, which is small. As emphasized earlier, homogenization requires the existence of two scales. At the larger scale, we denote the domain of size L (≡ 1 in the current case) as Continuum Approach 59 Ω, and use the coordinate system 0 ≤ x ≤ 1. At the small scale, characterized by the periodic cells of size  (≡ ε) (see Fig. 1.3.5 for a two-dimensional conceptualization), we denote the repeated domain as Y , and use the scaled coordinate y = x ε , (1.3.20) such that 0 ≤ y ≤ 1 in a Y -cell. With the above definition, we now express the coefficient aε, defined in (1.3.18), as aε(x) = a(y) = 1 1 + 2 sin2 πy . (1.3.21) For uε(x), we can express it as a two-scale function, u(x, y), and expand it into a power series in terms of the small parameter, ε, uε(x) = u(x, y) = u(o)(x, y) + εu(1)(x, y) + ε2u(2)(x, y) + . . . (1.3.22) This is known as the perturbation technique (Nayfeh, 2000). Substituting the above expression into (1.3.15), and applying the chain rule d dx = ∂ ∂x + 1 ε ∂ ∂y , (1.3.23) to the two-scale functions, we can expand and separate the resulting equation into several equations, each corresponding to the same power of ε: O(ε−2) : ∂ ∂y [ a(y) ∂u(o)(x, y) ∂y ] = 0, (1.3.24) O(ε−1) : ∂ ∂x [ a(y) ∂u(o)(x, y) ∂y ] + ∂ ∂y [ a(y) ∂u(o)(x, y) ∂x ] + ∂ ∂y [ a(y) ∂u(1)(x, y) ∂y ] = 0, (1.3.25) O(ε0) : ∂ ∂x [ a(y) ∂u(o)(x, y) ∂x ] + ∂ ∂x [ a(y) ∂u(1)(x, y) ∂y ] + ∂ ∂y [ a(y) ∂u(1)(x, y) ∂x ] + ∂ ∂y [ a(y) ∂u(2)(x, y) ∂y ] = 0, (1.3.26) and higher order equations. The boundary conditions (1.3.16), are assigned to the leading terms, such that the higher order terms take the null boundary conditions: u(o)(0, y) = 0, u(o)(1, y) = 1; u(1)(0, y) = u(1)(1, y) = 0; u(2)(0, y) = u(2)(1, y) = 0; . . . (1.3.27) 60 INTRODUCTION Also, the periodicity of the Y -cells requires that u(1)(x, 0) = u(1)(x, 1); u(2)(x, 0) = u(2)(x, 1); . . . ∂u(1)(x, y) ∂y ∣ ∣ ∣ ∣ y=0 = ∂u(1)(x, y) ∂y ∣ ∣ ∣ ∣ y=1 ; ∂u(2)(x, y) ∂y ∣ ∣ ∣ ∣ y=0 = ∂u(2)(x, y) ∂y ∣ ∣ ∣ ∣ y=1 ; . . . (1.3.28) A quick inspection of the O(ε−2)-equation, (1.3.24), shows that u(o)(x, y) = u(o)(x) (1.3.29) is an admissible solution. Based on a theorem of existence and uniqueness presented in Subs. 4.2.3, as shown in (4.2.67) and the bounding condition below it, (1.3.29) is, in fact, the only admissible (unique) solution of (1.3.24). With the condition expressed by (1.3.29), the O(ε−1)-equation, (1.3.25), can be simplified to ∂a(y) ∂y ∂u(o)(x) ∂x + ∂ ∂y [ a(y) ∂u(1)(x, y) ∂y ] = 0. (1.3.30) To solve this equation, we assume the existence of a Y -periodic function w(y), satisfying ∂a(y) ∂y + ∂ ∂y [ a(y) ∂w(y) ∂y ] = 0. (1.3.31) Comparing (1.3.31) with (1.3.30), and realizing that ∂u(o)/∂x is not a func- tion of y, it can easily be shown that u(1)(x, y) = w(y) ∂u(o)(x) ∂x + f(x), (1.3.32) where f(x) is an arbitrary function of x. Equation (1.3.32) can be differenti- ated with respect to y, yielding ∂u(1)(x, y) ∂y = ∂w(y) ∂y ∂u(o)(x) ∂x . (1.3.33) We can solve (1.3.31) with the periodicity condition w(0) = w(1), to obtain w(y) = ∫ y 0 dy a(y) ∫ 1 0 dy a(y) − y + c, (1.3.34) where c is an arbitrary additive constant. We now turn to the O(ε0)-equation, (1.3.26). Integrating it with respect to y over a Y -cell, i.e., from y = 0 to 1, we observe that due to the periodicity Scope and Organization 63 In Chapter 4, the homogenization technique will be used to derive the (macroscopic) Darcy’s law from the (microscopic) law governing viscous flow in the void space (Subs. 4.2.2), and the equivalent anisotropic hydraulic con- ductivity at the megascopic scale for a layered formation (Subs. 4.2.3). 1.4 Scope and Organization At the beginning of this chapter, we introduced the ‘problems’ that are of interest to groundwater hydrologists and managers: water flow, and contam- inant transport in both the saturated and the unsaturated zones. In both cases, models provide management with information that is essential to de- cision making, e.g., water levels and solute concentrations, to be expected if certain management decisions will be implemented. We have limited the presentation to isothermal conditions, although, under certain circumstances, (man-made, or naturally occurring) temperature variations may have a sig- nificant influence on the flow and solute transport regimes in the subsurface. With the above in mind, Chapter 1 introduces groundwater within the hydrological cycle, and outlines the role of aquifers within the framework of a water resources system. Since the models discussed in this book visualize the porous medium as a continuum, we also explain how the complex solid- fluid(s) domain, called ‘a porous medium’, is transformed into a continuum. Chapter 2 Presents the definition and classification of aquifers and intro- duces the hydraulic approach, based on the assumption of ‘essentially hor- izontal flow’, as an important and practical mode of modeling employed in many cases in practice. Chapter 3 reviews the regional groundwater balance and its components. Both natural and man-introduced components are considered. In each case, both the quantity and the quality of the water are discussed. Chapter 4 is devoted to the basic equation of groundwater motion— Darcy’s law. First, the equation for three-dimensional flow is introduced. Then, the integrated equations for flow in confined and phreatic aquifers are developed and discussed. In doing so, the transmissivity, as an aquifer parameter, and the Dupuit assumption, as a good approximation, are intro- duced in connection with the hydraulic approach to flow in aquifers. The basic motion (or Darcy’s) equation is extended to inhomogeneous fluids and to inhomogeneous and anisotropic aquifers. Chapter 5 leads to complete saturated flow models. It introduces the def- initions of specific storativity and aquifer storativity, and employs these co- efficients to construct mass balance equations for flow in three-dimensional domains. Then, models based on the concept of ‘essentially horizontal flow’ are developed for confined, phreatic, and leaky aquifers by integrating the three-dimensional flow model over the vertical thickness of the aquifer. It is shown how, by doing so, the phreatic and other upper and lower boundary conditions are incorporated as source terms in the equations. Following a 64 INTRODUCTION discussion of boundary and initial conditions, the structure of the complete mathematical statement of any groundwater flow problem is presented. Upon reaching this point, the reader should be able to state, correctly and completely, any problem of groundwater flow in an aquifer, in terms of a partial differential equation and appropriate initial and boundary conditions. Chapter 6 presents flow in the unsaturated zone. This is an important subject as aquifer replenishment takes the form of unsaturated downward flow, and because of the transport of pollutants with this downward flow. The discussion leads to complete well-posed flow models. Chapter 7 presents a comprehensive discussion on the transport and ac- cumulation of chemicals in the subsurface. The main new feature here is the phenomenon of hydrodynamic dispersion. The general mass balance equation for a chemical species, often a contaminant, is developed for both saturated and unsaturated flow. Following a discussion on boundary and initial con- ditions, the complete, mathematical statement of the problem of movement and accumulation of a chemical species is presented. The discussion includes the effects of sources and sinks of chemical species, due to such phenomena as adsorption, dissolution, volatilization, and chemical reactions, under both chemical equilibrium and kinetic conditions. Chapter 8 reviews numerical solution techniques for solving the flow and transport models presented in Chaps. 5, 6 and 7. The objective is to present an introduction to the fundamentals and methodologies of constructing nu- merical flow and transport models and solving them by using computer pro- grams (‘codes’). A brief introduction is presented to some of the more com- monly used codes. Chapter 9 deals with the important problem of seawater intrusion into coastal aquifers, as a consequence of over-exploitation. The discussion leads to complete models that describe this phenomenon. Two approaches are dis- cussed: one that visualizes seawater and freshwater as two immiscible fluids separated by a sharp interface, and another one that regards seawater and freshwater as a single liquid with variable concentrations of dissolved salts and, thus, having a variable density. Some comments are made on the issue of management of a coastal aquifer. Chapter 10 reviews the issues of uncertainty associated with modeling. Sources of uncertainty are examined. The basic definitions of a stochastic process and tools for its analysis, such as sensitivity analysis, kriging, Monte Carlo and perturbation methods, are introduced. Chapter 11 serves as an introduction to the management of groundwa- ter. It discusses the management issues involved, and then introduces tools that can be used by managers for quantitative predictions. Several optimiza- tion methodologies, including the constrained linear programming, the un- constrained gradient search, the genetic algorithm, optimization under un- certainty, and multiobjective optimization, are introduced. Inverse problems, particularly problems of parameter estimation, are also presented. Chapter 2 GROUNDWATER AND AQUIFERS 2.1 Definitions of Aquifers Subsurface water, or groundwater, is a term used to denote all the waters found beneath ground surface. However, groundwater hydrologists, who are primarily (but not exclusively) concerned with the water contained in the zone of saturation (Subs. 1.1.1), often use the term ‘groundwater’ to denote water in only this zone. In this book, we adhere to this definition, using the term subsurface water to denote all the water below ground surface. Practically, all groundwater can be regarded as part of the hydrological cycle (Fig. 1.1.1; or see any textbook on Hydrology). Very small amounts, however, may enter the cycle from other sources (e.g., magmatic water). Aquifer. Todd (1959) traces the term ‘aquifer’ to its Latin origin: aqui comes from aqua, meaning “water”, and fer, from ferre, meaning “to bear”. Thus, an aquifer is a term used for a porous geological formation that (1) contains water at full saturation (i.e., the entire interconnected void space is filled with water), and (2) permits water to move through it under ordinary field conditions. This means that whether a geological formation can be des- ignated as an aquifer, or not, depends on its ability to store and transmit water relative to other formations in the vicinity. Other terms often used are: groundwater reservoir (or basin) and water bearing zone (or formation). Aquitard. This is a semipervious geologic formation that transmits water at a very low rate compared to an aquifer. However, the term should not be used just for any low permeability formation. Instead, the term is restricted to describe a semipervious layer which (1) is thin relative to the thickness of the aquifer underlying or overlying it, (2) has a permeability that is much smaller than that of such an aquifer, and (3) extends over relatively large horizontal areas. In spite of its relatively low permeability, over such large (horizontal) areas, it may permit the passage of large quantities of water between the aquifers that are separated from each other by it. It is often referred to as a semipervious formation or a leaky formation. Aquifuge. This is a rarely used term that describes an impervious for- mation, which neither contains nor transmits water. 65 DOI 10.1007/978-1-4020-6682-5_2, © Springer Science+Business Media B.V. 2010 J. Bear, A.H.-D. Cheng, Modeling Groundwater Flow and Contaminant Transport, Theory and Applications of Transport in Porous Media 23, 68 GROUNDWATER AND AQUIFERS zone, or the root zone, close to ground surface, the intermediate zone, and the capillary zone (or capillary fringe), immediately above the water table. The upper part of the unsaturated zone is called soil water zone. It extends downward through the root zone. Vegetation depends on water in this zone, as plants require both air and water. The moisture distribution in this zone is strongly affected by conditions at ground surface, i.e., seasonal and diurnal fluctuations of precipitation, irrigation, air temperature, and humidity. It is also affected by the presence of a shallow water table. When the water table of the aquifer is deep, it does not affect the moisture distribution in this zone. Water in this zone moves downward during infiltration (e.g., from precipita- tion, flooding of ground surface or irrigation), and upward by evaporation and plant transpiration. Temporarily, during a short period of excessive in- filtration, the soil in this zone may be almost completely saturated. After an extended period of gravity drainage, without additional supply of water at ground surface, the amount of moisture remaining in the soil is called field capacity (see Subs. 6.1.9). Below field capacity, the soil contains water in the form of continuous thin films on the soil particles and menisci between them, held by surface tension. Water in these films is moved by capillary action and is available to plants. Below some moisture content, called the hygroscopic coefficient (= maximum moisture which an initially dry soil will adsorb when brought in contact with an atmosphere of 50% relative humidity at 20◦C), the water in the soil is called hygroscopic water. It also forms very thin films of moisture on the surface of soil particles, but the adhesive forces are very strong, so that this water is unavailable to plants. The intermediate zone extends from the lower edge of the soil water zone to the upper limit of the capillary fringe (see below). Its thickness depends on the depth of the water table below ground surface; it does not exist when the water table is high, in which case the capillary fringe may extend into the soil water zone, or even to ground surface. Temporarily, during replenishment periods, water moves downward through this zone as gravitational water. The capillary fringe (Subs. 6.1.8) extends from the water table up to the limit of the capillary rise of water. Its thickness depends on the soil properties and on the homogeneity of the soil, mainly on the pore size distribution. The capillary rise ranges from practically nothing in coarse material, to as much as 2 m to 3 m and more in fine materials (e.g., clay). Within the capillary fringe there is a gradual decrease in moisture content with height above the water table. Just above the water table, the pores are practically saturated. Moving higher, only the smaller connected pores contain water. Still higher, only the smallest connected pores are filled with water. Hence, the upper limit of the capillary fringe has an irregular shape. For practical purposes, some average smooth surface is taken as the upper limit of the capillary fringe, such that below it the soil is assumed practically saturated (say > 75%). In most regional groundwater studies, the capillary fringe is much thinner than the thickness of the saturated zone and is, therefore, disregarded. The phreatic surface then serves as the upper bound of the saturated zone. Classification of Aquifers 69 Within a homogeneous unsaturated zone, in the absence of infiltration or evaporation, the moisture content generally decreases gradually with height above the phreatic surface. Infiltration will cause saturation to rise temporar- ily as ground surface is approached. In many cases, the spatial variability in soil properties is the dominant factor in determining the distribution of mois- ture content, with regions of fine soils having a higher moisture content, while those of coarse materials having a low moisture content. Numerous complications are introduced into the schematic moisture dis- tribution described here by the large variability in pore sizes, by the presence of layers of different permeability, and by the temporary movement of infil- trating water. All these subjects will be discussed in detail in Chap. 6 Note that although, traditionally, water in the intermediate zone was re- ferred to as ‘vadose water’, in recent years, the term ‘vadose zone’ is often used to denote the ‘unsaturated’ zone, i.e., the zone between ground surface and the underlying phreatic surface. 2.3 Classification of Aquifers The term aquifer was introduced in Sec. 2.1. Let us now introduce the defi- nitions of specific aquifer types. The piezometric head and the piezometric surface will be defined in Subs. 4.1.1. At this point, the former will be defined as the water eleva- tion at a well that has a screened portion at a point within an aquifer. By measuring the piezometric heads at a number of spatially distributed obser- vation wells, tapping the same aquifer, especially with essentially horizontal flow, we obtain a contour map that defines a surface called the piezometric surface. The elevation of this surface at a point in the horizontal plane gives the piezometric head in the aquifer at that point. The phreatic surface is also a piezometric surface (just think of observation wells with screens just below the water table). In what follows, the location of the piezometric surface will be used in the classification of aquifers. A confined aquifer (Fig. 2.3.1a) is one that (1) is bounded from above and from below by impervious formations, and (2) the water pressure in it is such that the level of water in a well that is open in it will be at, or will rise above the upper impervious bounding surface. In other words, the piezometric surface of a confined aquifer is above the latter’s impervious ceiling. An aquifer that is bounded from above by a phreatic surface (as in Fig. 2.3.1b) is called a phreatic aquifer, or an unconfined aquifer. For an aquifer with essentially horizontal flow, the water table is also the piezomet- ric surface of the aquifer. A special case of a phreatic aquifer is the perched aquifer (Fig. 2.3.1c). This is a phreatic aquifer of limited areal extent, formed on a semipervious, or impervious, layer that is present between the persistent water table of a phreatic aquifer and ground surface. A perched aquifer may 70 GROUNDWATER AND AQUIFERS (d) Δ Leakage Leaky confined aquifer Leaky phreatic aquifer Δ (b) (Phreatic surface) Δ Water table Impervious Phreatic aquifer Δ (a) Well Ground surface Screen Confined aquifer Piezometric surface Δ ΔΔ Water tableImpervious orsemipervious lens Recharge Perched aquifer Δ (c) Impervious Impervious Figure 2.3.1: (a) Confined aquifer. (b) Phreatic aquifer. (c) Perched aquifer. (d) Leaky aquifer. exist only for a limited period of time, seasonal or ephemeral, as its water drains into the underlying phreatic aquifer. A leaky phreatic aquifer, shown in Fig. 2.3.1d, is a phreatic aquifer that is bounded from below by a semipervious layer, usually referred to as an aquitard. The latter was defined in Sec. 2.1. A leaky confined aquifer is a confined aquifer, except that one or both confining layers are aquitards, through which leakage may take place. Figure 2.3.1d shows the two kinds of a leaky aquifer. Figure 2.3.2 shows several aquifers and observation wells. Also indicated are the phreatic surface of aquifer A and the piezometric surfaces of aquifers B and C. The upper phreatic aquifer, A, is underlain by two confined aquifers, (B and C). In the recharge area, aquifer B becomes phreatic. Portions of aquifers A, B, and C are leaky, with directions and rates of leakage determined by the elevations of the piezometric surface in each of these aquifers. The boundaries between the various confined and unconfined portions may vary with time, as a result of fluctuations in the piezometric surfaces. Because, in zones with non-horizontal flow, the piezometric head may vary along the vertical, it is important that the screened portion of a piezometer be made relatively short. A well with an elongated screen, will provide a piezometric head that cannot be assigned to a specific point. In fact, when Solid Matrix Properties 73 100 100 90 100 90 80 80 80 70 70 60 60 60 Silt loamLoamSandy loam Clay loam 50 50 40 40 30 30 20 20 10 0 0 30 Sandy clay loam Silty clay Clay Silty clay loam Loamy sandSand Sandy clay 10 20 90 50 40 10 0 Silt 70 C la y (p er ce nt b y w ei gh t) Silt (percent by w eight) Sand (percent by weight) Figure 2.4.3: USDA soil textural triangle. they have been found to correlate well with various important physical soil properties, such as permeability. Various soil classification systems, based solely or partly on soil particle size distribution, have been proposed and used by various organizations. The USDA classification system is based only on the size distribution of particles finer than 2 mm (i.e., sand, silt, and clay). Depending on the fraction of these three size separates, the soil textural class is assigned as shown in Fig. 2.4.3. Another system, commonly used by geotechnical and geological engineers, is the unified soil classification (USC) method originally proposed by A. Casagrande in 1942, and revised in 1952 by the Corps of Engineers and the Bureau of Reclamation (Das, 1983, p. 35). 2.4.2 Porosity and void ratio Porosity, φ, at a point in a porous medium domain, is defined as the volume of void space per unit volume of porous medium at that point φ(x, t) = Uov(x, t) Uo(x) , (2.4.1) where Uo and Uov are the volume of the REV centered at point x, and the void space within the REV, respectively. The porosity depends on the texture and structure of the soil. Soil porosity varies over a wide range of values, from less than 30% to over 90%. Table 2.4.1 gives typical porosity values for a range of natural materials. 74 GROUNDWATER AND AQUIFERS Material Porosity Peat soil 0.6–0.8 Soils 0.5–0.6 Clay 0.45–0.55 Silt 0.4–0.5 Medium to coarse mixed sand 0.35–0.4 Uniform sand 0.3–0.4 Fine to medium mixed sand 0.3–0.35 Gravel 0.3–0.4 Gravel and sand 0.3–0.35 Glacial till 0.1–0.2 Sandstone 0.1–0.2 Shale 0.01–0.1 Limestone 0.01–0.1 Carbonate mud 0.4-0.7 Dolomite 0.001–0.15 Chalk 0.15–0.45 Fractured igneous rock 0.01–0.1 Karst limestone 0.05–0.5 Basalt 0.01-0.25 Table 2.4.1: Typical porosity of natural materials. Sometimes, the void space is made up of two parts: an interconnected portion, through which a fluid can move from any point to any other point within this portion, and a non-interconnected portion. A fluid present in the latter portion of the void space cannot leave it, except by crossing solid phase boundaries. Because we are interested primarily in the transport of mass of fluid phases and chemical species within the void space, unless otherwise specified, we shall use the term porosity to indicate only the interconnected portion of the void space. The terms interconnected porosity and effective porosity are sometimes used for this purpose. One should be particularly careful in using the term ‘effective porosity’, because this term has a number of additional interpretations. For example, in some porous media, the configuration of the (interconnected) void space is such that most of the flow of a fluid phase takes place through only part of the interconnected void space, with only a small fraction of the total flux taking place through the remaining portion of the void space. We often approximate the situation by assuming that in the latter portion of the void space, the fluid is (practically) stationary, or immobile. This may happen, for example, when pores have the shape of a dead-ends, or when very small pores, say between very small grains, are mixed with very large ones. We then assume that the fluid is practically immobile within the dead-end pores and in the very small ones. The term ‘effective porosity’ is used to describe that part of the total void space through which (most of the) flow takes place. The term void ratio, e, defined as Solid Matrix Properties 75 e(x, t) = Uov(x, t) Uos(x, t) = φ 1− φ, (2.4.2) is used mainly in soil mechanics, with Uos denoting the volume of the solid matrix within an REV. The bulk density of the soil, ρb (= mass of the solid per unit volume of porous medium), is defined as ρb = ms Uo = ρs Us Uo = ρs(1− φ). (2.4.3) 2.4.3 Specific surface The specific surface area, Σsv = Ssv/ms, where Ssv is the total surface area of the solid matrix (= solid-void interface) and ms is the mass of the solid matrix, is defined as the surface area of the solid matrix per unit mass of soil. A typical unit is m2/g. It is a very important soil characteristic, especially in connection with surface phenomena such as adsorption and ion-exchange (Sec. 7.3). Fine soils, e.g., clay, are characterized by a huge value of Σsv. To estimate Σsv, consider a soil made up of spherical particles of diameter d. For such spheres, the area per unit mass is given by 6/ρsd, where ρs is the mass density of the solid. For a soil composed of a number of fractions of particle sizes, with mi denoting the mass of solid in the ith fraction, we have Σsv = 6 ρs ∑ (i) mi ms 1 di . (2.4.4) For soil particles in the form of platelets × × b, the specific area is Σsv = 2(+ 2b) ρs  b , (2.4.5) where b is the platelet thickness. For very thin platelets, Σsv ≈ 2/ρsb. For soils made of spherical or cubical particles, with a distribution of sizes, we can make use of the relationship Σsv = 3(Cu + 7) 4ρs d50 , (2.4.6) where d50 is soil particle size of 50% passing (50% of the cumulative distri- bution curve), and Cu = d60/d10 is the soil uniformity coefficient. The above relation assumes that the soil cumulative grain size distribution can be ap- proximated by a straight line on a standard semi-logarithmic plot. For clayey, plate-like particles, the above equation can be modified to Σsv = (2 + β)(Cu + 7) 4ρs d50 , (2.4.7)