A SERIES OF FOUR PAPERS ON WATER FLOW IN DEFORMABLE POROUS MEDIA BY JACOB BEAR, PAUL GOEBEL VISITING PROFESSOR AND M, YAVUZ CORAPCIOGLUj VISITING SCHOLAR DEPARTMENT OF CIVIL ENGINEERING THE UNIVERSITY OF MICHIGAN ANN ARBOR, MI 48109 FEBRUARY, 1981

I I~ ~ ~~~{)n~nd-M

PREFACE This publication contains four papers dealing with the flow of water in deformable porous media. In carrying out the research reported in these papers, we had a practical objective in mind, namely to deal with the land subsidence and the consolidation that takes place in an aquifer as a result of pumping. Examples of significant land subsidence in many places in the world can be cited. As is common in dealing with regional problems of flow in aquifers, our objective was to develop models of consolidation based on the hydrauZic approach, in which the flow is assumed to be everywhere essentially horizontal. In such models, all dependent variables, e.g., piezometric head, horizontal displacement and vertical subsidence, are functions of the horizontal planar coordinates and of time. The justification for employing the (simplified and approximate) hydraulic approach for regional studies stems from the observation that horizontal lengths of interest are much larger than the thickness of the considered aquifer. When this is not so, the problem on hand is a local one and three-dimensional models have to be employed. The procedure for developing the aquifer models is made up of two steps. In the first step, we start with a description of the considered behavior (e.g. motion of mass, transfer of heat, distribution of stresses and strains) within the

continuum of each individual phase (and that includes the solid phase) present in the porous medium system. Then, we average the behavior over a representative elementary volume within the considered porous medium domain. The result is a model made up of a set of equations describing the considered phenomena in the porous medium domain, regarded as a continuum in three-dimensional space. We have demonstrated this step in the case of heat transfer (paper no. 4). In the other cases, this step is omitted and we have started from equations that are already written at the porous medium level. In the second step, the equations obtained in the first one are integrated over the vertical thickness of the considered aquifer, taking into account, in a rigorous manner, the conditions on the top and bottom surfaces bounding the aquifer. In this way, regional models are obtained which have to be solved as boundary value problems in the horizontal flow domain only. The various phenomena (flow of the water, subsidence, energy flux) are usually coupled so that the set of equations has to be solved simultaneously. This usually means that numerical solution techniques have to be employed. As the construction of the mathematical model was the main objective in this series of papers, no attempt has been made to develop these numerical solutions. For some simple cases of a single well, pumping from an unbounded aquifer, an analytical solution of practical interest could be obtained.

The first paper, deals with the problem of saturated and unsaturated filtration through a sample of soil placed in.a centrifuge. The objective of this paper is to demonstrate the treatment of flow in a deformable porous medium using the water balance (or conservation) equation together with the basic (equilibrium) equations of elasticity, as the mathematical model. We have also shown how consolidation affects the hydraulic conductivity. The second paper employs the basic assumption of vertical consolidation only, introduced by C. E. Jacob in the early 40's. A mathematical model in terms of subsidence as the dependent variable is developed for regional subsidence in an aquifer due to pumping. The model is obtained by first developing the equations for a three-dimensional space and then integrating them over the vertical thickness of the aquifer. In the third paper, both vertical and horizontal consolidation is taken into account. The result is a mathematical model in terms of vertical land subsidence and averaged horizontal displacement components, The same approach is extended in the fourth paper to the case of pumping water in a geothermal reservoir (or injecting water into an aquifer for energy storage purposes). In this case, water temperature becomes an additional dependent variable. Because they deal with a sequence of related topics, and at the same time, demonstrate a methodology of treating

a certain class of problems, we have gathered these papers under a single cover, hoping that in this form it will better serve both those interested in the subject of consolidation in aquifers and those interested in the methodology of constructing mathematical models for transport phenomena in deformable aquifers. The research summarized in this publication was carried out in the Department of Civil Engineering of The University Michigan, at Ann Arbor, Michigan, by Professor Jacob Bear, of the Department of Civil Engineering of the Technion, Israel Institute of Technology, Haifa, Israel, who spent his sabbatical as a Paul Goebel Visiting Professor in the Department, and by Dr. Yavuz Corapcioglu of the Middle East Technical University, Ankara, Turkey, who served as a Visiting Scholar. We would like to thank the Department of Civil Engineering, and especially its Chairman, Prof. R. D. Hanson, and Prof. E. B. Wylie, for making it possible for us to carry out this research in the department. To Professor S. J. Wright, to Professor V. L. Streeter for their hospitality and to the most efficient secretary, Suzette Madej who patiently typed and retyped these papers, our thanks and gratitude. M. Yavuz Corapcioglu Jacob Bear Visiting Scholar Paul Goebel Visiting Professor Ann Arbor, MI August, 1980

SECTION I CENTRIFUGAL FILTRATION IN DEFORMABLE POROUS MEDIA -Abs tract 1 -Introduction 2 -Equilibrium Equations for the Solid Matrix 7 -Continuity Equation for the Fluid 11 -Initial and Boundary Conditions 15 -The Complete Mathematical Statement 24 -Solutions and Discussion 25 -Appendix A - Equilibrium Equation for Solid in Multiphase Flow 28 -Appendix B - Effect of Water Compressibility 34 -List of Notations 36 -References 37 -List of Figures 39 SECTION II MATHEMATICAL MODEL FOR REGIONAL LAND SUBSIDENCE DUE TO PUMPING PART 1 - INTEGRATED AQUIFER SUBSIDENCE EQUATIONS BASED ON VERTICAL DISPLACEMENT ONLY -Abs tract 1 -Introduction 3 -Three-dimensional mass Conservation Equations 9 -Conditions on Top and Bottom of Aquifer Boundaries 14 -Integration Along the Thickness of a Confined Aquifer 17 -Integrated Equation for a Layered Confined Aquifer 21 -Integrated Equation for a Leaky Confined Aquifer 22 -Integrated Equation for a Phreatic Aquifer 23 -Subsidence in a Confined Aquifer 24 -Subsidence in a Layered Confined Aquifer 29 -Examples 30 -Summary and Conclusions 33 -References 35 -List of Notations 39 _T.A ci- F 4 o rf i; A'

SECTION III PART 2 - INTEGRATED AQUIFER SUBSIDENCE EQUATIONS FOR VERTICAL AND HORIZONTAL DISPLACEMENTS -Abstract 1 -The Integrated (water) mass Conservation Equation 3 -The Integrated Equilibrium Equations 11 -Integrated Equations for Radial Coordinates 21 -Averaged Equations for a Layered Confined Aquifer 24 -Example: Displacements Due to Pumping from a Single Well 29 -Coupling vs Uncoupling 38 -Summary and Conclusions 39 -References 42 -List of Notations 43 -List of Figures 45 -Comment 50 SECTION IV A MATHEMATICAL MODEL FOR CONSOLIDATION IN A THERMOELASTIC AQUIFER DUE TO HOT WATER INJECTION FOR PUMPING -Abstract 1 -Introduction 3 -Fluid Mass Conservation Equation 8 -Energy Conservation Equation 12 -The Equilibrium Equations 21 -Conditions on the Top and Bottom Surfaces Bounding the Aquifer 24 -The Complete Mathematical Model 37 -Temperature and Pressure Dependent Parameters 41 -Example: A Model for a Single Pumping/Injection Well 42 -Summary and Conclusions 47 -References 52 -List of Notations 55 -List of Figures 58

-1SECTION I CENTRIFUGAL FILTRATION IN DEFORMABLE POROUS MEDIA J. Bear1 and M. Y. CorapciogZu2 ABSTRACT Following a literature survey of studies in which a centrifuge was employed in order to drain a porous sample under enhanced drainage conditions, either for determining a medium's properties, or for a dewatering process, the complete mathematical statement of the centrifugal filtration problem, taking into account porous matrix deformability, is presented. The medium is considered as a combination of an elastically deformable solid material and an incompressible viscous fluid. The flow of air is neglected. The solid matrix forms the skeleton of the porous cake, while the viscous fluid (water) partially fills the pore space. Coupled transient flow (continuity) and quasi-steady deformation (equilibrium) equations are presented. In the mathematical descriptionofthetwo-phase medium considered here, expressions are used for the material variables relevant to the situation. Possible boundary conditions during Professor in the Dept. of Civil Engr., Technion-Israel Institute of Technology, Haifa, Israel. Presently on sabbatical leave at the Dept. of Civil Engr., University Michigan, Ann Arbor, MI 48109 Visiting Scholar, Dept. of Civil Engr., The University of Michigan, Ann Arbor, MI 48109. On leave from Middle East Technical University, Ankara, Turkey.

-lacentrifuging are stated and discussed. The problem is solved for a simple case in order to show flow and deformation of an unsaturated soil column placed in a centrifuge. Some numerical results are presented; they show the changes in pore-water pressure, saturation, permeability, and sample deformation that take place.

-2INTRODUCT ION Centrifugal filtration, through porous materials, both under conditions of saturated as well as unsaturated flow, plays an important role in various diversified fields, such as soil physics, sludge dewatering and sugar and paper pulp drying. To meet these varied needs, a number of studies have been conducted, most of them of an experimental nature. We shall start by reviewing some of the investigations as reported in the literature, first in soil physics and then in connection with the chemical industry. A literature review of the studies associated with soil physics shows that an early use of a centrifuge to determine the specific retention of soils is reported by Briggs and McLane (1907). In certain soil physics and groundwater problems, information is required not only on specific yield, but also on the areal and temporal distribution of water content in the soil following the initiation of drainage. Centrifuging of natural moist soil samples has been proposed as a promising method of estimating moisture distribution. Following the early works of Briggs and others (1907, 1910, 1912), many investigations have been made, and this experimental procedure has been accepted as standard by most investigators (Hassler and Brunner, 1945; Slobod et al., 1951; Miller and Miller, 1955; Marx, 1956). Prill (1961) and Prill and Johnson (1963) made a comparison of the water content of samples after centrifuging

-3and after drainage by gravity for various artificial (glass beads) and natural (Fresno sand) porous materials. They concluded that the values of moisture content obtained by the two methods show a striking similarity, thus indicating the applicability of the centrifuge technique for predicting gravity-drainage relations for sandy materials. In 1963, Johnson et al. reported a detailed study of the centrifugemoisture-equivalent method. In addition to their investigations on various factors affecting the centrifuge moisture content, the data. at different heights within samples of 0.120 and 0.036 mm glass beads and loess (Trenton, Nebr.) were given. For most of the materials tested, the measured centrifuge moisture content decreased with increased distance from the bottom of the sample. The centrifuge moisture equivalent is defined as the moisture content of a soil after it has been saturated with water and then subjected for one hour to a force equal to a thousand times that of gravity (Am. Soc. Testing materials, 1961). The centrifuge moisture equivalent may be adjusted by a correction factor proposed by Piper (1933). The adjusted value, considered to be specific retention, is then subtracted from the porosity to obtain specific yield. In Johnson et al's (1963) experimental procedure, a short column of soil is placed horizontally, with its long axis along a radius, in a centrifuge with a porous plug at the outer end. A similar procedure was previously used by Slobod et al.

-4(1951) in a different study. Corey (1977) reports the use of a centrifuge without the porous plug at the outer boundary and the measurement of saturation using gamma radiation attenuation. In all the above-mentioned studies, only laboratory techniques were developed; no attempt has been made to derive any mathematical expression even empirically. Therefore, the results of experiments cannot be extended to unmeasured physical properties, such as hydraulic conductivity, diffusivity, etc. Alemi et al. (1974) proposed two centrifugal techniques for determining the hydraulic conductivity of cores of natural moist soil. Experimental results are presented for one technique in which the change in weight of one end of the sample, previously centrifuged, is measured with a balance. The mathematical equations describing this redistribution process were developed and fitted to the data to ascertain the soil water diffusivity. The value of the hydraulic conductivity was obtained indirectly. The second technique for which a theory is presented,but no experimental values are given, depends upon the measurement of the outflow from the sample when the speed of centrifugation is suddenly increased. Their theoretical study contains certain assumptions and boundary conditions which are questionable. For example, pore pressure head distribution at the end of centrifugation was assumed to be the same as in the uniform rotation of a fluid without considering the porous soil structure. Water content equation during redistribution was solved by making use of this boundary condition.

-5Important applications of centrifugal filtration take place in the chemical industry (e.g., Purchas, 1971). Studies in chemical engineering have aimed mainly at the estimation of the rate of filtration. An analysis of the centrifugal filtration given by Grace (1953) is representative of several studies in this category. In Grace's study and in other ones (Storrow, 1957; Bingeman and Coates, 1960), experimental data on centrifugal filtration are correlated by a simplified flow rate equation. Grace (1953) concludes that, a general expression giving the compressive stress in a porous medium as a function of radial distance is needed for a rigorous analysis. Storrow (1957) derived his similar formula for a rigid saturated porous medium. Later, Bingeman and Coates (1960) reached the same expression by using a much simpler approach. In these studies, the filtration rate equations were obtained from simple balance equations rather than by coupling flow and deformation equations. In fact, no account was taken of the medium's deformability. Even the basic principle of mass conservation for the fluid was not utilized. Hence the expressions obtained were oversimplified. Since the material is deformed also by the forces associated with the flow, the flow and deformation processes are coupled; any analysis of such a problem requires the simultaneous solution of the flow and the deformation equations. Joo and

-6Lederer's (1974) formulation appears to be the only one considering the fundamental equation of flow in a compressible porous medium. However, a flow term due to the centrifugal force was not included in the formulation, nor do they consider the balance of stresses acting on the medium during rotation, coupled with the flow equation. An important part of centrifugal filtration is also the final dewatering or drying phase, following the cessation of liquid feed. Of particular interest is the level of undrainable (residual) moisture which can be achieved, and how it is affected by various factors such as the angular velocity and the thickness of the porous layer. Batel (1961) presented some studies on this subject, with particular reference to the removal of water from coal. From the review presented above, one may conclude (a) that most of the work has been experimental, (b) that porous medium deformability has not been properly included in the analysis and (c) that no serious attempt has been made to couple correctly the flow and the deformation of the porous matrix in the analysis. The main objective of this work, which is primarily a theoretical one, is to attempt to present a complete mathematical statement of the centrifugal filtration problem, taking into account porous matrix deformability.

-7Three cases will be studied: (A) Unsaturated flow in radial coordinates, (B) Saturated flow in radial coordinates, and (C) Unsaturated, one-dimensional (along a radius) flow. For the last case, a numerical solution is also presented. Its purpose is to obtain some numerical results that show an example of what happens along a column placed in a centrifuge. EQUILIBRIUM EQUATIONS FOR THE SOLID MATRIX Figure 1 shows a cylindrical sample of soil (or, in general, a porous medium) placed in a centrifuge. Its width is R2 - R1, where R2 and R1 are the outer and inner radii;respectively. The height of the sample is H. The outer boundary, at R2, is a perforated basket, or a porous plate, such that water can flow through it with practically no resistance. At a later stage we shall also consider a boundary condition with a specified resistivity of this porous plate. Initially, the soil sample may be either fully or partly saturated with water. In the unsaturated case, the initial moisture distribution is usually uniform (neglecting the effect of gravity, as H is relatively small). When initially saturated, we may continue to maintain saturated flow conditions by a constant feed of water, maintaining a layer of liquid (inner radius R0, outer radius R1) as shown in Figure 1. The centrifuge rotates at a constant angular velocity w around a vertical axis. The driving force, which varies with the radius, due to rotation, produces filtration through the

soil sample. The compressive stress acting on the solid matrix (effective stress, integranular stress; Terzaghi; 1923) at any point in the sample is affected both by changes in the pore-water pressure (as in the general case of flow through a deformable porous medium) and, in the particular case studied here, also by the centrifugal force acting on both the solid and the water. In Appendix A, following Bear and Pinder (1978), the total stress a at a (macroscopic) point in the porous medium domain is shown to be the sum of the intergranular (or effective) stress a' and the pressure, p, in the water that in saturated flow, fills the entire pore space (porosity, n) and only part of it (at saturation Sw) under unsaturated conditions a = a' - S p [1] In [1], I is the unit tensor, and each stress component is per unit area of soil sample. A change in effective stress produces a deformation of the solid matrix, resulting in a change in n. In unsaturated flow, water is released from storage in any unit volume of porous medium as a result of three processes that occur simultaneously: the decompression of water, the reduction of porosity due to an increase in effective stress and a reduction in saturation. Often the water is assumed to be incompressible. In saturated flow, water is usually taken as compressible, and the saturation remains constant (S = 1). Because we assume no shear w stress in the water, the total shear stress is carried by the solid matrix only.

-9 - The total stress acting at a point in a porous medium domain satisfies the conditions of equilibrium. The latter is the (macroscopic) momentum balance. When the inertial terms are assumed to be negligible, in view of [A-4] and [1], the equilibrium conditions may be written in the form V-(' - Sw P I) + F = 0 [2] where F is the sum of all body forces. For the special problem of the centrifuge, we rewrite [2], using cylindrical coordinates, in the form (Biot, 1941; Verruijt, 1969, except for the introduction of Swp instead of p) - (C S p) 1 aC, I a C' - a rr w - r rz rr F =,:+ - + + F 0 [3] 3r r Me 3z r r,a%' 1 a(a - SwP) az CCr + r + + 2 0 [4] 3r r 3e az r' +1 a'e C(a' - S wP) a' zr +_ z zz w zr__z + +z - r= [5] 9r r 38 3z - where we have neglected the body force due to gravity with respect to the centrifugal force due to the rotation. The latter is 1,000-3,000 times larger than the former. Hence, no body forces exist here in the 0 and z directions. Due to the conditions of axial symmetry that we have in the rotating centrifuge (a' = a' = 0), [3] through [5] reduce to the single equation a~' - a' - ~03 (S,p) rr rr (Sw) r + r FF = [6] arr r r r In writing [6] we have also assumed that stresses and deformations are independent not only of 0 but also of z, which leads to o' = 0. zr

-10The body force F per unit volume of porous medium, resultr ing from rotation, can be expressed as Fr = [nSwPf + (1 - n)ps]w2r [7] (compare with [A-9] where only gravity is considered). In [7] Pf and p denote the densities of the water and of the solid, respectively; n denotes the porosity of the porous medium, and w is the angular velocity. Under the above stated conditions, the strain components are related only to the displacement Ur in the r direction DU U E:r; r [8] rr Dr ee r This follows from the general relationship ~ = V*u. We now assume that the solid matrix is elastic, homogeneous and isotropic, and use Hooke's law in order to relate the components of the effective stress, to the components of strain o' = 211' + X' E [9] rr rr aoe = 2P' Ce + X' [10] a' =X' E [11] zz where X' and i' are the (macroscopic) Lame constants of the porous medium as a whole, i.e., not merely of the solid phase (see any text on the theory of elasticity). With these stress-strain relations, the volume strain C is given by DU UU rr1 ~ (r ir r rr ~ r +r r 9r r 12]

-11By substituting the strain relations [8] into [9] through [11], the following alternative form of Hooke's law is obtained U ar + r r C (X' + 2t') r + r [13] rr Dr r U a U'd = (X' + 2p ) r+ [14] DU U ar +, r + r [15] zz Or r Inserting this form of Hooke's law into the equilibrium equation [6], gives a2U 1 aU U a(s p) (2k' +' 2+ r - -r 2] - [nSw pf + (1 - n)p ]c2r = O:r 2 r Dr r 2 Dr w f s [16] or in terms of volume strain E as 1 D(SWP) w2r[nS p + (1 - n)Ps] 3r (2p' + X') r (2' + X') [17] Equation [17] contains four variables: ~, p, n, and S. We need additional equations to obtain solutions for these variables. CONTINUITY EQUATION FOR THE FLUID Although in unsaturated flow we have the simultaneous flow of both air and water, we shall assume here that the air is stationary and at a constant pressure equal to the atmospheric pressure (= 0); only the water flows. This constraint may easily be relaxed if so desired. We shall also assume that under the pressure conditions in the centrifuge, the water is practically incompressible, pf = constant,

and the solid is incompressible, i.e., pS = constant. However, the solid matrix as a whole is deformable, with n = n(r,t). One can easily introduce water compressibility (Appendix B). Under these assumptions, the continuity equation for the water is given by a(nS ) V-q + = 0 [18] where q is the specific discharge of water with respect to the fixed coordinates. It can be determined from Darcy's law, which expresses the relation between fluid pressure and specific flux of the fluid with respect to the solid. In the case of a deformable porous medium, the solid grains also move at a velocity V. It is the specific discharge of the water relative to the solid, qrek, qre = q - (nSw)V [19] that is expressed by Darcy's law. Inserting this relative flux expression in [18], we obtain d n dS V-qre + S + n + nS VsV = [20] -re wdt dt w -s d s 0 where dt = t + V.V denotes the total derivative with respect to moving solid. The fourth term on the LHS can be obtained from the equation of mass conservation of solid, assuming the density of the solid to be constant V-[(1 - n)V] + (1 - n) = 0 [21] ~S Dt

-13Equation [21] can be written in the form d n VVs - 1n dt [22] Therefore, [20] becomes S d n dS Vqre+ 1 - n dt dt[23] Another expression for V.V is V-V = [24] us 5 t where the right hand side expresses the rate of change of the volume strain E. By combining [22] and [24], and, as an approximation, replacing the total derivative of Sw and by a partial ones, i.e., aS w assuming V.VS << and V ~ - [23] reauces to -Vs w at -S to DS Veq + S n 0 [25] -re9 w t ht One should remember that the capillary pressure Pc = Pair - Pwater Since we assume Pair = 0, we have Pc Pwater E - p. The capillary pressure itself is a function of the saturation, Pc = Pc(Sw)' The graphic expression of this relationship is the retention curve, which varies from one soil to the next. This relationship exhibits hysteresis, a fact which for the sake of simplicity, we neglect here, assuming only drainage everywhere. The relative specific discharge is expressed by Darcy's law. Three driving forces (per unit volume) are acting on the water in the centrifuge: a pressure gradient w2r 2 - Vp, a centrifugal force - Vpf 2 and a gravity - pfgVz, which we neglect; lr is a unit vector in the direction of the radius. Hence k k qres = oreZ (Vp + Vpf 26]

-14where k is the medium's permeability at complete saturation and 0 kre kr e(S ) is the relative permeability (O < k < 1), I is rek re w re' the dynamic viscosity of water. We note that in a deforming porous medium ko varies with the variation in porosity. By inserting [26] in [25], we obtain ko krek 02r2 DS V [- k ret V(p +Pf 2 + S - n( ) D [27] The terms DSw/Dpc can be obtained from the retention curve of the soil. We have here three dependent variables: p, Sw and ~, and two variable parameters: n and k re. In order to solve them we have equations [17] and [27], and a relationship between p and Sw in the form of a retention curve for the considered soil. In a deforming porous medium, we also need, (a) a relationship n = n(p) (actually of effective stress which, in turn, is related to pressure and saturation, (b) a relationship kreR = k (Sw), and re9i w (c) a relationship ko = k (n). 0 0 The system of equations developed above describes the centrifugal filtration through a deformable unsaturated porous medium in any coordinate system. Filtration in a saturated porous medium is a special case of this process. When we replace Sw = 1 in [17] and [27], these equations reduce to 1 a P- w2r[npf + (1 - n) ps] +aE = 1 a p~~~ ~ [28 r (21' + X') r (2 + ) (2H' + X') k V(p + r = 29] which describe the filtration process in various industrial operations under saturated flow conditions. The dependent variables in

-15these equations are E and p. In a deformable porous medium n = n(p), ko = k (n). In laboratory experiments to determine the specific retention (e.g., Johnson et al., 1963) and hydraulic conductivity (e.g., Alemi et al., 1976) of soil samples, a short column of soil is placed horizontally along a radius in the centrifuge, and rotated (see Figure 2). In this type experiments, the flow and stresses are essentially one-dimensional, and [17] and [27] should be rewritten in one-dimension only, rather than in radial coordinates. For simplicity, we shall denote the single dimension here by r. Also it is more convenient to express ~ in terms of the displacement in one-dimension. We obtain kk 2U S D _ okrekp 2 rP DP re p + P wr)] + S - n( = 0[30] Dr' Dr f w Dr Dt 2Ur M1 a(SwP 2r[pfnSw + P (1 - n)] r _ 1 -w + — r[fnS - n)- [31] ar2 (2pu' + X' ) Dr (2p' + A' INITIAL AND BOUNDARY CONDITIONS Initial Conditions Initially, a porous medium sample may either be saturated or it may have a certain water content. In the case of saturation, the initial condition Sw(r,O) = 1 and p(r,O) = O [32] Unsaturated samples may have initially a uniform water content or a certain water content profile. Then Sw(r,O) = C1 and p(r,O) = C2 [33]

-16or Sw(r,O) = fl(r) and p(r,O) = f2(r) [34] where C1 and C2 are known conrtants and f and f2 are known functions. Also, we have to specify initial distributions of porosity, permeability and other matrix parameters, including the displacement Ur, since it is also a variable n = n or n(r,O) = f3(r) and U (r,O) = 0 [35] where n is the porosity of an initially homogeneous sample and f3 (r) is a known function. Boundary Conditions for Pressure Boundary conditions are defined by the particular physical conditions imposed during an experiment. In general, at any boundary, we have [q - eu]..VF = 0 [36] 1,o where [A]i.o = Ai - Alo denotes a jump from the inner side i, to the outer side, o, of the boundary; q is the specific discharge with respect to a fixed coordinate system, 0 = nSw is the water content, u is the velocity of the boundary, and F = 0 is the equation of the boundary. In the one-dimensional case considered here, for the inner boundary we have F1 = r -Rl(t) = 0. For the outer boundary F2 = r - R2 = 0. For saturated flow 8 = n. A similar equation should be written also with respect to the solids (for a constant p5)

-17[(1 - n)(V - u)] 0VF = 0 [37] If we assume that the boundary is a material surface with respect to the solid, i.e., (l1 )(V - u)I.VF = 0; (1- n)(V - u) F = 0 [38] s i s ~ we can combine [36] and [38] to yield an expression in terms of the flux of the water relative to the solid, qreq = q - 9V' which is the flux expressed by Darcy's law [q ~rei *VF = 0 [39] [qreQ]i,o [ When the outer side is completely impervious, i.e. 8outside= and q|outside = 0, the condition at an impervious boundary (which is also a material surface with respect to the solids) becomes qr VF = 0 [40] When an outer side contains no porous material, (Vsloutside 0) but an accretion at the rate qig = N-VF is supplied there to the boundary, we have: (q re -N) VF = t [41 where d- = - u -VF = 0. Here VF = + lr. dt 9t

-18The following are several examples of possible boundary conditions. A) At the outer end of a sample, r = R2, we may have an impervious boundary, i.e., no flux normal to the boundary, i.e., qrek = 0 at r = R2 or in terms of pressure qret - k -r PfWr) r=R =2 3Pr = - pfL2R2 [42] 3r r R2 f 2 which is a second type (Neumann) boundary condition. B) At the outer boundary we have a thin porous plate which offers some resistance to the flow of water through it. This may happen in some industrial processes and laboratory techniques where the outer boundary of the sample is a filter cloth or a porous membrane. At such a boundary we require that the outflow at the boundary r = R2 be equal to the (assumed saturated) flow through the membrane. Hence (P + c02r 2 3(P 2r 2 k L + 2r f'P 2 f) r=R2 P 2 Pf) r=R2+d' k ap + W2rpf ) k' (r f r r=R2 11 d' [43] where k' and d' are permeability and thickness of the membrane. This leads to the condition k r 2 k'2 _p l 2R2pf) = R2d2P 2- d' Wcp [44] 1 r=zf,'/k' r=R2 with d' << R2, k' << k, and d'/k' = membrane impedance,we obtain

-19Up r 2= _ R2 r kd' I/k' =- PfR2 r=R 2 which is a third type (Cauchy) boundary condition. C) Since in the absence of any semi-previous membrane, the water at the external side of the outer boundary is already at atmospheric pressure (p= 0), we have the condition PIr = R 0 [46] 2 which is a first kind (Dirichlet) boundary condition. Unlike conditions (a) and (b), this last condition is valid only for saturated flow. For unsaturated flow, we do not know the pressure just inside the boundary. We shall therefore employ the balance condition described below. D) At the outer boundary, if water can freely leave the sample as a result of centrifuging, then the flux qrek is equal to the total change in moisture content along the sample. Hence = d 2 rer =_ R 2 j J(r,t)dr; 6 = nSW [47 =2 R1 (t) or with Darcy's law (Pp 2 R2 = R ~ -R1 3-r rr = R2 f 2) (R2 1 3t [48] where 8(t) denotes the instantaneous average water content along the column

-20e R - RTP J(r,t)dr [49] 2 1 J i(P2 - I 1 and k - k(I ). r=R2 2 In developing [48], we have made use of Leibnitz rule of differentiating with variable boundaries. The term R 1/3t is, in fact, related to the displacement at the inner face of the sample. For a small finite period At, we obtain _ r (.r2 = R 2pR)At (R2 -R1) r = R r r = R (~R1/at)At = Ur=R [50] Equation [50] expresses a balance during a small period of time At. In numerical solutions, because of the various approximations involved, we should check that this equation is always satisfied. E) At the inner face, r = R1, we usually have a flux condition. As stated in the introduction an important application of centrifugal filtration is the drying and dewatering of porous materials. In this case there is no layer of liquid over the porous medium. This is also the situation after a liquid feed at r = R1 has ceased in the filtration through a saturated material. At the inner end, r = R1, energy is supplied to the sample, producing evaporation from the sample. When the rate of energy supply is known, the rate of evaporation, qQ, is also known. The boundary condition is then, from [41], we have qe = a-rp + 22pfR1); k = k( r=R ) [51] sa~~~~~~~~~rR

-20 a Equation [51] is valid provided the soil, under the condition of e(r) and 3p/3r that develop can actually supply qi to the boundary. Otherwise q. is an upper limit and the actual flux will be less than q; part of the energy will be used up to heat the soil surface or produce evaporation at the microscopic air-water interfaces and vapor flow (but this is beyond the scope of this paper).

-21Often we approximate [51], by assuming q = 0. Then k r I = R + ) P = o; k = k(l ) [52] sR Dra r = r = R1 or r r=R + Pf 2R1 = 0 [53] which is a second kind boundary condition. We recall that in the above equations we may replace ap/ar by (Op/3e) (Oe/3r), where ap/3e is a known function from the retention curve. Another possibility is that a certain flux qi is supplied at r = R1, but at a rate which is insufficient to produce ponding on the sample's surface, then qi (a r = R1+ WPfRl) k = k(e R 54 1' r=R11fr is the boundary condition at the inner surface. One should note that no information on p or Sw at r = R1 is available to enable the specification of a first kind boundary condition. When supply qi exceeds infiltration capacity, we have ponding. F) At the inner face, a constant (or time dependent) supply of water may maintain a layer of water adjacent to the porous medium. This is similar to the layer of water shown in Figure 1 for the radial case. In many industrial filtration problems, liquid is continuously added from the center of the centrifuge, maintaining a layer of liquid with constant (or time dependent) thickness over the porous medium. Then,assuming~ that the liquid layer like the soil also rotates about the vertical axis at the

angular velocity a, in a forced-vortex motion, we obtain P _ pf 2 (R - R2) >O [55] where R is the distance to the surface of the liquid layer from the center of centrifuge. This is a first kind (Dirichlet) boundary condition. At such a boundary, we also have Sw = 1. One should note that all second kind boundary conditions may also be expressed in terms of Sw since a ap r where DP is obtained from the retention curve, assuming independence w of variations in porosity. Boundary Conditions for the Displacement Ur G) At the outer boundary r = R2, there is no displacement due to the rigid basket, or fixed holder. Hence the boundary condition is Ur = R = 0 [56] 2 H) At the inner face, we have a moving boundary R1 = Rl(t) due to the deformation, or compaction, of the soil sample. Without a liquid layer on this boundary, we assume that the effective stress vanishes, i.e., a' = 0. Hence from [13], which for the rr one-dimensional case reduces to aU a' = t(2' +' ) [57] rrD we obtain the condition a ar = [5s]

-23This condition is valid for both saturated flow without a water layer and for unsaturated flow, if we neglect the effect of the (negative) pressure in the water on the effective stress. However, if in the unsaturated case we do take the pressure in the water into account, the total stress on the surface r = R1 is indeed zero, but because of the pressure in the water o' = Swp (Note rr w that this only is true if a = 0). This means that the boundary condition is Swp ur SWP [59] r (2i' + X') However, one should think of this detail as belonging to the microscopic level of consideration. At the macroscopic level, we have to assume for the exposed solid surface in the r direction that o' = 0. For the general radial case, the condition is rr derived from [13] in which we insert a' = 0. rr I) In the case where a water layer exerts an external stress (compression) on the water-solid surface, we obtain the effective stress (at the boundary S = 1 andp is expressed by [55]) 2 Pf rr 2 (R2 - R2) [60] rr 2 1 o Using [1] and [55], we obtain for the present case a' = 0 [61] rr The condition for U is therefore r r 0 [62]

-24This completes the discussion on the various possible boundary conditions. In most cases, the boundary conditions remain valid also for a soil sample between two concentric cylindrical surfaces. In other cases, the required modifications are obvious. THE COMPLETE MATHEMATICAL STATEMENT A solution is sought of the system of equations describing the flow through a deformable unsaturated porous medium as it occurs in certain laboratory experiments (e.g., Johnson et al., 1963). Since an analytical solution is not possible. A numerical solution is employed. The equations to be solved are [30] and [31], subject to the boundary conditions [50], [53] for p and [56] and [58] for Ur. In addition, we need information on the relations kO = kO(n), n = n(o'), and Pc = Pc(Sw). The medium's permeability, ko, at complete saturation, and the relative permeability, kret will be expressed in the example studied here by -7 n k = 1.06 x 10 [63] 0 (1 - n)2 S -S 3 k ( wo [64] WO Equation [63] is known in the literature as Kozeny-Carman equation; SWO is the irreducible saturation. The quotient (aSw/p9c) will be approximated from the following approximation of the retention curve

-25S = - 0.062 p + 1.00 [65] W C The relation between porosity and pressure will be given by n = n - (1 -n )(- po) (S + p dSw/dp)/(2' + ) [66] where no and po denote initial porosity and pressure, respectively. Initial conditions for p and Ur are given by [33] and the second equation of [35]. Equation [66] is obtained by making use of compressibility definition and [A-16]. As well known, the compressibility is defined as - 1 1 d(l - n) (2p' + X') 1 - no da' SOLUTIONS AND DISCUSSION Because the objective of the numerical solution presented below is just to get some feeling for the drainage process that takes place in a centrifuged soil sample, a very simple numerical scheme was employed. No attempt has been made to make use of a more sophisticated scheme and/or a more refined grid. Some of the equations are non-linear. In what follows, they have been linearized, to facilitate a simple solution. Equations [30] and [31] are replaced by the following algebraic finite difference equations, using an implicit scheme and linearization 1 (kkre)i+l (kokr)i t Pi+l Pi 2 t+ or i ] + P W ] -Ax Ax f i - [(kokre t Pi+l 2p + Pi.1 t+l t [ rez- [' ] + [(S ) t (Ax) 2 U - U. t+l U. - U t i+l i i+l i 1 t+l t [()St t Pi -p - [(n) i]t [(a) i At = 67]

-26Ui -2U. + U t+l t pi+l Pi t+l (2p' + C' ) I I * i i- =[( Sw)i] [ ax ] (Ax) 2 A + 2x[(nS ) p + (1 - n.i) ] [68] 1 w i f 1 s For the numerical calculations, the following hypotetical values were used: po = 5 gr/cm, n= 0.37, S = 0.69, X = 10 1/sec, R1 10.5 cm, R 15.5 cm, p = 2.5 grm/cm, and 1/(2p' +' ) = 4 x 10-6 cm2/gr. Spatial mesh size Ax and time step At were taken an 1 cm and 0.25 seconds, respectively. These values are realiable for soil samples. The results of the computer runs for solving [63] through [68], are presented in graphical form in Fig. 3. The changes of porosity are not shown in the figures since it was found that these changes are negligibly small during centrifuging (under the conditions considered here). For example, the porosity changes were less than 0.01% after one hour of centrifuging. As seen in Figure 3, in the course of centrifuging, water gradually drains from the soil sample. The reductions in p, Sw and k are larger during the early stages. Note that we assumed that w reaches its constant value X = 10 sec at no time. -1 In reality, a certain time may elapse until we get to X = 10 sec At large times, the decrease is smaller and asymptotic values are reached. As is stated earlier, eventually, the porosity changes are negligibly small, so that the changes in medium's

-27permeability caused by changes in relative permeability k -only. reZ About 85% of the displacement takes place at the very beginning of centrifuging, this is due to the immediate application of a large centrifugal force. The latter is independent of time, but varies with radius. This is shown in a dotted curve in Fig. 3. The rest of the displacement is due to the decrease in pore-water pressure, which in turn is balanced by an increase in radial intergranular (effective) stress transmitted through the grain-to-grain contacts of the porous medium. For the particular values of X' and I' selected here, the compressibility is very small during centrifuging. Nevertheless, it may be quite large for soft materials such as sludge and paper pulp in dewatering and drying processes. It may be concluded that the theory presented in this paper satisfactorily and correctly simulates centrifugal filtration in deformable porous media. The change of pore-water pressure, displacement, permeability, degree of saturation, and porosity can be predicted by this model. Numerical solutions were given to show what happens along an unsaturated soil column placed in a centrifuge. A more sophisticated, detailed numerical technique would improve the accuracy of the results. This will be the subject matter of a subsequent paper.

-28APPENDIX A EQUILIBRIUM EQUATION FOR THE SOLID IN MULTIPHASE FLOW Bear and Pinder (1978), in dealing with porous medium deformation in multiphase flow, developed the macroscopic equilibrium equations by volume averaging the microscopic equilibrium equation for each of three phases (solid, oil and water) present in a porous medium domain, and summing the resulting equations. They employ the following two kinds of averages for any property G(x') at a microscopic point x' < G > (x,t) = U G(x',t) dU(x') [A-l] o (U) called phase average, and G(x, t) G(x',t) dU(x'), < G > = e G [A-2] Oaa called intrinsic phase average, where UO is a representative elementary volume (REV) around a (macroscopic) point x, and U oca is the volume of a considered a phase inside Uo; e = U a/Uo is the volumetric fraction of the a phase. The selection of the kind of average (= macroscopic) value to be used in each case depends on the measurable quantity in that particular case. They use two averaging rules (Gray and O'Neill, 1976) (a) < VG > = V< G > + ( — GdA; o (A),e VG=v V {+- v GdA [A-3] rv ~ ~~~ v GdAr

-29<G ><G> o (b) <sa-t at U- (A dA; Uo (A) _at at Uo J Gu) v dA [A-4] where Aa is the total surface separating the a phase from all other phases in UO, u is the velocity of that surface, and v is the outward unit vector on that surface. By applying the volume averaging procedure to the (microscopic) momentum balance equation written separately for each of the three phases (a = s,o,w) present, and adding the three resulting (macroscopic equations), they obtain: <p f > + <p f > + < p f > + V- < T > V- < p I - V' < p I > <Polo ww sas -s on > + ~[T] -v dA + 1T v [] ]v- dA = o0 o A ) o,w -o U W o 0 )WA ow ws so [A-5] where f is the body force on the a phase, T is the stress in the a phase, pa is pressure (positive for compression) in a fluid a phase, Aa is the surface of contact between phases a and 8 in UO, [T ] T - T = jump in T across A and I is the unit second rank tensor. In writing [A-5], they have already assumed that (i) inertial effects are negligible and can be neglected, (ii) all interfaces separating phases are material surfaces. Assumina also that (i)UO I [T] V dA = 0, - [T] v dA = 0, OW WS Aw

-30(ii) in the absence of shear stress in the fluids, ~1I r vdA [P] Iv dA= I (p P I P vdA P. P...dA P= (W) - ( = [A-6] ow where Pc = p0 - PW is the capillary pressure and f'(0 ) is a function of the volumetric fraction of e; f'(8), dims.L1 is a vector that gives the total Ao area per unit volume of porous medium. This introduces the effect of capillary forces into [A-5]. (iii) Po, Pf' Ps are constants; f4= g With these assumptions [A-5] reduces to (0 p + 8 p + 8 p )g + V-(<T > - e80 I- WP I) + p( w)f'(ew) =0 o o w w s o ~- s woo c w.w In saturated flow, this reduces to [A-7] {npw + (1 - n)ps}g + V-(< % > - np wI) = O [A-8] For unsaturated (i.e., air-water) flow, where e8Pa<<~wwp and assuming Pa = 0, we obtain {nS W + (l-n)p }g + V(<T > - nSwPwI) + P (0 )f'(08) = 0 [A-9] were 8 E nS; S = degree of saturation. In this paper.we neglect the effect of capillary pressure. then F + V (< T > - nSwPwI) = 0 [A-10] where F represents the total body force on the porous medium and

-31(< Ts > - nSwPwI) represents the total stress. In the text < T > a o' is the effective stress and p p is the pressure in -S w the water. In spite of the above analysis, it seems that in granular materials, where the area of contact between grains in a very small traction of any cross-section, the common practice is soil mechanics is define the integranular stress, a' for saturated flow by a' = a p I [A-11] where Pw is positive for compression, rather than by a - nPwI as indicated by the above analysis. The common practice in soil mechanics (verified by numerous experimental and field evidence) can thus be interpreted as assuming that solid matrix deformation is caused only by the stress in the solid matrix < T > minus the isotropic effect of the water pressure surrounding (practically completely) each grain. It is assumed that this pressure is equal to the pressure in the water at that point. Hence; the effective stress is defined by' = < s > + (1 - n) pwI [A-12] -S wThis can be obtained from [A-10] by writing a = < T >- np I = (1 - n) + (1-n)p I - (l-n)pI - np I -n)(T - WI) - Ia' I [A-13 W n)( - p,I) I = a' - pI [a-13]

-32as in [A-ll]. The effective stress causing matrix deformation is thus defined as (1 - n) (T - PwI). Thus, in considering saturated flow, we shall adopt here the definition of a' as defined by [A-13]. In unsaturated flow, we may follow similar considerations. We shall assume that the average equivalent pressure acting on the entire surface of the grains is given by Swp I, where it assumed that ratio of solid-water surface area to total solid surface area is equal to S. Hence: the effective stress is w < T > - SwPwI. Neglecting air pressure, it can also be derived from a= < T > - np I = I T -e p I + (1 - n)S p I - (1 - n)S pI ~~ w s - s s w w w w ww (1 - n) (s - SwPwI) - pWI [9 + (1 - n)Sw = a' S I [A-14] where e = nS, = (1 - n), = n(l - Sw): e + e + e = 1 w w s air w' w s a Thus [A-14] defines the effective stress a' and relates it to the water pressure and to total stress in the case of unsaturated flow.Eq. [A-14] will be used in the present paper. Finally, we may add that some soil mechanics investigators (e.g., Bishop, 1960) propose a=, - XP I [A-15] s -w where X = xi(S ) and depends also on the percentage of clay in w the soil. Although this relationship is empirically found to

-33be non-linear, it seems that for sandy and loamy soils, X = Sw is a good approximation. For clayey soils, more accurate expressions have to be introduced for X(Sw). Lee (1968) reviews this subject in detail. The analysis presented in this paper for X = Sw is just an example and can be extended to other cases of X (Sw) recalling: — that Sw = Sw(pw). If the total stress is assumed to be constant, then [A-15] yields day dS da'- S + [A-16] dp wis determined from the retention curvedp The quotient (dSw/dp) is determined from the retention curve.

-34APPENDIX B EFFECT OF WATER COMPRESSIBILITY The basic mass conservation equation, for both saturated (S = 1) and unsaturated flow is w Vpq + 3(npSw)/Dt = 0 [B-1] where q = nS V = nS V and q q + nS V. Another form of ~ w w w r W S [B-l] is with V-V defined by [22] and d ( )/dt = a( )/at + V V( ), [B-1] can be rewritten as nS d p nd S S d n Veq + s w + w s [B-2 -r dt d 1 - n dt2] We now define water compressibility B' = (1/Pw)dwP w/dp and moisture capacity for a deformable porous medium C'(Pw) = dsSw/dPw, and obtain dw dp S dn V-qr + nSw Bdt + nC'(Sdt) + w n 0 [B-3] The last term on the L.HS. of [B-3] can be expressed either through a porous medium compressibility a', defined by d (1 - n) d a' d n iS as s 1 s 1- n do' dt - n dt or through the volume strain d n V ~ V - aE 1 s [B-4] -s -t = 1 n dt Recall that in the present work for do = 0, we have do' = d(Swp). With [B-4], [B-3] becomes: dp d p V * q + nS't w + nC'(S) dt + S a 0 [B-5 ~-r w dt w dt w t Various approximations are possible at this point. For example, (a) VVp << ap/9T, hence dwp/dt can be replaced by ap/at (b) V p << ap/t, hence d p/dt can be replaced by $p/$t. ~~~~~~Ss

-35(c) VsVSw < asw/at hence C' (Sw) can be defined by asw/Dp. With these approximations, [B-5] reduces to V q + [nS + nC'(S)] D + S 0 [B-6] r ww at w at For saturated flow, [B-6] reduces to V ~ q + nat -t P0 [B-7] where 3 = (1/Pw) Dpw/ap. Another possible form is to introduce (a) through (c) in (B-2) leading to V q + nBS P w + Sw [B-8] In [B-1 and [B- at wthus introduced the ef In [B-7c and [B-8i we have thus introduced the effect of fluid compressibility.

-36APPENDIX C List of Symbols d' = thickness of the membrane (filter) (cm) F' = body force per unit volume of porous medium (gr/cm3) I = unit tensor k ~ = medium's permeability at saturation (cm2) k = relative permeability (cm2) kr = k (cm2) o reQ 2, k' = permeability of the membrane (filter) (cm) n = porosity (cm3/cm3) no = initial porosity (cm3/cm3) 2 p = pore water pressure (gr/cm2 ) pc = capillary pressure (gr/cm2) Pc p0 = initial pore-water pressure (gr/cm2) q = specific discharge with respect to fixed coordinates (cm/sec) qrek = specific discharge relative to the moving solid (cm/sec) 1 = radial distance to internal surface of soil sample (cm) R2 = radial distance to external surface of soil sample (cm) Ro = radial distance to the surface of the liquid layer(cm) Sw = degree of saturation (cm3/cm3) t = time (sec) U = radial displacement (cm) _r V = velocity of soil grains (cm/sec) ~S r,e,z = cylindrical coordinates a = total stress tensor (gr/cm2) o' = integranular (effective) stress tensor (gr/cm2) ~ = strain tensor (cm/cm) ~ - = volume strain (cm/cm) p' =' Lame constants of the porous medium ] ~ = compressibility of soil sample (cm2gr) [=r- 1 /(2i' + 1

-37REFERENCES Alemi, M. H., Nielsen, D. R., and Biggar, J. N., 1976. Determining the Hydraulic Conductivity of Soil Cores by Centrifugation. Soil Sci. Soc. Am. Proc. 40: 212-218. American Society for Testing Materials, 1958. Procedures for Testing Soils: Am. Soc. Testing Materials, Philadelphia, Pa. Batel, W., 1961. Menge und Verhalten der Zwischenraumf lussigkeit in Kornigen Stoffen. Chemie Ing. Techn. 33: 541-547. Bear, J. and Pinder, G. F., 1978. Porous Medium Deformation in Multiphase Flow. Jour. of the Engr. Mech. Div., ASCE.: 881-894. Bingeman, J. B. and Coates, J., 1960. Centrifugal Filtration through Beds of Small Spheres. A.I.Ch.E. Journal 6: 58-62. Biot, M. A., 1941. General Theory of Three Dimensional Consolidation. J. Appl. Physics. 12: 155-164. Bishop, A. W., 1960. The Principle of Effective Stress, Publ. 32, Norwegian Geotech. Inst., Oslo, Norway: 1-5. Briggs, L. J. and Mclane, J. E., 1907. The Moisture Equivalent of Soils. U.S. Dept. Agr. Bur. Soils, Bulletin 45. Briggs, L. J., 1910. Moisture equivalent determinations and their application: Am. Soc. Agronomy Proc. 2: 138-147. Briggs, L. J. and Shantz, H. L., 1912. The Wilting Coefficient for Different Plants and its Indirect Determination. U.S. Bur. Plant Industry Bill 230: 1-82. Corey, A. T., 1977. Mechanics of Heterogeneous Fluids in Porous Media. Water Resources Publ., Fort Collins, Co. Grace, H. P., 1953. Resistance and Compressibility of Filter Cakes: Part III: Under Conditions of Centrifugal Filteration. Chem. Eng. Progr. 49: 427-436. Gray, W. G. and O'Neill, K., 1976. On the General Equations for Flow in Porous Media and Their Reduction to Darcy's law. Water Resources Research. 12: 148-154. Hassler, G. L., and Brunner, E., 1945. Measurement of Capillary Pressures in Small Core Samples. AIME Transactions. 160: 114-123.

-38Johnson, A. I., Prill, R. C., and Morris, D. A., 1963. Specific Yield - Column Drainage and Centrifuge Moisture Content,: U.S. Geol. Survey Water Supply Paper 1662-A: 1-60. Joo, G. and Lederer, P., 1974. Osszenyomod6 Lepenyu Centrifugalis Szures Szamitogepes Szimulalasa. Magyar Kemikusak Lapja. 29: 132-137. Lee, I., and Donald, I. B., 1968. Pore Pressures in Soils and Rocks in Soil Mechanics, Selected Topics. I.K. Lee, ed., American Elsevier Company, New York, N.Y.: 58-81. Marx, J. W., 1956. Determining Gravity Drainage Characteristics on the Centrifuge. Petroleum Transactions, A.I.M.E. 207, 88-91. Miller, E. E. and Miller, R. D., 1955. Theory of Capillary flow - I. Practical implications: Soil Sci. Soc. Am. Proc. 19: 267-271. Piper, A. M., 1933. Notes on the Relation between the Moisture Equivalent and the Specific Retention of Water-bearing Materials. Am. Geophys. Union Trans. 14: 481-487. Prill, R. C., 1961. Comparison of Drainage Data Obtained by the Centrifuge and Column Drainage Methods. U.S. Geological Survey Prof. Paper 424-D. Art. 434: D399-D401. Prill, R. C. and Johnson, A. I., 1963. Centrifuge Technique for Determining Time-Drainage Relations for a Natural Sand. U.S. Geological Survey Prof. Paper 450-E. Art. 236: E177-E178. Purchas, D. B., 1971. Industrial Filtration of Liquids. 2nd Ed. CRC Press, Cleveland, Ohio. Storrow, J. A., 1957. Hydroextraction: Flow in Submerged Cakes. A.I.Ch.E. Jour. 3: 528-534. Terzaghi, K., 1923. Die Berechnung des DurchlHssigkeitsziffer doe Tones aus dem Verauf der hydro-dynamischen Spannungserscheinungen. Sitzungsberichte, Wissenschaften, Abteilung 112, 132, Vienna, Austria: 125-138. Verruijt, A., 1969. Elastic Storage of Aquifers in Flow through Porous Media.R.J.M. De Wiest, Ed., Academic Press, New York, N.Y.: 331-376.

-39LIST OF FIGURES Figure 1 - Top View and Cross-Section of a Porous Medium Sample Placed in a Centrifuge (Radial Case). Figure 2 - A Horizontal One-dimensional Soil Sample in a Centrifuge. Figure 3 - Change of Variables Along a Soil Column Placed in a Centrifuge (numbers above the curves show centrifuging time in minutes).

40..cr' ~ L. -—: 0 ~,','., a, ~ I IC CL (f)I ~p e r for a't e d.,,,,;_\_.1_ L RZ~

vertical axis IL) soil porous sample plate R2 — R1. 2 11~~~~~~~

42 Pi I Sw 0 (%) 0 m 1 151 1 1.601 6 —-- - 6 0606 -7 120 50__________________2 1 2 0 Ur/0 instantaneous displa2- | cement due to centri) )~ (~c m) | fugation 40 3~~~~30 30 20 ~2 10=12 ot

SECTION II MATHEMATICAL MODEL- FOR REGIONAL LAND SUBSIDENCE DUE TO PUMPING Part 1: Integrated Aquifer Subsidence Equations Based on Vertical Displacement Only. Jacob Bear and M. Yavuz Corapcioglu2 Department of Civil Engineering, University of Michigan Ann Arbor, MI 48109 ABSTRACT A mathematical model for regional land subsidence has been developed by employing Terzaghi's concept of effective stress and assuming vertical solid compressibility only. First, the groundwater mass conservation equation in a compressible aquifer is developed for a three-dimensional space. Assuming that the problem justifies an assumption of a thin aquifer, relative to horizontal lengths of interest, this equation is then integrated over the aquifer's thickness, taking into account the conditions on top and bottom surfaces bounding the aquifer. The result is the commonly used aquifer flow equation in terms of averaged piezometric head. Introducing a relationship between changes in averaged piezometric head and corresponding changes in aquifer thickness, a single equation is obtained in terms of land subsidence. On sabbatical leave from Technion-Israel Institute of Technology, Haifa, Israel. On leave from Middle East Technical University, Ankara, Turkey. -1

-2The development is carried out for a single and a multilayered confined aquifer and for a phreatic one. An example is given of the application of the equation in radial coordinates describing land subsidence in a hypotetical six-layered aquifer.

-3INTRODUCT ION Land subsidence is often caused by the withdrawal of water from an underlying aquifer, or the withdrawal of gas, oil (and water) from an underlying reservoir. As a result of pumping, the pressure in the water (= pore pressure) decreases. This is accompanied by an increase in the integranuZar stress (= effective stress) in the solid matrix, and compaction of the solid skeleton of the aquifer ensues. The latter manifests itself in the form of land subsidence. Although this phenomenon is more significant in confined and leaky-confined aquifers, it may also take place in phreatic ones. In fact, changes in saturation, and hence in pore-pressure, in the unsaturated zone will also contribute to subsidence, though to a much lesser degree. When soft materials (sometimes as layers or lenses) comprise part of the aquifer's solid skeleton, or when aquicludes separate leaky aquifers from each other, appreciable compaction may occur. Actually, both vertical and horizontal displacements of the solid skeleton take place as a result of changes in the effective stress. Vertical displacements mainfest themselves as land subsidence. However, horizontal ground displacement, sometimes of significant and damage causing magnitudes, have also been observed.

-4Only vertical soil movement is considered in the present paper. In a separate paper, to be published as Part 2 of the present one, the authors consider both vertical and horizontal displacements. In order to determine land subsidence,we have first to determine' the (effective) stress distribution in the solid skeleton, and then employ some assumed stress-strain relationship to determine the distribution of strain. In determining the stress distribution, we rely on the relationship between the effective stress and the pore pressure. Based on the effective stress concept introduced by Terzaghi [1952], two basic approaches may be found in the literature on land subsidence. In the first one, originally presented by Biot [1941], a simultaneous solution is sought for the pressure in the water and for the strain in the solid matrix. Actually, Biot's theory describes the strain in the solid matrix in a three-dimensional space, i.e., with both vertical and horizontal displacements taking place. Verruijt [1969], who employs Biot's approach, shows that when only vertical displacement is being considered, Biot's formulation reduces to that presented by Jacob [1940] who assumed vertical consolidation only. In the second approach, following Terzaghi [1925], water pressure is first obtained by solving a simple partial differential equation (expressing water mass conservation) in terms of water pressure in the aquifer as the only dependent variable. In deriving this

-5equation it is assumed that pressure changes produce changes in the effective stress, which, in turn, produce changes in the porosity of the solid skeleton. Once the water pressure distribution has been derived, the effective stress and the resulting strain distribution is determined. Finally, the latter is used to determine land subsidence. Thus Terzaghi's theory is implemented as a two step procedure [e.g., Helm, 1975]. Pore pressure distribution is either obtained from field measurements, or calculated independently by solving the fluid flow equation. In this approach, the land subsidence is assumed to be one-dimensional (vertical) only. Gambolati and Freeze [1973] and Narasimhan and Witherspoon [1976] made use of this approach. The first approach stated above, based on Biot's fully coupled three-dimensional formulation in terms of pore pressure and displacements, has been further developed by several researchers. In a coupled three-dimensional model, one fluid-flow equation in terms of pressure or piezomettric head and three equilibrium equations in terms of vertical and horizontal displacements are the governing equations. Consequently, numerical methods for the simultaneous solution of the coupled partial differential equations have been utilized. Ghaboussi and Wilson [1973] and Lewis and Schrefler [1978] present finite element solutions. Safai and Pinder [1979] extend this approach to subsidence of saturated-unsaturated

-6porous medium. It should be emphasized again that by this approach, the lateral deformation can also be simulated. An early attempt by McCann and Wilts [1951] to arrive at a mathematical analysis of land subsidence constitutes another kind of approach. They investigated the consequences of two different models, both based on elastic continuum mechanics, called "tension center" and the "vertical pincer" model. A more general but analogous treatment was subsequently presented by Gambolati [1972] and Geerstma [1973]. In principle, the problem can always be treated as one in a three-dimensional space, determining the stresses and the corresponding strains at every point within the flow domain. The strains are then interpreted as vertical and horizontal displacements. However, in regional land subsidence problems, we are usually interested not in the variations of displacement, say in the vertical direction, along the vertical thickness of an aquifer, but rather in the integrated effect over the entire thickness of the aquifer. The integrated effect is a function of the horizontal coordinates and of time only. We have introduced here the term "regional" to emphasize that we are not interested in a "local" problem, where horizontal distances of interest are comparable in magnitude to the thickness of the considered aquifer. Instead, the considered problem of land subsidence, and horizontal displacements,

-7when these are also considered, is such that horizontal distances of interest,say between wells and points at which subsidence is determined, are much larger than the thickness of the aquifer. We then have the case where flow in the aquifer is considered essentially horizontal and the hydraulic approach is applicable [Bear, 1979]. In a similar way, we may integrate (or average) Jacob's [1940] model of vertical displacement only, or Biot's three-dimensional consolidation model, over the thickness of an aquifer and obtain a model in which the dependent variables: averaged piezometric head in the aquifer, averaged horizontal displacement, and land subsidence, are functions of plane coordinates and of time only. The two-dimensional problem is simpler than the three-dimensional one. One should note that the considerations of an aquifer system as a thin, two-dimensional slab has been reported by a number of authors [e.g., Corapcioglu and Brutsaert, 1977], but theirmodels are different from those obtained by averaging over the aquifer's thickness. Obviously the integrated approach should not be used for problems of a local nature. The present work is reported in two parts. In the present paper (Part 1), the regional land subsidence model is obtained by employing Terzaghi's concept of effective stress, and assuming vertical soil compressibility only. First, the three-dimensional equation of (saturated) water mass conservation is developed for compressible fluid and

-8and solid matrix. This equation is then integrated over the aquifer's thickness, taking into account conditions on the top and bottom surfaces bounding the aquifer. The result is a flow equation in terms of averaged piezometric head. Relating changes in head to land subsidence, a single equation is then obtained for land subsidence as a single dependent variable. In this way, land subsidence is obtained as a single step procedure. The development is carried out for a single and a multilayered confined aquifers and for a phreatic one. An example is given of the application of the equation describing land subsidence in a layered confined aquifer. In Part 2, to be presented in a subsequent paper, the starting point is Biot's theory of consolidation, but the objective is the same, namely to derive integrated equations for averaged pressure, land subsidence and this time also average horizontal displacements, all as functions of plane coordinates and time.

-9THREE-DIMENSIONAL MASS CONSERVATION EQUATIONS We start from the three-dimensional equation of mass conservation for a saturated porous medium [e.g., Bear, 1979] V ~ pq + pn) = 0, q = nV (1) where V and q are the velocity and the specific discharge of water, respectively, p is the density of the water, and n is the porosity of the porous medium. In deriving (1), we have neglected both the dispersive flux and molecular diffusion due to spatial variations in water density. In a deforming porous medium, it is the specific discharge relative to the moving solid, qr that is expressed by Darcy's law q = q - nV =- K * V* (2) ~ ~S where Vs is the velocity of the solid, * is Hubbert's [1940] piezometric head for a compressible fluid (see below)and K is the hydraulic conductivity tensor. We introduce the definition of a material derivative with respect to the flowing water d ( ) a( ) wdt - t+ v ~ V( ), V = q/n (3) and with respect to the moving solids

-10d ( ) a( ) S - - -+ V V( ) (4) dt 9t s By inserting (2) into (1) and employing (3) and (4), we obtain d p d n p V qr + n dt pn V ~ V = 0 (5) r dt dt -S Assuming that the solid's density, Ps, is constant, the equation of solid mass conservation can be expressed by V * (1 - n) V + 3(1 - n) (6) from which, we obtain 1 dsn d p V~ V = =-at - (7) - 1 - ndt dt where a' is the coefficient of matrix compressibiZity of a moving solid. Equation (7) actually serves as a definition for a'. In soil mechanics, a coefficient of soil compressibility a is often used v such that a e a a= a(l + e) v' 1 - n where e = void ratio = e~ + a (a'- a'0); e = n/(l - n), e, a' = initial values of e and a'. In developing (7), Terzaghi's [1925] concept of intergranular stress a = a' - p I (8) has been used. Where a is the total stress, a' is the intergranular stress, p is the pressure in the water (positive for compression), I is the unit second rank tensor, and it is assumed that no shear stresses exist in the water.

-11Actually, the total (macroscopic) stress a at a point in the porous medium can be shown [Bear and Finder, 1978] to be = (1 - n) < a > - n < P > I; < P (9) where the average stresses < a > in the solid and< Pw > win the water (positive for compression) are defined as intrinsic phase averages by s 1 < a > (x) (x',x) dU (x'), — S UO (U OS os) <p (= op (x,x) dU (x') (10) w ~U0 (U w ow ow where UOs and Uow are the volumes of solid and water, respectively, within a representative elementary volume around any considered point within the porous medium domain; x' is the microscopic point within the individual phases. However, in the practice of soil mechanics, Terzaghi's [1925] concept of effective stress, a', is defined by (8), i.e., - < P>w I (11) as if each solid particle is assumed to be completely surrounded by water at pressure < Pw >W. The latter produces a stressin the solid particle. However, as solid density is assumed to be constant, no strain is produced in the solid by this stress. Only the difference (1- n)[< as >s + < Pw >WI] produces solid deformation. By comparing (9) with (11), we obtain:

-12ar =)<a (1 - n < Pw I = a' - < Pw where the deformation producing effective stress is a' ~ >+ w I = (1 - n) [< s >s + < Pw >w I] This corresponds to Terzaghi's definition as employed in soil mechanics. Equation (11) will be employed in the present work. If we assume that du = 0, i.e., no change in total stress, then do' = d(pI) (12) Assuming verticaZ compressibility only, (12) is replaced by do' = dp (13) The compressibility of water (in motion) is defined by d p' = 1 wp (14) p dp By inserting (7) and (14) into (5), we obtain dp d p V * q + nV' + a' (15) r dt dt By assuming nap/3t>> q Vp and 3p/3t >> V ~ Vp, we may ~S approximate (15) by V q + (a' + nB') 0 (16) r 3t It can easily be verified that in writing (16), we have neglected the term [nB' V + a' V ]- Vp on the left hand side.

-13By assuming nap/3t >> q * Vp,and (as an approximation) defining the solid and fluid compressibilities in (7) and (14) by partial derivatives rather than by material ones (using the symbols a and B, respectively, to indicate this approximation), (1) becomes V q + S = 0; S = (1 - n)a + Bn (17) op 3t op where Sop may be interpreted as the specific storativity with respect to pressure changes (= volume of water added to storage per unit volume of porous medium per unit rise in pressure). From (1), by assuming n0o/3t >> q ~ Vo and 3p/3t >> V -Vp, -~S we obtain V q + (' + nS) at 0 (18) Thus we see how different approximations lead to different forms of the continuity equation. The relative specific discharge q in the continuity equation (18) can be written in the form qr = - K * V*= - (Vp + ogVz), * = z + dp (19) Po where k is the medium's permeability (or scaler k for an isotropic medium) and p is the dynamic viscosity of the water. Finally, (18) can be rewritten in terms of q* as V ~ q + S* (20) o 3t = (20)

-14where SO* = gp (a' + n3). Because in (20) we have V' qr and 0 not V * q, it would not be appropriate to refer to So as a 0 specific storativity (= change of water volume per unit volume of porous medium per unit change of p*), although often S* 0 is defined so in the literature. It is also often assumed that a' >> nB, so that nB is omitted. CONDITIONS ON TOP AND BOTTOM OF AQUIFER BOUNDARIES To obtain equations of flow and subsidence in an aquifer, we integrate the point (or three-dimensional) equation (18), or (20), over the aquifer:'sthickness [Bear, 1977]. By this procedure we obtain an integrated equation where the dependent variables (e.g., average pressure, or effective stress) depend only on the planar coordinates x and y and on time. The bottom and top boundary surfaces of the aquifer cease to serve as boundaries. The conditions on these boundaries become source/sink terms in the corresponding two-dimensional equations. In order to perform the integration along the vertical, we have to know the boundary conditions on the top and bottom bounding surfaces of the aquifer. We shall consider three types of top and bottom aquifer bounding surfaces: an impervious boundary, a semipervious one and a phreatic surface (with or without accretion). A more detailed discussion of the integrated aquifer equations, taking into account the conditions on the top and bottom bounding surfaces is given in Bear [1977]. To simplify the discussion, we shall

-15assume here that the shape of the lower aquifer boundary is steady (i.e., independent of time). Denoting the elevation of a point on this surface by b1 = bl(x,y), the shape of this surface may be described by the function F1 = F1 (x,y,z) Ez - bl(x,y) = 0 (21) We shall assume that the shape and position of the upper surface vary with time; the elevations of points on it are given by b2 = b2(x,y,t). The function describing this surface is 2 = F (xy,z,t) - z - b2 (xyt) = 0 (22) The thickness of the aquifer is given by B = b2 - bl = F1 - F2 (Figure la). For any moving boundary we also have dF _F dt = + u ~ VF = 0 (23) where u is the speed of displacement of the boundary. When we consider the flow of water, the condition to be satisfied at any boundary is [q - n uIu ~ VF = 0 (24) where [A]u, = A u - Al denotes a jump in A from the upper side (subscript u), the lower (aquifer) side (subscript Z) of the boundary. An impervious boundary is a material surface with respect to the water (i.e., no water passes through it). Hence

-16(q - n u)I VF {qr + n(Vs - u)} I K VF = 0 (25) If we also assume that F is a material surface with respect to solid, then (V - u)! ~ VF = O (26) Hence, (25) written for the impervious boundary F2 becomes qr * VF2 = O (27) An alternative form of boundary conditions can be obtained in terms of q from (23) and (25) r -n 2 for a moving surface q 19' VF = t (28) ~ K~ VF0= o for a stationary surface where nIl = n is the aquifer's porosity at the boundary. In a Zeaky aquifer, let F2 = 0 be the surface through which leakage takes place. If the surface is also assumed to be a material surface with respect to the solids, the boundary condition on it is VF O VF = ru VF KS * - V2IVF ~r uP 2 2 2 u 2 (29) where K is the hydraulic conductivity of the semipervious formation overlying the aquifer. When the leakage takes place through a relatively thin layer of thickness Bs and resistance a = B /KS, we may write the condition in the form

-17qrlQ VF2 S (30) where (1j - p*) is the difference in head across this layer. When a phreatic surface serves as the upper boundary of a flow domain; we have on it the condition aF VF { - N Vz VF2 + [n]u, at2 ~1P 2 V aF (31) -N + [n] 2 where qlu N = - N Vz is the rate of accretion, and -[n] = n - ni should be interpreted as effective porosity or specific yield. INTEGRATION ALONG THE THICKNESS OF A CONFINED AQUIFER As stated in the introduction, our objective is to obtain a field equation to be used for predicting subsidence as a function time and of the plane coordinates x and y only. To achieve this goal,we start by integrating (20) describing saturated flow in a three dimensional space, along the vertical thickness of a considered aquifer, taking into account the various boundary conditions on its top and bottom surfaces. This will yield equations written in terms of dependent variables which are averaged values. The latter are functions of time and the plane coordinates X and y only. Following the procedure outlined by Bear [1977, 1979, p. 522], the integration of (20) means b2(x,y,t) | (V* q + So t )dz =0 (32) bl (x,y)

-18By performing the integration, we obtain q VF2 qr VF IF F1 + S [B + ( mq _rr F2 at at 2 aF + - — F2 2* 1 = 0 (33) where the tilde (-) symbol indicates an average over the vertical thickness B, e.g., ~*(xyt) = | *(x,y,z,t)dz (34) B =(B) and B = B(x,y,t). The primed symbols indicate vector components and derivatives in the x,y directions only. For a confined aquifer by assuming verticaZ equipotentiaZs, i.e., ~* -P| *' F2 F1 and using the condition given by (27), (33) reduces to V' ~ Bq + S* B * 0 (35) o 0 t In (35), all dependent variables and coefficients are now averaged values. The average relative specific discharge, qr,may be expressed in terms of a, assuming vertical equipotentials, by q' (x,y,t) = b2,'(yz,t)dz = - K' V'F*dz qr B ~ r B K K' rb2 K' b B [b2 V'*dz B [V 2 *dz + * VF2 VF1 1] b b 2 K' B F[VB4* + P* 2 F2 V F VF 1] K B [BV'* + (*V'B +;*F VF2 - * VF1 ]-K' V'* (36)

-19where K' (x,y,z) and K' (x,y) are the hydraulic conductivity tensor and its average over the vertical, respectively. Hence, (32) becomes V~ (BK' -V'*) = S* B * (37) Let us introduce the following approximations for the spatial and temporal derivatives of %* V'v* V' V'p + V'z, % = z + p (38) 0 g go and as* as 1 2p + ~z (39) Dat t t a t go where z = b1 + B indicates the average elevation of the (moving) mid-thickness of the aquifer. By integrating (17) instead of (20), and making use of (28), we obtain DF 3F 2 1 + S __g B ~* V' ~ Bq' - n2 t - + n1l at + SOP (40) By assuming n n F1 nlF (or accurately so, if the aquifer is homogeneous), (40) becomes -V' ~ Bq' =n + S P g 3t (41) Since we have assumed that the bottom surface at z = bl(x,y) is stationary, and redefining a' in (7) for vertical compressibility only, as 1 3B 1 B 1 2 B a3-W B ~ B (2 ap ap

We obtain from (39) = (~~~~~~~1 + 1 DB (43) at 2 ~ ) at'B~p g Or, since usually'Bp.g/2 << 1 (e.g., for sandy-clay aquifer, pg = 100 N/m3r, -7 = 10 m /N, l'Bpg = 1000, B = 50 m) 3B BP, a pg''Bp' g (44) at = ( a B at at 1 + 2 Hence (41) becomes - V'*Bq = SO B (45) 0 at where S ** = pg (ns' + So). If we use (37) and approximate 0 op Sop = (1 n) a + n(3 n)a' +, then S)* = pg(a' + ni8), op o i.e., with these approximations q' qr', and both (41) and (45) may be approximated by V' * (K'BV'4) = S B (46) 0 at where So E S* S** = average specific storativity. o 0 0 If net withdrawal takes place at a rate of Qw = Qw(x,y,t) (in terms of volume of water per unit area per unit time), we should add - Qw(x,y,t) on the L.H.S. of (46). This comment is valid also for the other aquifer equations developed below. We note that in (46), the transmissivity of the aquifer T = K' B and its storativity S = S B vary with B. In fact, K and So also depend on the porosity n-which varies continuously during the consolidation process. Usually this effect is neglected in aquifers (not so in clays).

-21Equation (46), usually neglecting the affect of B on T and S, is the basic continuity equation for aquifer flow, i.e., for predicting the piezometric head distribution ( = ((x,y,t) in an aquifer. INTEGRATED EQUATION FOR A LAYERED CONFINED AQUIFER Let a confined aquifer be comprised of m layers of thickness B, B = EB [ Fig. l(b)]. Each layer is homogeneous, but m (M)m the various layers may have different hydraulic (e.g., K) and elastic (e.g., a') properties. We assume here that some of the layers are semi-pervious, or even impervious, they are not continuous, so that the whole complex may be regarded as a single aquifer with essentially horizontal flow. By integrating (20) over B, we obtain M b (Vq at ) Z m 1(V qrm + S* dz + Sdz bl bm_l, U M -' (VBmq' +JqV F q I VF m=l r rmF m rmF m-lu + S* Bt) = (47) where we assume that equipotentials are vertical throughout the thickness B(i.e. m = fl4F = 1 *IF = (* for all m); subscript u and Z denotes upper and lower side of F. If we assume that all the surfaces separating adjacent layers, as well as top and bottom impervious bounding surfaces are material surfaces with respect to the solids, then for every surface Fm, m= 1 through M- 1, we have [q r] VF = O. Hence (4) reduces to rum m'

-22M M V.mZE (B q') + ( mZ So* Bm) - (48) m=l m-rm m om m t or, using (37) M M V E( K' Bm)' V'E* = ( m S Bm) t(49) ~m=1 m om (49) Equation (49)can be used to define the transmissinty (T = K' B ) (m) -m m and storativity (S = (m) S* Bm) of a layered confined aquifer. (n) om m INTEGRATED EQUATION FOR A LEAKY-CONFINED AQUIFER For a leaky aquifer, where F2 is a semipervious boundary, by following the procedure outlined above and using on F2 the condition expressed by (30),we obtain V' Bqr' +q VF + S * B a (50)''- Bqr ru 2 0 where qrIu VF2 may be expressed by either (-Ks V4*)u VF2 or, for a semipervious layer which is a relatively thin (B') membrane, by Ac*/as where aS = BS/Ks is the resistance of the membrane and A~* is the difference in ~* across the membrane (positive for leakage from the aquifer). Inserting q' as expressed by (33) in (47), leads to V' (KB V ) * VF S* B (51) ~ u 2 o -( We may also apply to (51) the same approximations leading to (46). When pumping takes place, we add the term -Qw(x,y,t) on the L.H.S. of (51).

-23Equation (51), with transmissivity T - K'B and storativity S = S * B are the basic aquifer continuity equation used for 0 determining ~* = %* (x,y,t). Usually the effect of changes in B on T and S are neglected. INTEGRATED EQUATION FOR A PHREATIC AQUIFER To emphasize the fact that in this case the upper boundary is a phreatic surface, we shall replace b2(x,y,t) by the elevation h(x,y,t), so that the thickness of the saturated flow domain will now be given by B = h - b1. Following the procedure and assumptions (e.g., vertical equipotentials) outlined above for a confined aquifer, we first obtain (33), rewritten in the form V'.Bq + qrF2 2 rF1 VF1 + S B (52) r -F o 2t = 0 From (28) it follows that on the impervious bottom, qr F1 -rJF! VF1 = 0. For the phreatic surface, we obtain from (31) aF qr F2 YF2 t s F2 2 or qr F2 2 ( 53) -r1F2 VF2 = +Syat - (nV)F2F * VF2 (53) where S = - [nIu, is the specific yield of the phreatic aquifer. We often neglect the third terms on the R.H.S. of (53), assuming nV ~ VF << S h/~t. Then, with* c- h, (52) becomes ~ s2 y

-24V' ~ Bq - N= (S + S*B) (54) y 3t (54) t In general, S *B << S o y In this case, the total flow in the aquifer is given by h(x,y,t) Bq = f (-K'' V4*)dz - K'- (h - bl)Vh (55) mb _ 1(xy) where the approximation is due to the assumption h*P*Ih b The resulting equation is therefore V' * K'(h -b )-Vh - N = {S + S* (h - b (56) 1 Y (5 6) When the phreatic aquifer is composed of several layers, we replace in (56) K(h - b ) by K B, (m) m m and - where Bm = h - b (57) B~1)~~~ (h-(m) S *(h - B) by E S* B 0o 1 (m) om m One should note that actually within a phreatic aquifer, the total stress at each point varies as the water table fluctuates, although in our basic assumption expressed by (9), we have neglected this effect. SUBSIDENCE IN A CONFINED AQUIFER At this point, the usual procedure for determining land subsidence 6(x,y,t) in an aquifer due to pumping, is to start by determining the piezometric head distribution f= l(x,y,t) or the drawdown s = d(x,y,0) - %(x,y,t) in the aquifer, as

-25produced by the pumping distribution described by - Qw(x,y,t), using (46) for a confined aquifer, (51) for a leaky aquifer, or (56) for a phreatic aquifer (always adding - w on the left hand side). Then, the total subsidence,say for a confined aquifer 6(x,y,t) = B~(x,y) - B(x,y,t) (58) is determined by using 6(x,y,t) = B a' pg(s - Az) -BB' pgs = B 1 Ie (59) where, usually, B and e are taken as their initial values B~ and e~. One may also adjust B and n (or e) as subsidence and soil compaction progresses. Actually (59) is valid only for our present assumption of1 = bl(x,y). Otherwise, the land subsidence should be defined only as 6(x,y,t) = b2 (x,y) - b2(x,y,t). In this procedure, the effect of change in porosity on permeability and on soil compressibility, as well as the continuous change in B, is usually neglected; initial values B~, nO, k0 are usually employed. Instead, let us try to state the prob: of 6= 6(x,y,t) as the dependent variable. We shall continue to assume that b1 = bl(x,y), i.e., independent of time. We shall assume, as is common in consolidation studies, that some initial steady state exists, and that pumping produces incremental effective stresses and pressures which cause subsidence. Accordingly, with ~ replacing ~*, we write:

-26n = n~ (x,y) + (x,y,t) f n (x,y), n~ >> n = 0( (x,y) + ~e(x,y,t) - P0 - s (x,y,t) p = p~(x,y) + p (x,y,t) a' = aY'(x,y) + o'e(x,y,t) (60) B = B~(x,y) - 6(x,y,t) -e z = z~(x,y) - z (x,y,t) (60) Since B B O z - z~ = (b + ) - (b + ) = ze 12 ( 2 2 = a'o - p0 = a'O(x,y) + oae(x,y,t) - [po(x,y) -e + p (x,y,t)] = const. We have,e -e >e =ep (61) Hence'!B, He B pe t (62) since 6(x,y,t) =' B[p(x,y,t) - pO (x,y)] = a' Bp (63) Inserting (60) through (63) in (46) we obtain for the undisturbed initial steady state V' ~ (K' B ~ V'~~) = 0 (64) and for the deviation from steady state produced by pumping V' ~ (K' B e) -Qw(x,y,t) S B (65) where usuall we neglect the effect of subsidence on T-K'B in (64).

-27Obviously some steady state pumping may also be included in (64). Similarly, steady state recharge and natural replenishment may also be included. Now, ae 1 aze a a S t -- at (2) pg pg c'B 1 1 (6 1 1 1 + ( - C + +, 6 (B) (66) pg'B pgc' Pgeg'B i.e., the effect of changes in the aquifer's axis is negligible, and v~e -e -e 1 ) + oVb C z - V(+-) + V( ) Pg pg O B - ( + )V + 1 -6V( 1 ) Vd (67) pga'B'g V'B pg'B where we have assumed that &'B may vary in the xy plane, yet we have assumed 1 3B 1 3&' 1 6 1 1 1 < t ) << -- og' V B + 1 With these approximations (which amount to linearizing the effect of the gradual change in B) we obtain from (46), with pumpage K' K'Af~ = 11 69) V V6) - Q (x,yt) = S (69) pg'w pg' or, with K' = k' Pg/p, k' V ~ ( V6) - Qw(x,y,t) = $at t

-28S where = =(1 + ) Sv -- - - pga' a' In Soil Mechanics, the coefficient of consolidation for an isotropic medium is defined by Cv = k/p' = k(l + e)/pav k (1 + e~)/paa For clay fi << a'. Then S = (1 + n ) 1 and (70), with C = k'/a'P, may be written in the form a ~ 6(71) V * (Cv V6) - Qw= at In order to solve (71) for 6=6(x,y,t), with 6(x,y,O) = 0, in a given aquifer domain, we have still to provide boundary conditions on the boundaries of the flow domain in the xy plane. At a sufficiently large distance from the zone of pumping, subsidence vanishes, i.e., 6= 0. At the circumference of fully penetrating pumping wells, we may use the condition of specified pumping rates. For an isotropic aquifer, we obtain: Q)(t) 5s 9 1 -e 6 2 = rBK ar r= r rBK r (P 1 ) r=r w pg w.___ ~ rk M6 ~ ) rBK( + 1)- - =rC (72) atpgB r w l r=r w W W which is a Neuman Type boundary condition. It is of interest to note: (a) By the various approximations proposed here, we have actually linearized the problem by neglecting the effect of subsidence on Cvi although we may still have a nonhomogeneous aquifer, with C =C (x,y) ~V ~V

-29(b) Due to the linearization, superposition is permitted. This means that the total subsidence is equal to the sum of subsidences produced by the individual wells, each with its own pumping schedule. (c) We note that the thickness B is eliminated from both the p.d.e. (71) and the boundary conditions. Puzzling as this may seem at first, this is a reasonable consequence of the linearization and of the fact that for the same Qw' drawdown is inversely proportional to thickness while subsidence is proportional to drawdown. Equation (71) and the various initial and boundary conditions are similar to those commonly used for determining drawdown and piezometric head distribution in aquifers. The advantage lies in determining 6(x,y,t) in one step rather than by first determining drawdown. SUBSIDENCE IN A LAYERED CONFINED AQUIFER Let 6 denote the contribution of the mth layer to 6. m Then, with vertical equipotentials, and hence with the draw-e down s (and p ) common to all layers, we obtain the total subsidence 6 from 6 = E 6 = pg B am (s - Az m) pg( m B ct m) (m) m (m) mm m -e = ( E B )p (73) (m) or a 6 = [ Z B( e Ap (74) (m) m i With the nomenclature of the previous section applied to each layer separately, and with p* ~ we now obtain V' ( Z K' B ) 7V' ~ = O 0 (75) (m) =m m

-30for the initial steady flow without subsidence producing pumping and V ~* [ Z(K' B )VI'e] - Qw(x,y,t) = (S B) e (76) (m) m m wOm Bm) t (76) for the subsidence causing piezometric head changes produced by the pumping. By applying (66) and 67) to each of the individual layers, (76) becomes approximately *v ~ (C V6) - Qw(x,y,t) = S(77) vat where Z Km Bm Z SomBm _ (m) T (m) S v =; Sv = pg (m' B pg mE CmBm Pg o()Bm pg Z a'B (m) (m) (m) ('mm (78) where Cv and Sv may vary in the xy plane. Again, if water compressibility is negligible (i.e., nm8 << am), then S 1. It should be noted that in Cv, the contribution to T is primarily by the aquiferous layers while the contribution to ac'B is m m primarily by the clay, or aquaitard, layers (see example below). Table 1 gives some typical values of Cv EXAMPLES Let us demonstrate the one-step approach proposed here by examples of a single well pumping from a (practically) infinite homogeneous confined aquifer. For this case, (71)

-31reduces to a 2 1 ad ad (79) + - a] = a The boundary condition at a sufficiently large distance from the pumping well, say R = r is 6(R,t) = 0. At a fully penetrating pumping well, which is pumped at a constant rate of flow, Qw' the condition is obtained from (72) 36 Qw lim [r D-] = - w (80) 2 rCv r = r + 0 w where we have assumed an infinitesimally narrow well. The solution of (79), subject to the above condition, is Q 2 = w w(_ ) (81) 4 TCv 4Cvt Equation (81) is analogous to the classical Thesis solution for drawdown in a single pumping well;W( ) is the well function for a confined aquifer. As a second example, we consider a region (approximated as a circle of radius R) in which more or less uniform pumping takes place through a large number of wells. In this case we may use the solution derived by Hantush [1964 for drawdown in an aquifer which contains a large number of pumping wells within a circular area of radius R with more or less uniformly 2 and rR distributed pumping rate. For this case, for t >0.4r /Cv and rR

-32-'Q-w R2 4Ct r 2 47T e ( )[1 - exp )1(exp )} (82) v 4C t R 4Ct R4C t for t > O.4R2/ and r > R:ZQw j2 2 2 6= (w{W( r.5) + cR exp (- r ) (83) 4 rCv 4Cvt 4C t 4Cvt EQN is the total pumping within the circular area. In order to demonstrate the application of this solution, let us assume a layered confined aquifer, for which the data is summarized in Table 2. For this layered aauifer system, C is calculated to be 6 x 10 cm7sec. For EQw = 500 lt/sec, the amount of subsidence calculated at a distance of 3 km from the center of pumping by using (83), is shown in Figure 2.

-33SUMMARY AND CONCLUSIONS By averaging the three-dimensional equation of water mass conservation, over the thickness of an aquifer, and introducing a relationship between changes in averaged piezometric head and land subsidence, a single equation in terms of land subsidence has been obtained for pumping from a confined aquifer. This was stated in the introduction as the objective of this paper. The development has been based on (1) Terzaghi's concept of effective stress, (2) an assumption of essentially horizontal flow in the aquifer, and (3) an assumption of vertical aquifer compressibility only and no horizontal displacements. Similar equations were also derived for a layered aquifer. Obviously, due to the various approximate assumptions involved, the resulting land subsidence should be viewed only as an estimate of the true land subsidence. Such estimates should be sufficient for-most engineering purposes. The justification for averaging over the vertical stems from the fact that we are interested here only in regionaZland subsidence problems, i.e., problems in which horizontal lengths of interest are much smaller than the aquifer's thickness. Otherwise, we have to treat the problem as one in a threedimensional space.

-34An advantage of the resulting equation is that the problem of determining land subsidence has been reduced to one in two-dimensions (x and y) only. The examples demonstrate the simplicity of the proposed approach. In a subsequent paper, the authors present another model for regional subsidence, based on averaging Biot's three-dimensional model. Such a model provides also estimates of averaged horizontal displacements.

-35REFERENCES Bear, J., On the aquifer's integrated balance equation, Advances in Water Resour., 1, 15-23, 1977. Bear, J. Hydraulics of Groundwater, McGraw Hill, New York 1979. Bear, J., and G. Pinder, Porous medium deformation in multiphase flow, J. Engrg. Mech. Div. Amer. Soc. Civil Eng., 104 (EM4), 881-894, 1978. Biot, M. A., General theory of three-dimensional consolidation, J. Appl. Phys., 12, 155-164, 1941. Corapcioglu, M.Y., and W. Brutsaert, Viscoelastic aquifer model applied to subsidence due to pumping, Water Resour. Res., 13, 597-604, 1977. Gambolati, G., A three-dimensional model to compute land subsidence, Hydrol. Sci. Bull., 17, 219-226, 1972. Gambolati, G., and R. A. Freeze, Mathematical simulation of the subsidence of Venice, 1. Theory, Water Resour. Res., 9, 721-733, 1973. Geertsma, J., Land subsidence above compacting oil and gas reservoirs, J. Pet. Tech., 734-744, 1973. Ghaboussi, J., and Wilson, E.L., Flow of compressible fluid in porous elastic media, Int. J. Numer. Meth. in Engrg., 5, 419-442, 1973. Hantush, M.S., Hydraulics of wells, in Advances in Hydroscience, edited by V.T. Chow, pp. 281-442, Academic Press, New York, 1964. Helm, D., One-dimensional simulation of aquifer system compaction near Pixley, California, 1. Constant parameters, Water Resour. Res., 11, 465-478, 1975. Hubbert, M.K., The Theory of groundwater motion, J. Geology, 48, 785-944, 1940. Jacob, C.E., The flow of water in an elastic artesian aquifer, Trans. AGU, 21, 574-586, 1940. Lewis, R.W., and B. Schrefler, A fully coupled consolidation model of the subsidence of Venice, Water Resour. Res., 14, 223-230, 1978. McCann, G.D., and C.M. Wilts, A mathematical analysis of the subsidence in the Long Beach-San Pedro area, California, report, 119 pp., Calif. Inst. of Technol., Pasadena, 1951.

-36Narasimhan, T.N., and P.A. Witherspoon, Numerical model for land subsidence in shallow groundwater systems, Proc. Second Int. Sym. Land Subsidence, Anaheim, Int. Ass. Hydrol. Sci. Publ., 121, 133-143, 1976. Safai, N.M., and G. Pinder, Vertical and horizontal land deformation in'a desaturating porous medium, Advances in Water Resour., 2, 19-25, 1979. Terzaghi, K., Erdbaumechanik auf bodenphysikalischer Grundlage, Franz Deuticke, Vienna, 1925. Verruijt, A., Elastic Storage of aquifers, in Flow Through Porous Media, edited by R.J.M. DeWiest, pp. 331-376, Academic Press, New York, 1969.

-37TABLE 1. Range of typical values for K, a',and C for different soils. Material K(cm/sec)' I(m2/N) C (cm 2/sec)...Gravel..1 -.10|2 -10 -7 7 12 Gravel 1 - 10 1 10 -10 10 -10 -7 -9 -5 -2 9 Sand| 10 - 1 10 - 10 2 10Clay 10 - 10 - 10 - 10-5 10 - 104 -4 -2 -4 Peat 10 - 10 2 x 10 0.5 -50

TABLE 2. Data for a hypotetical aquifer system. Thickness Thickness Aquifers Aquitards (m) K(cm/sec)' (m 2/N) K(cm/sec) a' (m /N) 5 6 x 10 1 x 10 40 0.010 1 x 10-7 10 3 x 106 7 x 106 30 0.005 3 x 10-5 -6 7 1 x 10 5 x 106 -8 50 0.008 9 x 10

-39LIST OF NOTATIONS a = a coefficient of soil compressibility. v blrb2 = elevation of bottom and top aquifer bounding surfaces. B = b2 - b = aquifer thickness C = consolidation coefficient v ds( )/dt = total derivative of ( ) with respect to the moving solid. d ( )/dt = total derivative of ( ) with respect to the w moving water. e = void ratio F = 0 = equation of a surface, where F(x,y,z,t) = z - b(x,y,t). h = elevation of water table I = unit second rank tensor K = hydraulic conductivity tensor n = porosity N = rate of accretion o,e = (as superscripts) denote steady initial values and incremental, or excess, unsteady ones, causing consolidation. p = pressure in water < P >w = pressure in water w q = specific discharge (vector) of water = specific discharge (vector) relative to solids QW = rate of withdrawal r = radial coordinate R = radius of circular pumping area s = drawdown S = storativity (= SoB) So, So = pg (a +nB)

S** = pg(na' + So) = pg(a + nB) o op S = specific storativity with respect to pressure op changes [= (1 - n)a + nB] T = transmissivity S = specific yield of phreatic aquifer y u,f = (as subscripts) denote the upper and lower side of surface F(x,y,z,t) = 0 u = speed of displacement of the boundary U = solid displacement U = volume of solid os U = volume of water ow V, V = water and solid velocities, respectively us W( ) = well function x,y = horizontal cartesian coordinate z = vertical cartesian coordinate z = average elevation of the mid-thickness of the aquifer = coefficient of matrix compressibility ca~ = coefficient of matrix compressibility of a moving solid = water compressibility = land subsidence E = volume dilatation = dynamic viscosity of water p = water density Q, Q' = total and effective stresses, respectively ~s* = z +gpp) = Hubbert's potential < a > = solid stress b* = B (B) %* dz (and similarly for other variables and parameters).

-41= piezometric head [= z + p,/pg] = over a vector or an operator (V',V'.) denotes vector components or operators in the xy plane only (i. e., not in a', B' and o')

-42LIST OF FIGURES Figure 1 - Nomenclature for a single layered (a), and multilayered (b) aquifer. Figure 2 - Subsidence at 3 km from the center of pumping. LIST OF TABLES Table 1 - Range of typical values for K, a', and Cv for different soils. Table 2 - Data for a hypotetical aquifer system.

(q) L=WuJ r Z. p 1 * *'*. r IrA x q: q.z.= A' ":'.,,.: -:' wn 0 \Q \\\\\\\\ o ~i n ~$s~;K~ l b Ja;!nb::'.'.':''),u:

B= 142 m 3 2 Czone 610 cm /sec pumping EQ = 500 It/sec r= 3.0 km I \ R= 1.5 km observation j / point t (years) 0 5 10 25

SECTION III MATHEMATICAL MODEL FOR REGIONAL LAND SUBSIDENCE DUE TO PUMPING Part 2: Integrated Aquifer Subsidence Equations for Vertical and Horizontal Displacements. 1 2 Jacob Bear and M. Yavuz Corapcioglu Department of Civil Engineering, University of Michigan Ann Arbor, MI 48109 ABSTRACT A mathematical model for regional subsidence due to pumping from an aquifer is developed on the basis of Biot's work on coupled three-dimensional consolidation. Following Biot's work on three-dimensional consolidation, with coupling between mass conservation and equilibrium equations, a mathematical model for regional subsidence due to pumping from an aquifer is developed by averaging the three-dimensional model over the thickness of the aquifer and assuming conditions of plane stress. Both (vertical) land subsidence and horizontal displacements, as functions of plane coordinates and time, can be estimated by solving the model equations for a given confined or leaky-confined aquifer. An analytical solution is presented for the special case of a single well pumping from an infinite homogeneous isotropic aquifer. The solution provides O10n sabbatical leave from Technion-Israel Institute of Technology, Haifa, Israel. 2On leave from Middle East Technical University, Ankara, Turkey.

estimates of changes in averaged (over the vertical) values of piezometric head, vertical subsidence and horizontal displacement. The results indicate that under the conditions of the studied case of radial flow, the solutions for piezometric head is identical to the one obtained by non-coupled models. Furthermore, half the volume strain is produced by vertical subsidence, while the other half by the horizontal displacement. Hence, the vertical subsidence is only half the value obtained in noncoupled models which neglect horizontal displacement. A numerical example demonstrates these conclusions.

3 INTRODUCT ION The author's objectives in Part 1 [Bear and Corapcioglu, 1980] of this paper were to develop a mathematical model to determine the regional subsidence, 6 = 6(x,y,t), produced by pumping from an aquifer. As a first step, a water mass conservation equation was developed for a point in a threedimensional flow domain using Terzaghi's concept of effective stress, and taking into account the compressibility of both the water and the solid matrix. However, only verticaZ aquifer compressibility was considered. This equation was then integrated (or averaged) over the thickness of the considered aquifer, employing a procedure which takes into account the conditions on the top and bottom surfaces bounding the aquifer [Bear, 1977, 1979]. In this way, the commonly used aquifer continuity equation in terms of averaged piezometric head = -(x,y,t) was derived. Then, assuming certain relationships between changes in % and changes in aquifer thickness, B, a single equation was obtained, in which the dependent variable was land subsidence 6 = 6(x,y,t). This equation enables a direct determination of 6(x,y,t), for a given distribution of pumping from an aquifer. Certain averaged consolidation coefficients have to be known. Because only vertical aquifer compressibility was considered, often observed horizontal displacements, which

4 are due to the nonuniform distribution of drawdown, could not be taken into account in the model described obove. Accordingly, the main objective of the present paper is to develop a mathematical model for regional aquifer consolidation due to pumping, in which both vertical land subsidence as well as horizontal displacements will be considered. Obviously it is possible to achieve this objective by solving the problem as stated among others by Verruijt (1969), in a three-dimensional space, following the classic work of Biot [1941]. This, no doubt is the most accurate approach. Analytical solutions of such problems were presented by de Josselin de Jong (1963), de Leeuw (1964) and Verruijt, 1969). An example of a numerical solution of this problem, employing essentially the same formulation, but for a saturatedunsaturated domain, is presented by Safai and Pinder (1979). Numerical solutions of this problem were also presented by Ghaboussi and Wilson (1973) and Lewis and Shrefler (1978). However a basic assumption underlying both the earlier paper on this subject (Bear and Corapcioglu, 1980) and the present one is that we are dealing with a "regional" problem, i.e., one in which horizontal length of interest, say between a pumping well and a point at which both the piezometric head and land subsidence are observed, are much larger than the aquifer's thickness. In such cases, the aquifer may be considered as a relatively thin slZab, with averaged (over the vertical)

5 properties and with averaged dependent variables (e.g., piezometric head) which are functions of the plane coordinates x,y and of time only. The reduction from a problem in a three-dimensional space to one in two-dimensions only, is achieved by the procedure of integration over the vertical, already mentioned above. This procedure is also employed here. In this way, the resulting subsidence model contains averaged piezometric head (or water pressure), land subsidence and averaged horizontal displacement as dependent variables which are function of x,y, and t only. To a certain extent, the work follows that of Verrujt (1969), but w'ith a more generalized averaging approach. To avoid repetition, the basic methodology of averaging over the vertical, taking into account the conditions on the top and bottom boundary surfaces, is presented only in Part 1 and is not repeated here. Although in this work the porous medium is assumed to be perfectly elastic, other types of media, described by different stress-strain relationships may be employed, following the procedure presented here. We shall consider only the cases of a confined and a leaky confined aquifer. However, in the latter case we shall assume that the semi-previous formation is a relatively thin membrane and shall not take its compaction into account, separately.

6 THE INTEGRATED (WATER) MASS CONSERVATION EQUATION In Part 1 (Bear and Corapcioglu, 1980), we have developed the following equation for water mass conservation equation in saturated flow, where both the water and the solid skeleton are considered compressible d p d n w S pV qr + n dt + Pdt + pnV * V = 0 (1) where all symbols are defined in Part 1 (see also list of symbols at the end of this paper). With U denoting the solid's displacement and ~ denoting the volumetric strain, we have au au au au d 1 d n = V U = + + V V= V = 3x ay - - s at ~s dt 1-n dt (2) By using (2) to express V *Vs and d n/dt, and assuming a >> V -VE and ap/~t >> V * Vp, (1) becomes V r + + n P - 0; B 1 I (3) V ~qr +at at+ n = ap or, in terms of %* P V qr + + pgn =* 0~; + gdp (4) Po We now integrate (4) over B = b2 - b1, using the procedure presented by Bear (1977, 1979). We obtain

b2 (x,y,t) I CV *~r + at Pgn = t (5) q + rF2 VF2 -rF VF1 + pgnSB t B 2 ]aF= 0 (6) agj[Ft F1+ For impervious top and bottom boundaries which are also material surfaces for the solid, the second terms on the L.H.S. of (6) vanishes in view of the condition qr * VF = 0 satisfied on an impervious boundary described by F(x,y,z,t) = z - b(x,y,t) = 0, where q = - K- V~*. -r t By assuming vertical equipotentials, i.e., f* = F = *F, the last term on the L.H.S. of (6) also vanishes. 2 We obtain r- 0 V' -Bq + B - + pgnBB t (7) where Bq' is expressible by Bqr = -BK'"V' (8) In the presence of sinks of magnitude Qw (x, y,t), (7) becomes V' Bq + B + pgn6B + Q = 0 (9) V + gn t Q

8 In a leaky aquifer, q r F VF2 expresses the leakage through F2. We shall express this leakage by AC*/L where L = resistance of semipervious layer = ratio of the thickness to the hydraulic conductivity of that layer; A%* denotes the difference between the averaged piezometric head in the aquifer (4*) and that in an overlying one. The continity equation then becomes V B'r + B + Q = (9a) 3t 3t L + Qw Because we are interested only in subsidence, we separate the flow into (a) steady flow, including possibly steady pumping, without subsidence, and (b) excess flow producing subsidence. Denoting parts (a) and (b) by superscripts o and e, respectively, we have p(xy,zt) = p0(x,y,z) + pe (x,y,z,t) q (x,y,z,t) e= (x,y,z) + q (x,y,z,t) -r = r + r.*(x,y,z,t) = %*(x,y,z) + *e(x,y, zt) etc., and the corresponding averaged values p(xy,t) = PO(x,y) + pe(x,y,t)'r(x,yt) =+'e(x,y,t), etc. The continuity equation (9) can now be separated into V'Bq o + Qw = 0 (10) r ++w V'B're +B -~ + pgnB= 0+ O (11) + a giif3B a w

9 For a leaky confined aquifer, we add A~*/L on the L.H.S. of (11). In separating qr into its two parts, we have neglected the effect of consolidation on k,p, etc. Equation (11) may be rewritten in terms of pe by employing VI-~ P~ b +b V'r* vp + V' = - V'p + = z +, z= 2 (12) Pg Pg and I + (13) at =g - ~ at at The averaged temporal rate of change of s, a/at can be obtained as follows B = B VV = V' BV' + V VF- V VF (14) at -S -s s ~sF2 2 S F1F 1 Since the top and bottom surfaces are assumed to be material surfaces with respect to the solid, we have on them (VS - ) F VF = 0 (15) Hence: Vs F VF = U VF = - aF/at (16) and (14) reduces to Bt- = V BV'- + aB +

10 With Vs defined by (2), we obtain Da u a3F a at at F at F BV' = j V'dz =dz BU) + U' f t 2: - F2 2 (18) At this point, following the integration procedure, we need information on U' IF and U. This information is not available. One could use an assumption with respect to slip or to shear stresses on the boundaries. Instead, we introduce a simplified assumption, namely that of practically no variation in horizontal displacements along the vertical. U U IF UF (19) ~ ~ ~F2 F1 In view of (18) and (19), we now obtain: aB v' B + B' U B =-t -t B'* U' +(20) The continuity equation (11) can therefore be rewritten for a confined aquifer in terms of variables C* (or pe) and U' (or components Ux and Uy ) as xy V' ~ Bqre + V' - B t —+ at + pgnBB + = 0 (21) ~r t Qw where Bq' can be expressed in terms of o*e following (8). For a leaky confined aquifer, we add At*/L on the L.H.S. of (21).

11 Equation (21) may be linearized by introducing B(x,y,t) = b2(x,y,t)- bl(x,y,t) = b (x,y) UIF -(b (x,y) + UF - (by) + U 2 1 = B~(x,y) + AZ; A = UF - UIF (22) and neglecting second order small terms. THE INTEGRATED EQUILIBRIUM EQUATIONS The componentsof the total stress tensor, a, at a point within'the flow domain satisfy the following equilibrium relationships (Biot, 1941 Verruijt, 1969). x + xy + Xz + f =0 "x 3y 3z x yx Yy yzf =(23) + + a + f = (23) zx +zy + f = 0 3x 3y 3z z where f represents the body force and the inertial force a2Ui components p Ui, i = x,y,z, have been neglected. Using at2 ='- pI to express the total stress in terms of effective stress and pressure (positive for compression) in the water, and (following Verruijt, 1969) separating both a' and p into steady initial values ao' and p0, and consolidation producing incremental effective stress,'e, and excess pressure pe, we obtain for the initial steady state equations

12 1j + f DpO ~7+ f?- _ = 0,i,j = x,y,z (24) 9x. i ax. J i and for the incremental effective stress and pressure Wa'e x, aP 0 = i,j = x,yz (25) ] 1 The summation convention is to be invoked in both (24) and (25) as in subsequent equations written in identical notation. In writing the last two equations we have assumed (as a good approximation) that the body force remains unchanged, i.e., f = o, = x,y,z. We now assume that the solid is isotropic and (for the relatively small deviations considered here) perfectly elastic. For such a solid, the stress-strain relationship are expressed by: aU. aU. Du a'e = G( + ) + (Xk)ij (26) ij!x x ax 1D i J k where E = aUk/aXk, G = E/2(1 + v) is the shear modulus, E is the modulus of elasticity, v is Poisson's ratio and X = E/(1 + v)(1 - 2v); X and G are the Lame constants. Relationships other then (26) may also be used. Following the methodology of deriving integrated aquifer equations, before continuing to integrate (25) and (26), we have to introduce the conditions on the top and bottom boundary surfaces.

13 The conditions for total stress at any boundary F = 0 is [a] ~ VF = O (27) [ u,' or, in terms of effective stress and pressure [C' - pI] u, ~ VF = 0 (28) Separating into initial steady state and displacement producing excess effective stress and pressure, we obtain for the upper surface (F2 = 0) [oa + a } VF2 = {(o'O - poI) + (ate _ peI)} * VF2 where, as before u and Z denote the upper (or external) and lower br internal) side of F. Hence, C~ u VF2 =(,O0 - pOI) ~ VF2 eu~ VF2 = (ale - peI)L * VF2 Assuming now that the total stress remains unchanged, i.e, =e = 0, we obtain the condition (le _P I) F VF2 = 0 (30) In a similar way we obtain for the lower surface (F1 = 0) (ae pei) VF1 = 0 (31) Next, we integrate the first equation in (25) i.e., for i= x

14 I a'e ale ae b2 + r' b x ay + x Z -aX ) dz = 0 (32) DX ay Dz Dx bb 1 a B ( ( ) + ae F 2 e jF, ax e + y'el y ax Bxx f xx y xy F2 y xx F2 ~x xx F1 2x aF + p F 1 = (33) However, from the first equation of (30) and (31), we obtain: al |~Fi aXFl aDF. a atF. a-e...I i I xx F1. CY y'e + o-'e = 0 i=1t2 i ax xyFz F. ~X (34) Since aFi/xz = 1 for i = 1,2, by inserting (34) into (33), we obtain a (Ba5le) + a(Bea) e a (Be = 0 (35) In a similar way, from the two remaining equations of (25) we obtain a(BYx + a (Be) a (e = (36) (B er) + a (Bple ) = 0 (37) ax Xz axy yz in which all averaged values are functions of xy and t, only.

15 We now express the average excess stress components in terms of averaged displacements, making use of (26) and the assumption expressed by (19). We obtain: au au au au au A FE ++ E +. = + - + + + x y z ax ay z ax +y B au au A ale X- + 2G~ = (X + 2G) + X( Y + xx x DX ay B _cr_ + x yy y ay ay B au au G =e G( Y+ X) - xy Dx ay yx oe = G( + 2G = G -y + [U + u =z z a aZ ax B x x + G axB + xzy +n nz + te G eo B iF1 ZF1 xlz G(f t G n zof x, u..u ay = G a - Bp-+= G (3 Gy a a Z- F yz + [ F x'U aU DU A al X x= Y) + (X + 20) (38) zz 3Z ax ay B By inserting these expressions in (35) through (37), we -e obtain 3 equations in the 4 variables p UxI U and u, all functions of x,y and t. D aux 3u Aza u au {B[(X + 2G) + X( ++ BG( ax Dx ay B ay ay ax n_ (m~e) = 0 (39) r%_

16 au au aU aU A + y { B[A+ (A. + 2G) -Y + 3 X ay 9 x' y B a (Bpe) = 0 (40) u DB F 2 u 1 BG ax + G (Uz Fx z F2 x Uz F 1 au DB 9F 2F + — { BG yZ + G U) = 0 (41) 9y ay z 9yD Y z F2 ay z F1 a- = For constant X and G (actually Xand G), and with B(x,y,t) = B(x,y) + Az(xyt), Z << Bo, we obtain by linearizing (39) and (40) (A /Bo) ae GV' Ux + (A + G) - G x (42) 9'A /Bo) -e GV' U + (A + G) - G = 0 (43) or U 9U X (A /Bo) De GV' Ux + (A + G)(a + ay + ax ax GV' U + (A + G)( + +)+ — 0 (45) y Dx + y ) y. y The second and third terms on the L.H.S. of (21) may also be written in a linearized form a +B (46) VI B~+ (46)

17 where B = B + A A << B and we have assumed V' ~ V'B << 3B/at. 0 z o s In (41) we have used U F = U F and UF = U F as in (19) We note that in (41) we have only Uz. In principZe, together with (21) we have 4 equations for the four dependent variables. However, in these four equations, we have UzIF and Uz F (and Az = Uz F - UF, B = BO + Az). These are actually conditions on the two boundaries F1 = 0 and F2 0, for which we have 2 no information. In fact, the land subsidence U IF is the very z 2 unknown for which a solution is sought in most subsidence problems. At this point we may continue by introducing certain simplifying assumptions instead of the missing information. For example, we may assume that the bottom of the aquifer is stationary i.e. Uz = 0 and that IU varies linearly with z, i.e., U = U z F! z z 2 z F2 - 6/2 where 6 is the land subsidence (positive downward). Then we end up with 4 equations for Ux' Uy~ 6 and pe(or P*e) Another approach is to assume that consolidation occurs under the condition of plane incremental totaZ stress, as suggested by Verruijt [1969, p. 347]. This means e e e e e e o s = a = 0, CT =- = 0 (47) zz xz zx yz zy As indicated by Verruijt [1969], this assumption is justified when the aquifer is between two soft confining layers (e.g., clay) which cannot resist shear stress. Furthermore, this assumption also justifies (19) since in a relatively thin

18 aquifer as implied by the plane stress assumption, lateral deformations is, more or less, uniform throughout the relatively small thickness. From (47), it follows that (25) reduces to,DO e _ - -a_ o-0, i = x,y (48) Dx. ax. ] 1 with the boundary conditions (30) and (31) also written the xy coordinates only. Following the integration procedure which led above to (35) through (37), we now obtain (35) and (36) only. Equation (37) drops out. We could have obtained the same result from (35) through (37) by assuming that deformation occurs under conditions of plane incremental averaged total stress described by (47). Accordingly, (41) vanishes, but (39) and (40) or (42) and (43) or (44) and (45) remain unchanged. We now have to solve (21), (44) and (45) for pe,, U and A. The needed x y z fourth equation is now obtained from the first condition in (47), which leads to and e f (49) and therefore, from the expression for oe in (38), au au az A pe X (-x + ) + (X + 2G) = + 2G (50) ax ay ~~~B

19 This completes the formulation of the problem where -e p (x,y,t), U (x,y,t), U (x,y,t) and A:(x,y,t) are the sought x y z for unknowns. Usually we assume UzlF = 0, hence A = UIF 1 2 = - 6(x,y,t) = land subsidence. Initially, the values of these variables is zero. As external boundary conditions, we shall usually assume vanishing values of all these variables at a distance sufficiently remote from the zone of pumping (Q (x,yt)). By differentiating (42) with respect to x, (43) with respect to y, linearizing them and then adding the two equations, we obtain 2 X 2 (X + 2G) V' ( + y + Vz 2 V e or, for constant A, G and B~ V2 [(I + 2G) (3X + y ) + X B =0 (51) Following Verruijt (1969), we integrate (51), obtaining 3U 3U A A x( z z =e (a +2G)( + y X= (A + 2G)6 - 2G B p + f (x,y,t) (52) -2 Where f satisfies V f = 0 for every t, comparing (52) and (50), obtained by introducing the plane stress assumption, we find that au au A A f = 2G( X 3y 2G(6 - 2 z) (53) where i defined by B(38). where E is defined by (38).

20 If we assume no average lateral displacements, i.e., ir' = 0, (21) reduces to V'B4qe + + pgnSB t BB + A (54) and from (50) we obtain be Az p= ( + 2G) B (55) which are two equations in p and a z. We maov combine (54) and (55) to yield 1 e V:Bq' + 1 2Ge+ = 0 (55a)'Bqr + B[+ 2G + n] t where we have assumed A << B. z For a leaky-confined aquifer, we add + 8~*/L on the I.H.S. of (54). Equations (54), (55) and (55a) are obtainable also by assuming vertical compressibility only, with cl' = 1/(X + 2G), as assumed by the authors in Part 1 of this paper. If we now compare (53) with U' = 0 with (50), obtained by the horizontal plane stress assumption and vertical displacement only, we obtain A f = - 2G B As pointed by Verruijt (1969) the function f "describes the deviation of the simplified Terzaghi-Jacob Theory from the Biot Theory" where the former assumes vertical consolidation only, while the latter considers the three-dimensional nature of consolidation. Here f expresses the deviation

21 when integrated aquifer consolidation equations are used. It is of interest to note that here in the comparison between (52) and (55), f does not vanish, while Verruijt (1969) shows that in the comparison between Biot's (1941) three dimensional considerations and (55), f vanishes. INTEGRATED EQUATIONS IN RADIAL COORDINATES When we consider consolidation in the vicinity of a single pumping well, we start by writing all flow and equilibrium equations in cylindrical coordinates, assume radial symmetry and integrate over the aquifer's thickness. Equation (21) reduces to 1 D oae 1 D r B Yh e ~-r trB (q (rB + - +pgn B t 0 (56) r ar rBr r ar at at at For a leaky confined aquifer, we add + AP*/L on the L.H.S. of (56). Under the conditions of plane incremental total stress (i.e., ar = re = e 0), the only remaining zr'z zz equilibrium equation reduces to D'e o e _Cle r r r 0 ad ape + rr 0 (57) r r Tr The constitutive equations reduce to au aDU J aU d'e - 2G r + r + r z rr 3r Dr r Dz U ae - 2G z + X zz 2G + XE (58) zz +a

22 By averaging (57) over the vertical thickness of aquifer, using (30) and (31), we obtain rr F e r,e -e) a (Be) =0 (59) ar r rr ee0 ar Averaging the constitutive equation yields au u A au r r z r ore (2G= (2G ] + X r (60) rr -r r B ar U au A U r+[ r z,)e += (2G + ) r(61) 688 r B r a U A ~ Az r +_r= ] z,e = (2G + X)- + - + -] 2G + x (62) zz B ar r B where ~ = 3Ur/3z + Ur/r + Az/B. By substituting (60) and (61) into (59), we obtain r B au r _BU 3 2G~(B 0 (63) r [B(2G r- + X%)] + - 2G( (B ) = (63) rr B r r ar r 3r From (62) and the assumption of plane incremental total stress, we obtain oa' = 0 - &ze pe. Hence (to be compared with (50) above). A au u A ~e z ( r z pe + (2G r+ r+ ) = X( + 2G (64) "B ar r B Equations (56), (63) and (64) constitute now a set of equations to be solved for Ur, Az an d p By inserting Ur = 0 in (64) we obtain the relationship for the assumption of vertical compressibility only. By linearizing (63) we obtain

23 au U -U 2G r + - = (65) r Dr r ar -r r Integration of (65) with constant G and X yields _U U A r r "e _ ~e 2G(-+ ) + x+ pe (2G +,)~ -2G - p = g(t) (66) r r B where g(t) is the arbitrary function of t. Equation (66) has also been derived by Verruijt (1969) without the vertical integration. However, for the particular set of assumptions underlying the analysis here (plane stress assumptions and vertical integration), f in (52) and g in (66) are identical. A comparison of A au U = G + - r r p = 2G + = (2G + X)e - 2G[ —-+ ] (67) ~B Dr r (67~r obtained from the plane stress conditions (49) and (62), with (52) leads to au U A A r r zz f(r,t) = 2G(- + r- -) = 2G(C - 2 -) (68) Dr r BO B= which is similar to (53). The term Az/B, representing the boundary conditions on F1 = 0 and F2 = 0, is introduced by the vertical integration procedure. It is of interest to note that the term au au V'U' - (X+ Y) is replaced in the radial case by 1 X (rU + V' U', where in the latter case U' = U lr. r r r r r -r r

24 By inserting (64) into (63) we obtain 2G( + r zG ) - 2 - 2 ) = g(t) 2Gr r B B From (68) it then follows again that here g(t) E f(r,t). In comparison to Verruijt's (1969) non-integrated equations, au we observe that in Verruijt's case f(r,t) - g(t) = 2G z' he difference is due to the fact that he does not introduce the plane stress assumption in the equations written in casterian coordinates. AVERAGED EQUATIONS FOR A LAYERED CONFINED AQUIFER Figure lb gives the nomenclature for a layered confined aquifer, comprised of M layers of thickness B B = B m (m, m Each layer is homogeneous, but the various layers may have different hydraulic (e.g. K) and elastic (e.g., G,X) properties. We shall assume that equipotentialZs are vertical throughout the aquifer. This means that when some of the layers are semipervious or even impervious, they are discontinuous, so that the entire complex behaves as a single aquifer, with essentially horizontal flow.

25 One way of solving such a problem is to treat it as one in a three-dimensional space, employing appropriate boundary conditions at the surfaces separating the individual layers from each other, as well as on the top and bottom bounding surfaces. Obviously, this approach is the most accurate one. Another possibility, which yields estimates of subsidence for problems of a regional nature, is to apply the averaging procedure employed above, to the individual layers. However, that will require information about conditions at the interfaces between the layers and also about the rate of pumping from each layer. The layers are coupled by the transport of water and stress across these inner boundaries. To circumvent the lack of such information, we shall consider the entire aquifer as a single unit with averaged behavior. By integrating (4) over B, we obtain b M b I 2 q(V +q + pgn t- )dz =(V bq rm (br m=l bm u 1 m-lu + + gndz= VB + q IF VFm q m-VF +~t m m Dm=l m rm rm I Fm m-l -+ ( grBmt (70) + Bm at + ( Pm)B = 0

26 where we have assumed = for all m; m-1 m subscripts u and Q indicate upper and lower sides, respectively. If we assume that all surfaces separating adjacent layers, as well as top and bottom impervious bounding surfaces, are material surface Fm, m = 1 through M - 1, [qr ]u, VFm = 0, where [A]uk = AIU - A[, indicates the jump in A from the upper side of Fm to its lower side. On the top and bottom impervious surfaces we have qrFM VFM 0 and rF VF 0, respectively. Hence, the second and third terms in the sum vanish. The contribution to the total horizontal flow EBm q m ~rm will come mainly from the aquiferous layers. By following the procedure outlined in (14) through (20), we obtain for Bm Cm/*tassuming that within each layer Ul - Um I U m- 1 m m m - m VLB B +t V B -V'' + m (71) m -- =t m m t Um (t If we now assume that all averaged horizontal displacements are equal, i.e., Um = U' for all m, then M m a -U' B B. -= -VB + B - V'U' U+ 3B (72)

27 Altogether (70) becomes: V' Z 3U'B + B V'-U' + - (m) (Bmqrm) at t at + nmBm at (73) to be compared with (21). For a leaky confined aquifer, we add + A/*/L on the L.H.S. of (203), expressing the leakage qrFM VFM. Equations (23) through (26) are also valid within each layer. With the assumption that the total stress remains unchanged, i.e., ae - 0, (30) and (31) are valid for the top and bottom bounding surface. We may therefore rewrite them for FM and Fo, (Ole peI) VF = 0; (o'e _ pel) - VF = 0 (74) M o On all intermediate surfaces, we have the condition expressed by the second equation of (29), namely (e -peI) IQ. VF = ('e peI) VF (75) With these conditions, let us now integrate the first equation in (25) i.e., for m = x M b 3c, e mae Cre A i me, xx + ( xy + )dz = 0 (76):ml (~x Yy +d =( m-l u + e | aFm aFm + Mxy F - - _ F { (B aee ) +,ei m,u 2y xv P Iy x x

28 ax mm elF ax e (77)'x (B Pm) - Fm, m + IF amx 0 (77) m,9. M-lu In view of (74) and (75), (77) is reduced to Z ae) (B +e ) _( m 0 (78) { Bax n xxm (Bm xym ax (Bmm m=l and in a similar way, from the remaining equations of (25), we obtain M ~,e -e {- (B -Be.)}(B -0 (79) ax M yxin (Bmv M -(B - m=l M a{Y(Bme) + e Bm)} = 0 (80) m=l yzm in which all averaged values are functions of x,y, and t. Equations (38) are valid for each individual layer. Hence, for the layered aquifer, the equilibrium equation (39) becomes Mau aU A au 1 [B m( + 2Gm) ax + BX m(_y m m [B x m=l ax mm m 9x mm vy B ~e au asB p + y-) ] - m = 0 (81) _y ax Equations similar to (40) and (41) may be derived in a similar manner. At this point we introduce the simplifing assumption (consistent with the plane stress assumption to be invoked below) that

29 U = U and U = U for all m; the total land subsidence is xm x ym -Y given by ZA = 6. Under this assumption, (81) may be approxi(m) mated by r3ff D2G) 9U DU A DU DU CB-( + 2G+ X + [B{G(+X +k ax [B( + 2G) X ay B +ayLLa +x ay (Bpe) = 0 (82) where BX = EBmXM BG = EB mG and we have assumed that >> Zi i.e., changes in averaged piezometric head are much larger than vertical displacements of the mid-elevation of the individual layers. All averages in (82) are over the entire thickness of the aquifer. Thus, with the above assumptions and definition, (40) remains valid, if X and G are replaced by X and G, respectively. Attempting to continue in the same manner and develop the equation corresponding to (41), will end up with the same difficulties as encountered above when dealing with a single layers. We therefore introduce at this point the assumption of plane incremental total stress. This means that (47) may be written for each individual layer. Equations (37) and (41) vanish. Equation (50) with X and G replaced X and G and remains also valid. The problem actually reduced to one in a single layer with averaged properties. EXAMPLE: DISPLACEMENTS DUE TO PUMPING FROM A SINGLE WELL Consider a single fully penetrating well of radius r w pumping at a constant rate e = Qw from an infinite single layered confined aquifer. We wish to estimate the subsidence 6(r,t), the

3Q ~e horizontal displacement Ur(r,t) and the pressure drop -p (r,t). Initially, the aquifer's thickness is B~(x,y) = constant. The coefficients k,X and G are assumed constant. (a) The continuity equation (56), which with e k 3e aB ~ Az r r B( -r)r- B P B t(), and (r B B (r(r —-) 1 iAz < 1B aB B (i.'e.,<< 3t and - << ), reduces to 1 a k 2p 1 aT(r b + +a t: 0 (83) aSr ) at Actually, assuming a hydrostatic pressure distribution along B, e ae e ar Dr Dr where bl + b2 z = 2 Here and in (88) below we have assumed e >> P r 2r a

31 With k~, ino, fi denoting the respective initial values or values at some intermediate pressure, we obtain the linearized equation 1 p + e - (83a) r- (.r n- + + n (83a) (b) The definition of dilatation aU U A r + r + z (84) (c) Equation (66) rewritten as z -e (2G + XA) - 2G P + 2g(t) (85) where the 2 was added for convenience. (d) Equation (64) p B~ 2+ %E (86) BO There are 4 equations in the variables p, E, Ur and Az The boundary and initial conditions are: t < O r > r Pe U rz' (87) t<O r>wr pr z Q~O W 3r t>O r=r w (88)w 2~rr BO'k0

32 r Xp, = O (89) z r = r, U = 0 (90) r + oo U 0 (91) Following Verruijt (1969), we combine (85) and (86) to yield (A + G)~ = p + g(t), g(O) = O (92) By inserting ~ from (87) in (83), we obtain a(r k + - - 1 eg(t) 0 (93) r (r )r g - (X + G t x - + G t or 2ee 1 e e g(t) where a" = l/(G + A) and C = = = T p(a + i) (a" It+ =f) where S is the aquifer's storativity and T is its transmissivity. Verruijt (1969) already indicated that (94), (except for the last term on the L.H.S., which will be shown below to vanish as g(t) = 0), is identical to the commonly used aquifer continuity equation, derived by assuming vertical compressibility only. However, in (94) the aquifer's compressibility is defined by a" rather than by'= 1/(X + 2G) (compared with (55)). The Laplace transform of (94) leads to +... = -1 gs (95)

33 where an overscope denotes the Laplace transform of a variable. Equation (95) is a Bessel equation. Its general solution is = C1Io(r/~Cv/sI + C2KO (r/VC /s) - (96) k where I and K are zero order modified Bessel functions of o o0 first and second kinds, respectively. From (88) and (89), written in their transformed forms = 0 at r + and w obtain 3r r=r 2-rr Bok, we obtain Qo, /C /s r = (97) 2 2=rr wBOOS K(rw/ /C /s) w 1 v Actually, from (85) and (89) it already follows that g(t) - 0. Hence eQw V~C v/s K (r/ /Cv/S) Pe = pe ~ ~ J~~~.(98) 2Trrw BOkK 1 (rw/ s) Making use of the approximation x K1(x) + 1 for x + 0, be and taking the inverse of the resulting expression p, we obtain - ewWu r2 Sr 2 e p_ = W(u), u (99) 47TT 4C t 4Tt which is the usual equation used for drawdown in a confined aquifer. For a leaky aquifer, Verruijt (1969) obtained the same result except that W(u) is replaced by the Well Function for a leaky confined aquifer. In writing (99), p(p) was taken a constant as for as it affects T. Also, we have neglected the change in z (see (13)) in the relationship between e and p.

34 By summing (85) and (86), we obtain -e P = (X + G)~ = ~/e" (100) Hence: Qw~pga ~= - W4~TT W(u) (101) If we assume o" >> nB, we may neglect fluid compressibility, and (101) may be rewritten as: Qw 2 -0 E - 4CBOW(u), t Cv c (102) 4 C B 4C t (v 1,o 4 Equation (85) and (86) may also be combined to yield au U A r r z - (103) 3r r B 2 which is a very interesting result. Under the assumption of vertical compressibility only, we obtain (55) while here, the relationship between p and ~ is given by (100). We note also the difference in the coefficient as indicated already by Verruijt (1964). Accordingly the subsidence 6 = - A is given by Q W(u) (104) =8'TCW(u) (104) v This is half the value obtained in Part 1 where it was assumed that only vertical consolidation takes place. The latter assumption leads to an overestimation of subsidence.

35 From (103) and W(u) =- Ei(-u), we obtain Q 2 w r r-Yrr 8 a -C B~W ( 4Ct) (105) rUr = 8C BEi( 4C t)dr + C3 V V Qwt C3; r rU = 4'~B| Ei(-x)dx + C3; x = 4C t where W(u) = - Ei(-u) is the exponentiaZ integraZ. Using a table of integrals, we obtain Q t wQ -x rUr 4rrB- [xEi(-x) - 1 e + C3; or -Qwt - C Ur 4Br [uW(u) + 1 - e ] + ---- (106) r 47rBOr r Using the boundary condition (90), we obtain Q t -u r2 C3 4 [UwW(uw) + 1 e ], =Ct 4TrB~ w w w 4C t V Hence Q t -- u — U At this point we approximate, assuming uwW(uw) 0 and e w 1. Then -Q t -u U 4 r [uW(u) - e + 1] r 47 B Or

36 or Qwr 1 e-u Ur 16c B~[W(u) + 1 u (107) As r +-o, u -+, Ur - 0, thus satisfing (91). Since Ur = 0 at both r = r and r =, it may be of interest w to find its maximum value at r = rm. We calculate r 1 r (rUr) _ - [W W(u) 1 e D- rDr r r r 8-C BO 2W(u) - --— 2u For DUr/r = 0, we have 1 - e-u W(u) = -, u = 0.323 r = 1.1367 JT- (108) m v (108 Figure 2 gives the plots of 4~T.e 4rC B~ 8TC 16-C B~ 4r e v - v6 and v U Q c' Qwr r wQw Q Qw Qwr r Given T, B, Qwt C v a", r and t, we calculate u = r /4Cvt and determine e,, 6 and Ur from these curves. For a leaky aquifer, the analysis is very similar to that given above for a confined aquifer. In (83) we add + A4/L on the left hand side. This case was solved by Verruijt (1969). He obtains for p an expression similar to (99) except that the Well Function for a confined aquifer, W(u), is replaced by the

37 Well Function for a leaky aquifer, W(u,r/ATI-Y). /T-L is the leakage factor. This means that the same replacement should be introduced also in (101), and (104). For the horizontal displacement Ur, we obtain 2 u r Q~r 1 - e 4TLu Ur 16'rC B0 [W(u,r//TL) + 1 e (109) v r u +4 Numerical Example (see Table 2 in Part 1 for T and Cv). Data: Confined aquifer with = 3 T = KmBm= 821 m /d 95 cm /s, B~ = 142 m, w = 180 m /h 4 3 3 2 - 5 x 10 cm /s, C = 6x 10 cm2/s; For this case, r/U = 51093. rmax e Figure 3 shows p 6 and Ur as functions of r after 3 years of pumping. ~e Figure 4 shows variation of p, and Ur with time at 3 km form the pumping well.

38 COUPLING VS. UNCOUPLING In Part 1 of this paper [Bear and Corapcioglu, 1980], the-problems of determining the pressure distribution and that of determining the distribution of land (vertical) subsidence, were actually uncoupled. Horizontal displacements were neglected. The uncoupling was introduced by the -coefficient' defined by ds 1dn dp n t - -p- a t'. Then, once the p-distribution is dt 1-n dt dt known, the subsidence could be determined b a'B p /Dt = - ~/at. w 2 For the radial case example we obtained 6 = W(t) and ~e Qwg 2 4rrC rC P =w W(r) P 4- W(4rt) In the present paper, we left the two problems coupled d c by leaving t — as an approximation of dt flow equation (3). However, in the radial case example we note that since the special conditions of the problem led to g(t) - 0 in (92), we could relate E to p and replace:/3at in (83) by ~e e,G,=73t e a /at= ( /p) (a p e/at) =(X + G-p/at 7= a" p/at In this way (83) would have become a single equation in p identical to the one solved in Part 1 and indeed, the solutions for pe are identical in both solutions. Then, once pe is known, we can use (92) to determine ~ and from it use (103) to determine Az and U Thus, the interesting point is that also in this case, under conditions which lead to g(t) = 0 in (92) the two problems may be uncoupled.

39 SUMMARY AND CONCLUSIONS The three-dimensional coupled Biot equations are averaged over the vertical thickness of an aquifer, assuming that the problem justifies an averaged approach, an assumption of shear free boundaries and plane stress. The development is based on (1) Terzaghi's concept of effective stress, (2) an assumption of essentially horizontal flow in the aquifer, and (3) elastic stress-strain relations. The resulting model consists of a set of averaged mass conservation and equilibrium equations, with averaged pressure, volumetric strain, and vertical and horizontal displacements as dependent variables for which a solution is sought. Closed form analytical solutions are derived, by using the Laplace transformation, for an example of radial flow to a single pumping well in an infinite homogeneous aquifer. The expressions obtained for averaged changes in piezometric head, land subsidence, and horizontal displacement are functions of x,y and t only. The solution for changes in piezometric head (or pressure) obtained for this example in Part 1 [Bear and Corapcioglu, 1980], where the flow and equilibrium equations are uncoupled by introducing a coefficient of aquifer compressibility and the solution obtained here, where flow and strains are coupled, are identical. In fact, this result is obtained due to the boundary conditions of the problem which leads to g(t) = 0. In other words, these conditions cause uncoupling. Therefore, by (10~, the contributions

40 to total volume strain by the horizontal and vertical strains are equal. This leads to values of land subsidence which are half of those obtained by the uncoupled approach presented in Part 1. Also an expression is derived here for the horizontal displacement with a universal function R(u) for the single well case. The results show that the horizontal displacement has a maximum value some distance, rm from the well. An expression is given for calculating this distance. When numerical solution techniques are employed, the graduate change in B, n, k etc can be taken into account Whenthe equations are linearized, superposition is applicable. For example, the total effect of wells pumping simultaneously can be calculated simply by adding the effect of the individual wells. In this way, the total excess pressure drops (or draw-downs), the land subsidence, and the horizontal displacement at an observation point in a multiwell system can be easily estimated. In order to arrive at an analytical solution for the radial case example, we have linearized the equations in the present work. In this work, perfectly elastic stress-strain relations have been assumed. The results of many uniaxial tests indicate an approximately linear void ratio versus effective stress relation on a semi-log scale, and in many cases the effect of hysteresis due to swelling and reloading is neglected. This indicates that compaction of soils is not elastic in nature. This is especially

41 exhibited by various clay and sand types. When rebounds occur different coefficients of compressibility may have to be introduced for the rebounding process which are much smaller than those for compaction. The land subsidence occuring during pumping is mostly permanent and cannot be recovered by termination of pumping or by injection. Whenever a leaky aquifer has been considered here, the consolidation of the semipervious layer has not been taken into account. We have considered only the consolidation of the pumped aquifer, assuming that it is made up of a deformable material, or that its average behavior is that of a deformable material due to layers and lenses of soft materials The general model presented here can be used, using numerical solution techniques for nonhomogenous aquifers and for other than elastic stress-strain relationships.

42 REFERENCES Bear, J., On the aquifer's integrated balance equation, Advances in Water Resour., 1, 15-23, 1977. Bear, J., Hydraulics of Groundwater, McGraw Hill, New York, 1979. Bear, J., and M.Y. Corapcioglu, Mathematical model for regional land subsidence due to pumping, Part 1: Integrated aquifer subsidence equations based on vertical displacement only, submitted to Water Resources Research, 1980. Biot, M. A., General theory of three-dimensional consolidation, J. Appl. Phys., 12, 155-164, 1941. Churchill, R.V., Operational mathematics, Second ed., McGrawHill, New York, 1958. De Josselin de Jong, G., Consolidatie in drie dimensies, L.G.M.Mededelingen, 7, 57-73, 1963. De Leeuw, E. H., Consolidatie in drie dimensies, L.G.M.Mededelingen, 8, 17-48, 1964. Ghaboussi, J., and Wilson, E.L., Flow of compressible fluid in porous elastic media, Int. J. Numer. Meth. in Engrg., 5, 419-442, 1973. Lewis, R. W., and B. Schrefler, A fully coupled consolidation model of the subsidence of Venice, Water Resour. Res., 14, 223-230, 1978. Safai, N.M., and G. Pinder, Vertical and horizontal land deformation in a desaturating porous medium, Advances in Water Resour., 2, 19-25, 1979. Verruijt, A., Elastic storage of aquifers, in Flow Through Porous Media, edited by R.J.M. DeWiest, pp. 331-376, Academic Press, New York, 1969. Verruijt, A., Private communication of a summary on Horizontal Displacements in Pumped Aquifers, published in EOS, Trans. AGU, 51 (4), 284, 1970.

43 LIST OF NOTATIONS bl'b2 = elevation of bottom and top aquifer bounding surfaces B = b 2 - b1 = aquifer thickness C = consolidation coefficient sv d ( )/dt = total derivative of ( ) with respect to the moving solid d w( )/dt = total derivative of ( ) with respect to the moving water F=0 = equation of a surface, where F(x,y,z,t) = z - b(x,y,t). g(t) = a function of time G = shear modulus I = unit second rank tensor I Ko K = Bessel functions k = permeability K = hydraulic conductivity tensor L = resistance of semipervious confining layer n = porosity o,e = (as superscripts) denote steady initial values and incremental, or excess, unsteady ones, causing consolidation p = pressure in water q = specific discharge (vector) of water gqr = specific discharge (vector) relative to solids Qw = rate of withdrawal r,O = cylindrical coordinates R(u) = a universal function s = Laplace transform variable S = storativity v? = transmissivity u,i = (as subscripts) denote the upper and, lower side of surface F(x,y,z,t) = 0

44 u = speed of displacement of the boundary U = solid displacement Ur = radial displacement,VsY = water and solid velocities, respectively W( ) = well function x,y = horizontal cartesian coordinate z = vertical cartesian coordinate a' = coefficient of matrix compressibility [( + 2G)-1-'" = coefficient of matrix compressibility [(X + G)] A = B - B z 0 = water compressibility = land subsidence 6.. = Kronocker delta 13 ~ = volume dilatation X = a Lame constant 1P -= dynamic viscosity of water p = density of water a, Q' Z = total and effective stresses, respectively = Z + gpdp = Hubert's potential -= B (B) * dz (and similarly for other variables and parameters). = piezometric head [= z + p/pg] = over a vector or an operator (V',V'.) denotes vector components or operators in the xy plane only (i.e., not in a', 3' and a')

45 LIST OF FIGURES Figure 1 - Nomenclature for a single layered (a), and multilayered (b) aquifer. Figure 2 - Universal functions for solutions in a radial flow to a single pumping well. Figure 3 - Spatial change of average piezometric head, subsidence, and horizontal displcaement after 3 years of continuous pumping. Figure 4 - Temporal change of average piezometric head, subsidence, and horizontal displacement at 3 km from the pumping well.

46. F2.~~~~~~~~~f JIB 2.,|: xyt#-q e c0,Fl=O 10~~~~~~~ _-F~~ M m =M~~~~:'.., ~~~~~~~~aurn (a W r n- M Mo4.....v.-..',.... -,.-,:.. F -- bm x yt) m=3 2 o o,.*,,,ar,, I m, t*''-',......' *@".. m3 = _rO~~ ~ (b)~ ~ r ~~ F r r + F2 (b)

10 I. I W(u) for u>10, R(u)-/u R(u) 41TT T e for u\OD1, R(u) —W(u)+l Qwr 8 TT CvK R (u)=W(u)+-(1 - eU Qw,,TT C Bo' for fe and Suse W(u) curve L. for ir use R(u) curve 0.01 0.0001 0.001 0.01 0.1 1 10 100 u- r2/(4 Cvt)

25 Q=501t/sec T=95cm/sec e/ 10 w n3 2 3 C= 6%'1O cm /sec t=3years (cm) B0= 142m 8564m 15~ 0 j'm 10 15 2(k r km)

- I V I (year!)2 0 10~~~~~~~~a1 50O (cm) r = 3.0 km

50 COMMENT After submitting the paper for publication, we found by private communication with Prof. Arnold Verruijt of Delft Technological University, The Netherlands, that in 1970 he studied the problem of vertical subsidence and horizontal displacement in steady flow produced by pumping from a leaky aquifer. It is interesting that he also concluded that in a pumped aquifer, horizontal displacements of the same order of magnitude as vertical ones may occur (Verruijt, 1970).

SECTION IV A MATHEMATICAL MODEL FOR CONSOLIDATION IN A THERMOELASTIC AQUIFER DUE TO HOT WATER INJECTION'r PUMPING J. Bear and M. Y. Corapcioglu Department of Civil Engineering, University of Michigan Ann Arbor, MI 48109 ABSTRACT A mathematical model is developed for the areal distribution of fluid pressure, temperature, land subsidence and horizontal displacements, due to hot water injection into a thermoelastic confined and leaky aquifers. The underlying assumption is that the aquifer is thin relative to horizontal distances of interest, and hence all dependent variables of interest are average (over the thickness) values. The solid matrix is assumed to be thermoelastic. Following the development of three-dimensional conservation of mass and energy equations, and equilibrium equations in terms of horizontal and vertical displacements, the mathematical model is derived by averaging the three-dimensional model over the vertical thickness of the aquifer, subject to conditions of plane total stress. The effects of viscous dissipation and compressible work have been included in the formulation. The O10n sabbatical leave from Technion-Israel Institute of Technology, Haifa, Israel. O20n leave from Middle East Technical University, Ankara, Turkey.

resulting averaged coupled equations are in terms of porewater pressure, temperature, and vertical and horizontal displacements which are functions of x,y, and t only. The equations are non-linear and have to be solved simultaneously due to the coupling that exists among them. Equations and appropriate boundary conditions in radial coordinate have also been presented for an example of a single injecting (or pumping) well.

3 INTRODUCT ION There has recently been a growing interest in two groups of problems: land subsidence and heat (and mass) transfer in aquifers, the latter connection with geothermal energy production and with proposals to store energy in aquifers by hot water injection. Only a relatively small number of papers have treated the problem of land subsidence and horizontal displacements due to pumping from a geothermal field or due to the injection of hot water into an aquifer. This paper is an attempt to contribute to this important, and interesting topic. Its objective is, therefore, to present a model for estimating land subsidence and horizontal displacement with non-isothermal flow of a compressible fluid in a deforming aquifer A rather large number of papers on heat and mass transfer in porous media can be found in the literature, especially in that of reservoir engineering. Bear (1972) reviews some of these papers. In recent years, with the growing interest in geothermal energy, and with the development of powerfull numerical techniques, renewed interest in the subject has led to the development of a variety of models simulating the flow of heat and water in geothermal (hot water dominated and mixed watersteam) reservoirs under various assumed conditions. A study, reviewing and comparing the various balance equations appearing in earlier works was presented by Corapcioglu

4 and Karahanoglu (1979). Mercer and Faust (1979) present a descriptive review of numerical geothermal models. In most works, the mathematical model consists of a fluid mass conservation equation and an energy balance equation, both for a (single, or multiphase) compressible fluid and deformable medium [e.g., Huyakorn and Pinder, 1978; Horne and O'Sullivan, 1978], but do not take into account the complete effect of compaction and consolidation. In a large portion of the studies, porosity is assumed to be a function of pressure only. The pressure and temperature related equilibrium equations are not included in the model. For example, Faust and Mercer [1979] present a three-dimensional model with coupled flow and energy equations for a single and multiphase fluids, and the derived two-dimensional equations by averaging over the reservoir's thickness. The reservoir's compressibility is taken into account through changes in porosity produced by changes in pressure. Lippman et al [1976] introduced a model which simulates the effects of geothermal production, as well as reinjection, on the deformation of a liquid dominated geothermal system. Their approach combines a numerical model for themass and energy equations with the nume~rical solution of Terzaghi's consolidation equation. Since the temperature field varies much more slowly than the pressure, much smaller time steps have been taken in the flow cycles than in the energy cycles. The full coupling of the behaviors of the porous solid matrix and the

5 fluid in geothermal reservoirs, including momentum and energy transfer and the dependence of porosity and permeability upon- fluid and solid stresses,are discussed by Brownwell et al [1977]. The ground surface subsidence history of Wairakei has been examined by Pritchett et al. [1976] in'terms of the calculated two-phase fluid flow and the local geology. A similar study for a geopressured dissolved methane reservoir is presented by Garg et al [1977]. Aktan and Farouq Ali [1978] study the thermal stresses induced by water injection. They consider a single phase compressible fluid, and calculate the stress distributions in the reservoir. Similarly, Ertekin [1978] developed a two-dimensional, two-phase fluid flow, a threedimensional heat flow, and a two-dimensional displacement model to simulate the subsidence-compaction phenomena in a,) hot water flooded oil reservoir. In a more recent study, Derski and Kowalski [1979] derived a set of linear equations based on irreversible thermodynamics. The local temperatures of the fluid and solid components have been assumed identical. Coupled equations of conservation of mass, conservation of thermal energy and equilibrium have been given for a thermoelastic porous medium saturated with a compressible thermal fluid. The assumption that both the solid and the fluid are at the same (macroscopic) temperature serves as a basis for most of the studies mentioned here.

6 As demonstrated by some of the works described above, the problem of non-isothermal flow of a compressible fluid (or fluids in multiphase flow) in a deformable porous medium can, in principle be stated as a model consisting of the following equations: - mass conservation equation and an appropriate motion equation for the fluid - energy conservation equation for the porous medium - equilibrium equations for the porous medium - equations of state for the fluid and stress-strain relationships for the porous medium. However, in many cases, in treating problems of flow and subsidence in aquifers, horizontal distances of interest are much larger than the thickness of a considered aquifer. Under such conditions, a somewhat simpler approach although this approach is an approximate one, yet yields estimates which are sufficiently good for most engineering purposes. According to this approach, the three-dimensional model is reduced to a two-dimensional one, in the horizontal plane, by integrat-. ing the former over the thickness of the aquifer CBear, 1977, 1979). In the resulting model, all dependent variables (and coefficients) are averages over the aquifer's thickness.

7 Bear and Corapcioglu (1980a and b) demonstrate the application of the approach to regional subsidence resulting from pumping under isothermal conditions. In the present paper, this approach is extended to the case of non-isothermal flow conditions, It is obvious that the approximate approach employed here should not be used for cases of a strictly local, or three-dimensional nature. We shall refer mainly to a liquid flow in a confined aquifer, with comments relating to the extension of the results to a leaky aquifers, and to multiphase flows. A large number of assumptions are introduced along the development (e.g., that the solid is a thermoelastic one), in order to simplify the presentations. Also cross-transport phenomena, like the Soret effect or the Dufour effect have not been considered in this study. However, the methodology presented here can also be applied to other sets of perhaps less restrictive assumptions.

8 THE FLUID MASS CONSERVATION EQUATION The discussion in this section is at the macroscopic level. Without any special notation, all variables are at the macroscopic level (i.e. averaged over an REV of a porous medium). Our starting point is the fluid's (macroscopic) mass conservation equation for saturated flow in a three-dimensional space V.pfq + = 0, q = nVf(1) where we have neglected mass dispersion due to microscopic variations in fluid density and microscopic velocity and where both the fluid and the medium are compressible. In (1), pf = fluid's density, q = specific discharge, Vf = fluid's velocity, n = medium's porosity and t = time. Equation (1) may also be rewritten as 1 dw(Pfn) (2) V-V + -- = 0 (2) of pfn dt where dw( )/dt = ~( )/Dt + VfV( ). Because eventually we are interested in expressing the fluid's motion by Darcy's law, we rewrite (1) in terms of the specific discharge qr relative to the moving solids d n d p n ds n dwf +r l n r V dt r -n (3)q

9 where V is the solid's velocity and ds( )/dt = -( )/at + V * V( ). The change in porosity is related to the volume strain (dilatation) s of the porous medium by d ~ d n s 1 s V (4) dt 1 - n dt's obtained from the solid's (macroscopic) mass conservation equation V'p (1 - n)VS + a P(l - n) = 0 (5) in which we insert our assumption of constant solid density ps, and the definition of ~ (see below). In the present study, Pf = P f(p,T), where p is the fluid's pressure and T is its (absolute) temperature. Hence, wPf apf dT (6) dt aP jT dt. T p dt Introducing the coefficient of fluid's compressibility (at constant temperature>,B and fluid's coefficient of thermal (volumetric) expansion (at constant pressure-) BT, defined by 1 Pf 1 3Pf P Pf T Pf ~Tp=const. (7) equation (3) can be written as

10 d ~ d p d T VCr + n n (8) r dt p dt T dt In general,p = B (T). P P The last term on the L.H.S. of (8) introduces a coupling between the continuity equation (8) and the fluid's temperature T. We shall see below that a further coupling is introduced by assuming that the dilatation e is also temperature dependent, in addition to its dependence on effective stress (i.e., c = ~(alj, T); see below). Note that a fluid sink is not included in (8) and in other continuity equations in this section. It will be introduced at a later stage. Let us approximate (8) by assuming: V ~V<< /t, V << ap/at, V sVT << aT/3t (9) Then, (8) reduces to -r at + np t T t T T For vertical consolidation only, we may assume c = s(p,T). Then d a d P d T d P d T P+ T dt T ()dt dt ap T dt - T p dt p dt T dt

11 and (8) becomes: T Voqt + (aT + nBT) + (a - +(a V + nV ~* VP p p3t T T Dt P-S + (OTVs - nBTV VT = 0 (12) With the approximations of (9), (12) reduces to Vrq + (a + n ) + (aT - nBT) T - q *VT = 0 (13) By introducing (9) and the assumption of e = (p,T), with the coefficients ap and aT, we have actually removed p the coupling with F. Equation (13) includes only p and T as dependent variables. In (10) and (13) we express q by Darcy's law: qr (VTp + p gVz) (14) where k is the medium's permeability which depends on n p and p are the fluid's temperature and pressure dependent density and dynamic viscosity. We need expressions for k = k(n) Pf = Pf(p,T) and p = p(p,T). Usually the dependence of p on p is mtch smaller than on T, and may be neglected.

12 THE ENERGY CONSERVATION EQUATION In this section, we will develop the macroscopic energy conservation equation. We start from microscopic considerations and derive the macroscopic ones by averaging the former over an R.E.V. of the porous medium. All symbols with asterisk (*) denote microscopic values. The microscopic equations of energy balance for a point in a fluid continuum may be written as (Bird et al, 1960, p. 315) aT* PC(t+ Vf *VT) = - VJT* -* f - f *Vf) f f (15) where Cs is the heat capacity (per unit volume) of the fluid at constant volume and Jf is the flux of thermal energy by conduction. The last term of the R.H.S. of (15) expresses the irreversible rate of internal energy increase per unit fluid volume by viscous dissipation. In many cases, this term is neglected. We may do so also here, or note that if we assume no shear stresses to be present in the fluid, then -(T*.VV*f) can be replaced by +(pfVV;) as an expression for viscous dissipation. The second term on the right hand side of (15) is the reversible rate of internal energy increase per unit volume by compression. We recall that d p* d E * 3 V' = -* w (16) dt dt -

13 where ~* is the fluid dilatation. f From (7) it follows that (ap/T) fP* appearing in (15) is equal to the ratio BT/8P. This ratio depends on the nature of the specific fluid and is obtainable from the appropriate pressure-volume-temperature relationships. For a certain range of temperatures, the following linearized relationship is often used with constant Bp and BT P = p*[l - T(T* - To) + S (P*- P* )] (17) where p* is the density at P* and T*. o o fo' Combined with the mass conservation equation apf t + V.p*V* = (18) 3t f-f equation (15) without the dissipation term becomes ~ d C a3 (p*C T) + V.p*C V** = - V.J* - T* ( )V- VV- T*p* w v (pTf v)f —ff -f f T f dT* f f f dt (19) where the last term on the R.H.S. drops if we assume C = const. We note that p~*C T could be replaced by pfh*, where h* is the fluid's enthalpy. For an isotropic thermoelastic solid continuum, the corresponding energy balance equation is given by (Nowacki, 1975, p. 12) aT* DE* PSC (t+ V YVT*) = - VJ* T*y (20) where C is the specific heat of the solid (per unit volume) at constant strain, and

14 ~Yij (0y* / */aTS) i = consti s = stress in the solid. ij s13 cs We note that a thermoelastic body is one in which we have a reversible elastic process and an irreversible thermodynamic one. The external energy sources or sinks have been omitted in all the balance equations presented above. They will be introduced at a later stage. Equation (20) can also be combined with the solid's mass balance equation, leading to a ~ s t (p*C T*) + V (p *C V*T*) V.J* - T*y (21) t (s 5 s s E-S s ~ s s at where we have assumed C = const. We note that the terms related to energy change by dilatation are often overlooked in the energy balance equations for both solid and fluid. In order to obtain *the macroscopic energy balance equations, we have to average (19), assuming dwCV/dt = 0 and (21), over a Representative Elementary Volume (REV) of the porous medium (Fig. 1). The methodology of averaging is well known and will not be discussed here [see, for example, Bear, 1979; Hassanizadeh and Gray: 1979 ]. Use is made of the following averaging rules for an intensive property Ga of the a phase.

15 < t > = - G > U J( )GcUa v dS (22) < VG > = V < Ga > + U1J (S- -- a (23) where <a fu (U0)GadU < Ga > = 1 ( G dU < G > =8 < G > ac a a 8 = U /Uo, S is the total area of contact between the a a o o c phase and all other phases in UO, va is the outward unit vector normal on S and ua is the velocity of S. Here u u e - n, S - S. Subscripts f and s denote the fluid and the a a fs solid phases, respectively. By applying (22) and (23) to (19), we obtain: D ~1 ( < *CT* >- D *CfT *u f dS + V' < p*C V*T* > t f V f Ui J(SSf)M f f-Uf -f f f f + I p*C T*V*v = -V.V < J* >- -vdS sf ~f U0 J(Ssf) f Tvf-Vf vf = V -f UO J(S f )J v dS T - < Tf p*V'V > (24) where C is assumed constant. If we do not neglect the dissipation term, - < T*-VV* > has to be added on the R.H.S. of (24). -f -f

16 Let us now make the following assumptions and approximations: f f f (1) That < p*C T* > = n< p*CvT* > = n < CT > < >, i.e., f v f f V f V f i neglecting the average of the product of their deviations from the average, (i.e., assuming < p C > << p >f< CTf >) (2) That the solid is a material surface, i.e., on it (V* - u). = 0. f (3) That the thermal dispersive flux of heat through the tortuous fluid filled void space is much smaller than the average convective flux through it. Then < pC vT*V* > = n< p*C T*V* >f f — f f f't f f n T<p C T> >f< V >f n <p* >f< C T* >f< V >f. The neglected f v f -f v f f thermal dispersive heat flux through the fluid is the average of the product of deviations of both temperature and fluid velocity from their respective averages (= <(C Tf*)f >). Eventually, we may include this flux in the macroscopic equation, if so desired. (4) That the conductive average flux through the fluid < J* > = < X* VT* >, where X* is the thermal conductivity of -ff f f the fluid, may be expressed by 1 f < J* > X> < T* l < TVT* > > + - Tv> dS] f f f ftV<Tf U (S of)f n - T -7 < T* > E - nX V < T* >f, f=f f f f where the second rank symmetric tensor Tf is the tortuosity of fluid occupied void space [e.g., Whitaker, 1967]. The coefficient X is the thermal conductivity of the fluid occupied void space.

17 Underlying the introduction of a tortuosity is the assumption that it is a geometrical property, independent of the heat flow pattern. (5) That < T ( n) VT* >V<T > n < V* > p f fp f =- <f> V-y> < > < >- P*> <ff p f p f (6) That < T~'VVf > = < p V-*V > f f f n < p* > < v.v > < P* > < V-Vy> f f - <p* > V-n < Vf > With these assumptions and approximations, and with a dissipation term included, (24) can be rewritten as >f M*><Vf fV<T> f an < P < C > + n < > < CT >< V-nXf V<T> 1 f T f f f o(S )Jfvf S + [<T> T < =0 > ]*= 0 0 s) ff f ~ff f pf (25) Brownwell et al [1977] relate the viscous dissipation term to the (macroscopic) fluid velocity relative to the moving solid. To obtain the averaged equation for the solid, we apply (22) and (23) to (21). We obtain

18 < P*C T* LJ -) T *Uv dS + V- < P*C V*T* > + t P C VsT sdS Vs < J* > J s - < T*y.. >- (26) We now introduce assumptions and approximations which are similar to those introduced above for the fluid phase: 01~ ~ (S s s (1) < p* C r* > << < p* > < C T* > S C S S s s s 5"C ~T"s 5 < p*C T* > " (1 - n)< p* >.< C T* > (2) (V - u)*v = 0. (3) The thermal dispersivity flux throuc the solid skeleton is negligible with respect to the convective flux, i.e., < (C T*)V* > S << < p*C T* >s< V* >, hence F sS SS -S < p* T*V* > =(1 - n) < p*C TV* >V (1- n) < p*C T* >s< V S ~ S-S s ~ S-S S 5 S -S (1 - n) < p* >< < V* >s s ~ s -S (4) < J* > (V < T* >, with similar comments about the approximation involved in assuming a constant solid skeleton tortuosity. (5) < T'y V > (1 - n)< T'* > Y> (1 - n) < T* > Y s s s't

19 By introducing these expressions into (26), we obtain (1~ - n)< p* > < C T* > + V (1 - n)<ps* > < C T* > < V* > at ( 1 - n s S $ S -S - V-(1 - n)X -V < T* > + J dS ~s s U0(S (S) *dS a < ~* > s S + (1 - n)< T* > (27) s y t We obtain the thermal energy conservation equation for the porous medium as a whole, by adding (25) and (27), noting that = - J = J*-V on S and assuming f ~, sf ~f -S -S sf' < T* T > = T, i.e., the fluid, the solid and the porous 5 f medium as a whole have the same averaged temperatures. Removing now the intrinsic average symbol i.e., < p* > Ps, I S < * > pf < V* >f V < < ~* > = ~, the resultff~~~~ -,ff' -~s ~s s ing equation is t I(npfC + (1 - n)ps C)T] + V [(npfCVVf + (1 - n)psC~Vs)T - V [nXf + (1 - n)X ]VT + [T( P]V nV + (1 - n)Ty O (28) We note that the second term on the L.H.S. of (28) may also be written as [- npfC (V - Vf) + (npfCV + (1 - n)pSC )VS]T; n(V s-Vf)qr Since we have assumed that C and C are constants, we may combine (28) with the mass conservation equations

Dnpf a(1 - n)Ps + VfV = 0; + V-(1 - n)pV = 0 (29) at f f at ~-S We obtain (P T V A VVT+ [(pC) V + P C q VT mat Mrn f T r BT E: + [T - P]V-nVf + (1 - n)Tyt= O 0 (30) where we have introduced (PC)m = nPfCv + (1 - n)psC = heat capacity per unit volume of the porous medium as a whole, and Am = nX + (1 - n)X = coefficient of thermal conductivity for the porous medium as a whole. We note that (PC)m depends on p and n, which in turn depend on pressure and temperature; A depends on n which also depends on pressure and temperature. Equation (30) may further be simplified by (T s aT (1) Noting that (pC)m(at + V VT) - (pC)md- (PC)m at where we have V *VT << aT/at. (2) Assuming that n>VV >V.Vn and therefore VnV - V (q + nV ) V-q + nV-V Vq + n nf = r -(qr -S qr at where we introduced the assumption V sVe << ~c//t, We then obtain UT ST (PC) V T p q T + pT - P].Vq at fvqr + + {n[T -- - P] + (1 - n)Ty} - = 0 (31)

With the approximations of (9), (31) reduces, for vertical consolidation only to {(PC)m + aT[n(T - P) + (1 - n)Ty]} t +V -VT p T T = + pfC VqrVT + [T - P]Vqr + a [n(T - P) + (1- n)yT] O pf~rr + 3 + [ Bp p p p (32) Equation (32) coupled with (13) are the governing equations for the temperature and pressure fields in a porous medium with vertical compressibility only. Obviously, other simplifications and approximations are also possible, depending on an estimate of the order of magnitude of the various terms. In all these equations, we have a clear coupling between pressure, temperature and dilatation (or changes in porosity). The last two terms in (30) and in (31) express the source of heat due to the internal energy increase per unit volume of porous medium by viscous dissipation and by compression. They may be relatively small in many cases of practical interest. THE EQUILIBRIUM EQUATIONS All considerations in this section are at the macroscopic level. Bear and Pinder (1978) develop the macroscopic equilibrium equations by volume averaging the microscopic ones. We assume here that the solid skeleton of the aquifer is a thermoelastic body. For such a body, the mechanical and thermal state at a given instant is completely described by the distributions of temperature, T, and strain, s... Two 1j

22 processes occur when changes in temperature take place: a reversible elastic process and an irreversible thermodynamic one. The main difference between the problem of thermoconsolidation considered here and that of consolidation under isothermal conditions lies in the stress-strain relationships which in the former case include also the effect of temperature. Neglecting the inertial terms, the total stress a at a point within the aquifer satisfies the following equilibrium equations a + + f. = 0, i = 1, 2, 3 (33) xj where f represents the body force, and the summation convention is employed. Using aj. = a! - F-i to express the total stress in 13 13 13 terms of the effective stress o!. and the pressure p (positive for compression), and separating both a!. and p into initial 1J! O steady values a.. and p and consolidation producing incremental 13 eo ao~e al,e ones a!( and p, we replace (33) by (a) + fO _ apO O, (b) 1 () 3xj i x ax. -— x = 1 J 1 where f? - f.. 1 1 For a thermoelastit porous medium, the stress-strain relationships are given by the Duhamel-Neumann relations (Nowacki, 1975) ~,~= C.. e _ 3 5) 13 ijk9 k - Bij where e where k are components of the incremental strain, T= T -e

is the incremental temperature, and the material coefficients'e ~C..' e e ) I, = ( /DT)I Cijk - ij kc T = const. are calculated in the isothermal state. The development of (35) is based on the assumptions (1) that strains are small, (2) that Te/T << 1, and (3) that the free energy depends on both the strain and the temperature. For the sake of simplicity, we shall henceforth limit the discussion to an elastically isotropic porous medium. We shall also assume that the medium is also isotropic with respect to permeability. For an isotropic elastic body (35) reduces to cje9 = 2GE. + (XE + yT e), E =k + E + (36) 13 1 kk xx yy zz where G and X are the Lame' constants of the porous matrix and y6. = (./Tt. The coefficient y is also 13 1'/ E = const. related to the coefficient of volumetric thermal expansion aT by 1 2 T P 2 = (A + - G)aT, = t (37) P + - G p Using the usual relationships between s.. and the components of displacement U. au. au. 2( + x3) ~n = V.U (38) ij 2X. s x. ] 1

24 we may rewrite (36) in the form au. aU. au e-=G(j 1 + + - YTe]6ij (39) 3J — ++k 1( Equation (38) together with (34b) into which we insert (39) provides four equations in the six variables Ui, e Te and C. The mass conservation equation and the thermal energy conservation equation provide two additional equations in terms of p, T and ~. Altogether we have six equations for the six variables. We have not mentioned qr as a variable as it is easily related to p. In addition, we need information on k = k(n), with n related to e, p = v(p,T)= p(T) and pf = pf(p,T). CONDITIONS ON THE TOP AND BOTTOM SURFACES BOUNDING THE AQUIFER According to the methodology of averaging (or integrating) over the vertical thickness, B, of an aquifer, we have to know the conditions that exist on the top and bottom surfaces bounding the considered aquifer. Let F = F(x,y,z,t) describe the shape and position a surface which bounds the aquifer from above or below (Fig. 2). The unit vector normal to this surface is given by ln = VF/IVFI. With b = b(x,y,t) denoting the elevation of this surface above some datum level, we have F = F(x,y,z,t) = z - b(x,y,t) = 0 (40)

25 For this surface dF -F + u.VF = 0 (41) where u denotes the velocity at which this surface is displaced. Water Flow Boundary Conditions Bear (1977, 1979, 1980), starting from the general boundary conditions on F for both the moving fluid and solid [q - nu]. VF = 0, [(1 - n)(V - U)] u, VF = 0 (42) where [A]u AIU - Al. indicates the jump in A from the upper to the lower side of F, shows that (a) On an impervious boundary which is also a material surface with respect to the solids, the condition (on the aquifer side) is q *VF = 0 (43) (b) For a semipervious boundary qr aquifer side.VF = qr outer side VF (44) where the R.H.S. gives the rate of leakage through F. Heat Flow Boundary Conditions For the flow of heat by convection and conduction, assuming Tf = T = T, the general condition on F is: (a) Through the fluid phase (per unit area of porous medium) [pfC T(q - nu) - nXf VT] uVF = 0 (45)

26 (b) Through the solid phase (per unit area of porous medium) [p C T(1 - n)(V - u) - (1 - n)X -VT] - VF = 0 (46) Together we have for the porous medium as a whole: [p C T(q nu) + C T(1- n)(V - u) A *VT] VF = 0 [f T n)+s E ~ s ~ m u,z (47) For a boundary which is a material surface with respect to the solid, (Vs - u) VF = 0 for both sides of F. Hence the second term in the square brackets vanishes and the first term reduces top pfCTq'VF which vanishes in view of (43) for an impervious boundary. If the upper side of F is assumed to be impervious to fluid (n = 0), then X X+ X and (47) reduces 5 s to - A VTIZVF = - X VTI vF (48) where subscript Q denotes the aquifer side of F. If the solid on the external side of the boundary is assumed to be impervious to fluid and an insulator with respect to thermal conduction, (48) reduces to A -VTI.VF = 0 (49) If the boundary, say the upper one, is not an insulator, and heat can leak through it, we have to use (48). We may replace the R.H.S. of (48) by some estimate of the rate of heat loss by conduction.

27 Equation (47) may also be rewritten as [p C Tqr + (PC) T(V - u) - Am VT] VF = 0 (50) fV r nm s ~m u, or [PfC Tqr + (pC) TV A'VT] ] VF f m ~s - ~m u, R [(pC)mT] uVF = - [(pC)mT] at (51) For an impervious insulator boundary, (51) reduces to {P fCTqr + (PC)TV A -VT} FF = - (PC) T Fat (52) f r m s m where subscript F denotes the aquifer's side. Stress Boundary Conditions On any boundary F = 0, equilibrium requires that the total stress a satisfies [J]U *VF =, l *VF = aItzVF (53) Separating the stress in to an initial steady one and an incremental, consolidation producing one, we may write [~]uj VF = 0 (54) e I] VF = 0 (55) Assuming that a remains unchanged on the external (or outside the aquifer) side of F, i.e., 0 then a VF = 0, and the condition on F becomes

28 aC|F *V]. (g'e Ipe) FVF = 0 (56) where I is the unit tensor and subscript F denote the aquifer side of the boundary surface F = 0. THE INTEGRATED EQUATIONS The dependent variables in the fluid mass conservation equation (10), energy conservation equation (31) and equilibrium equations (34b), are all functions of x,y,z and t. By integrating (or averaging) these equations over the vertical thickness of the aquifer, these variables, as well as the various aquifer properties, will be replaced by averaged dependent variable and properties which are functions of the plane coordinates x,y and of t only. With the nomenclature presented above in the discussion on the conditions on the top and bottom surfaces, the integration is carried out employing the following typical rules (Bear, 1977, 1979). (a) For any vector or J = J(x,y,z,t) and bl = bl(x,y,t), BV'J Jdz = VBJ + JI. F - JIF VF1 (57) b1 F2 2 1 (b) For and scalar 4 = i(x,y,z,t) b aF aF B-~-~ 2~4~ D2 1 B - I = dz = (B$) + F - at at at F 1 at (58) and =t2 d=VBV2 1

29 where the average is defined by' T b ~(xyt) = b 2 B (x,y,t)dz, B =b2 b1 = F1 - F2 (60) and the prime (') symbol indicates a vector or a vector operations in the xy plane only. J F1 etc. are conditions on the surfaces F1 and F2 which bound the aquifer. They have been discussed above. If we assume (58) and (59) reduce to aWC aWi v = Vp at at t The Integrated Fluid Mass Conservation Equation By integrating (10) over B, we obtain b2 (x,y,t) 2 (Vq + na at - ~qVT)dz = 0 (61) br +t T at Tqr bl(xy,t) Since we have assumed here that both F1 and F2 are impervious to water and material surfaces to solids, we obtain in view of (43), for the first term in (61) BV~q = V'.B + rF 2 q (62) qr'F2 VF -r F 1 r

30 For the second term, Bear and Corapcioglu (1980b) showed that B at VB at+ at (63) where U' (components U,Uy) is the vector of averaged horizontal xy displacement. In developing (63), it is assumed that U' U F U iF1 and V.VcE << ~ /3t F2 F1 ~ s For the third term in (61), we obtain { 7~B2E+ PB (64) np at pat PI F2 atp F1 t(64) Assuming that the pressure distribution in the aquifer is hydrostatic, then _,pgB'+ gB IF2 2' P 2 (65) Hence nf3 _ap n B(+ +ig a)B (66) P at p -t 2 at The fourth term is integrated with the assumption of no variation in temperature along the vertical thickness of the aquifer (i.e., T T TF). Hence 3T aT nn B (67) n T at T= nTB at Finally we assume that VT q rVT + qr T q; qr VT (68)

31 i.e., neglecting the average of the product of deviations, qr VT, of qr and VT from their respective averages. r r Then, with BVT = J (VT)dz = V' (BT) + T F VF2 B'T 169) 2 1 T 1 b we obtain q.VT Bq'.V'T (70) r r Accordingly, the integrated form of (61) is WVB at +B i+ B +B DB'Bq + V'B + B + B Bfi - r at at +fp at p z t - 3 T ~Bq~'IT = 0 (71) nfTB at - TB r If we assume vertical consolidation only, i.e., U' = 0, with 1 3B 1 3B a and a = p B 3- T = ap equation (71) reduces to + [lq 1[ 23 + c T] + fB B — V- Bqr + [1 + nSpp ] [ap at T at p at dT _- nTB T mBq'.VT = 0 (72) T at In general -SpB/2 << 1. By neglecting in (72) the term fip B/2, we obtain an equation which could also be obtained by integrating (13).

32 In (72), we shall use (65) to obtain the approximate expression k B Bq' - B' - B - (V + gV); z = b + (73) i p 1 2 P~ f where z = (b1 + b2)/2 and a hydrostatic pressure distribution has been assumed along B, with a fixed bottom at b1 = bl(x,y). The Integrated Energy Conservation Equation We shall integrate (31). Assuming T T - F T F1 I we obtain: b2 T aT aT (PC)m t dz = B(pC)m t - B(PC)m T (74) b With (49), we obtain for the second term b2 b2 BV(A ) V'-J (4.V'T)dz f bV (Am VT)dz = B-v(mvT) =b m T) IF2 VF2 m ( VT) F.VF1 = Am -V'IBA *VT (75) The integral of the third term yields the following approximation (2 ( -VV) dz B V'T (76)

33 where, as before we assume T = T IF T I and we neglect 1 2 averages of products of two or more deviations (^) from their respective averages (i.e., we neglect faCqrVT, pfCvqr VT, pfqrCv VT, etc). For the fourth term we obtain the approximation 2T T Jb {(T - p)Vq r} dz (T() - r)'B (77) 1 P The last term on the L.H.S. of (31) is estimated by 2 {n(T n)aT b fn(T T- P) + (1 - n)Ty) - dz 1T {( ) - ) + (1 - i)T} (V + t (78) where 3a/at is expressed by (63). Hence, the integrated form of (31) is ~T T B(PC) _ V*,B&_ V'T +p C BqI'V'T + ( -B T, ~ mam fVr + {n(T() - p) + (1 - n)Ty}(V'-B t + a) = (79) p Equation (32) can also be integrated in a similar way. The averages Cv, m etc. can take into account variations of these parameters along the vertical, say in a layered aquifer.

34 The Integrated Equilibrium Equations The point equilibrium equations are (34b). Making use of (57) and (58) and following the detailed presentation of the authors in an earlier paper (Bear and Corapcioglu, 198Qc), we now integrate (34b) over the vertical thickness of the aquifer. In view of the boundary conditions expressed by (56), we obtain BI(e + BID'e - I Bpe 0 (80) xx xx'ay =x / Bale + _a Bp 0 (81) ax yx ay yy ay Br e + Ba'e 0 (82) -x xz ay yz in which all averaged values are functions of x,y and t only. Next, the thermoelastic stress-strain rleations given by (36) are expressed in terms of average excess stress components and average displacements. Assuming, that U' U'F = U' I i.e., practically the same horizontal displacement along the vertical, we obtain au au A~=E +~ + + + =83y xx yy zz 3x y B(83) where B b2 b = (bz 0 + UF ) -(b (b + U-Fb) ( ) + (UZIF - UZIF ) = B~ + A, B~ is the initial value of B; 2 1 A is the change in B due to the movement of both F2 and F1. If bl is fixed, -Az is the land subsidence.

35 r-_ (au au xx xx + F~- T= (X + 2G) X + A ( Y + (84) ae = 2GE + X - Te = x + (( + 2G) x — Z Te (85) au au xy ax a + [Z a aG 2 U | 1 a'e-G( (+ 2 G ( 8) y_ yz y yz)B zy - ---- By UzIF 1 D,ye e= G(._ +- ay4 + x (R + (88),e+ z -e zx G z -,, 2, yT =X + + (X + G) B yT ZZ Dz ax ay B (89) By inserting these expressions into (80) and (81), and linearizing the resulting equations, assuming constant X, G, and y(in fact A, G, and y) and B(x,y,t) = B (x,y) + Az (x,y,t), AZ << B~, we obtain Gz a-G ax -o a/P-)e =~9 GV' U + (U + G) G (z/Y B = O (90) ~' ay ay ax+y or au au a3O (A /B ) e aT GV2' U + + G) + Y + z (92) (+ a G+(A /BO) - e _ e a y 3y ayz z! 2-~~~ ~~~ B [ \

36 Equation (82) becomes D z aB 2 -U DaX BG aTX + G(Uz ax) + ZIF x2 F1 ax {BG + UF U y {BG -y + G(Uz y + UzIF2 ay ZIF1 (94) Studying (94), we note that it contains information which is usually unavailable, (i.e., UZIF and UZF ). Therefore, 2 1 we have to introduce additional (albeit, simplifying) assumptions. One assumption could be that the bottom boundary is stationary (i.e., b1 = b (x,y)). Then UZIF = z = -d, where 2 6 is the land subsicence which is considered positive downward. In this way, the number of variables in these equations is reduced to four: 6, U, U pe Another approach (Bear and Corapcioglu, 1980, following Verruijt, 1969) is to introduce the assumptions that the stress distribution satisfies the condition of plane incremental total stress. This means e e e e e Czz =, xz = Czx = 0, cyz = czy = 0 (95) zz xz zx yz zy Bear and Corapcioglu (1980c)discuss this condition and its justification in detail. The variables in (95) may be replaced by averaged ones. The first of (95), combined with the definition of effective stress, leads to Me -= p (96) zz

37 and hence, (89) yields ~e x z _ p = ( x++ i) + 2G) e (97) Equation (82) now vanishes due to the second and third equations of (45), but with (97) replaces (82). We have again four variables (Ux U, A and pe) and four equations: (71), (92), (93) and (97). If we assume no average lateral displacements, i.e., U' = 0, (97) reduces to A -e z -e p (X + 2G) - yT (98) COMPLETE MATHEMATICAL MODEL The complete mathematical model consists of the following equations: * The fluid mass conservation equation (71), with q' expressed by (73) The energy balance equation (79) * The equilibrium equations (92) and (93) * Equation (97), selecting p to U' A and se These equations have to be solved for the variables * e ~e ~ ~jy -e -y * p, T, Ux, Uy, z However, we note that these equations contain also p and T as dependent variables as well as pressure and temperature dependent fluid and solid properties, e.g., p = p(p#,), u = j(T), n = n(p), k = k(n), and B, so that we have to add the p, T, and B as additional variables to be determined from

38 B = B~ + A (99) T = T~ + Te (100) P = PO + P (101) where po and T~ are the initial steady state distributions of pressure and temperature. They may be obtained from field observations or by solving V'Bq'O - TB-qo VTo = 0 (102) Bq'0 = - B (V'p~ + p~gVz~) (103) -V' ~BOI' 0 V'T~ + p[C BoqV'T + [T () - p~]V-Bq~ = 0 (104) f v r In order to follow the continuous change in i, as consolidation occurs, we integrate (4), making the usual assumptions with respect to the bounding surfaces F1 and F2, and n!F nlF. We obtain 1 2 B t = (1 - n) {a_ (B )105) =i) {-y- (B t)+ (B +j(105) where we have assumed 3n/3t >>Vs Vn. Using this additional equation, we may consider ni as an additional variable. By using (99) and (105), say in a numerical solution, ni(x,y,t) and B(x,y,t) can be determined at every time step. We observe that the equations comprising the model are non-linear, and therefore simple superposition of incremental changes and initial conditions is not possible. This nonlinearity may be even more significant if the considered problem involved

39 a large range of temperature and pressure variations. For example, the effect of temperature on the dynamic viscosity of the fluid (appearing in (73)) may be very significant, say in in hot water injection operations. The number of variables combined with the model's nonlinearity preclude any analytical solution. Numerical methods have to be employed. We may attempt to simplify the problem by introducing certain linearizations, say with respect to B. For example, since (A z/at)/Az >> (aB/aT)/B, we may assume aB/aT - Ba(Az/B)/3t. Also we may use the approximations zU' 3U' V-BqB ~BV-q'; V-' B -- B Equation (71) then reduces to 3U'A Vi r' + V't + (1 + nSp) at - nT at V r 3t _B__Z -+ nqp - = 0 (106) p 9t T r In a similar way we may linearize (79) with respect to B. We obtain Im ~f / (PC)m - m V'm T + f qrV'T [T( T - q p T ~ x Az + {i[T( ) p] + (1 - n)T yt ( + + -) =0 pat ~x a~3y B0 (107) In an earlier paper (Bear and Corapcioglu, 1980b) the 3U DU A authors show that ~= X + Y+ ax 9y B0

40 3- auv A e 2- x aa+)_ z _ GV'2 U + (X + G) [ +ax y ax Y (108) ~ 2 z e- 3T GV' U + (X + G) [ + yYi] + aT 0 (109) y ax ay ay BO.Y = The plane stress assumption leads to au au A -e x z -~e (-x + ay + (X + 2G)BZ - yTe (110) A similar linerization of (105) leads to = (1 - (axX + = (1 - ) (111) (I - 7t) ff = + n ) -e -e Altogether we have six equations to solve for p, T U, U, A and ni. x V z To complete the statement of the mathematical model, we require information on initial and boundary conditions (in the xy plane!) of pe, Te Us and A The above model has been developed for a confined aquifer where the confining layers are assumed also to be insulators with respect to heat flow. There is no difficulty in replacing, say the upper boundary by a boundary through which leakage of water and loss of heat may take place. This means that in the mass balance equation (71) we have to add a source of water representing leakage say, from an overlying aquifer (expressed in terms of the pressure difference across the semipervious layer). Similarly, we have to add in the energy balance equations,

41 terms representing sources of energy due to convection (by the leakage) and conduction through the semi confining layer. Actually, we do not have to introduce them as an after thought. They will appear automatically in the integration process. TEMPERATURE AND PRESSURE DEPENDENT PARAMETERS To solve the model for a particular case of interest, information is required on the parameters a1, aT Bpp BIT' ~i pf, k,', Cm', A, and G which appear in it. Some of them are sensitive to temperature and pressure changes. For example PfI A' Cv',T change considerably with temperature. Some like p vary also with pressure too. These effects further complicate the solution of the model. The equation of state for the fluid (say water) can be expressed as a function of temperature and pressure in various forms. For example, [Fernandez, 1972] pf = p exp[-T(T - T~) + Sp(P - P0)] (112) where pO is the fluid's density at some reference state (TO, PO). By retaining the first order terms of the series expression for (112), we obtain Pf = P~[l - T (T - T~) + S (P - PO)] (113) Sorey [1978] suggests that the pressure dependence of density in liquid (single phase) geothermal reservoirs can be neglected. An expression given by Wooding [1957] could then be used

42 pf p [1 - T) (T - T~) - (114) In equations (112) through (114), T is in degrees Celcius. The liquid's viscosity is also strongly dependent on temperature. Bird et al. [1-960, p. 327] suggests D = exp[A(- )] (115) T T~ in which A is a constant. Huyakorn and Pinder [1978] used an expression to compute the viscosity of water. s= 10-6 x 10[248.37/T + 133.15] (116) where i is given in gr/(cm.sec) and T is in degrees Celcius. (Recall that everywhere else in this paper, T is the absolute temperature.) The fluid's thermal expansion coefficient ST can be estimated by an expression given by Sorey [1978] =0 (llV) =T p= (T - TO) (117) where f is given in (gr/cm3) and T is in degrees Celcius. Heat capacity of the fluid at constant volume C-, is also temperature dependent, and is given in tabular form by Dorsey [1968]. Other parameters of the model can be approximated as constants to a certain degree of accuracy. EXAMPLE: A MODEL FOR A SINGLE PUMPING/INJECTION WELL When we consider thermoconsolidation in the vicinity of a single pumping/injection well (rate Qw' radius r ) extracting W

43 hot water from an infinite (liquid) geothermal reservoir, or injecting hot water into an infinite cold water aquifer for storage, or cold water flooding of a warm reservoir, and we assume a situation of complete symmetry with respect to the well, we write all linearized averaged flow, energy and equilibrium equations in radial symmetry. For this case, we obtain from (65) k Bp p k ap B k p _ a Br~rr r [r + ] (118) where z = (b1 + b2)/2. Assuming here and in (126) below that ~ >>,p - - (and also for Pe), (118) reduces to (qr)r ~-~-2(119) rrr B~ In (106), < Pf <<1. Hence this equation becomes k ~p. ai~p ~ T - k D_~ aT -[r + + n t- r Dr Dr Dt P -t T at T Dr Dr (120) where we have introduced the dilatation defined by A DU Ur A 1 (U) + r + z r Dr (r B- = r r B~ (121) The equilibrium equation is DU ar U _e ~ e r + ) 1 r r _ T _ (2G +r - 2G( + r (122) r Dr r Dr Dr where e e =T T-T0

44 The plane stress assumption leads to A au U A = (2G + X) + X( r T = 2G Z + E - e (123) Bear and Corapcioglu (1980b) show that from (122) we obtain A z =e +, (2G + X)E - 2G 0 = p + T + 2g(t) (124) where g(t) is a function of time to -be determined by the boundary conditions. A comparison of (124) with the total plane stress assumption (123) gives A g(t) = 2G[V' U' - z (125) In radial coordinates:DU U A A r r z (2Z g(t) = 2G[ + 2G[ - 2 BO(126) Equation (107), can be rewritten in the form a" _ 1 9ii (127) The energy equation reduces to aT 1 a aT k a9 a (PC) - [r - C m at r 3r -m1 r f V_ ar ar - 1 a k_ - [T() - (r r +T - + - 0 (128) + i[( + ( - )T (128) p~~~~~~~~~~~~~~18

45 Equations (120) through (123), and (128) constitute now a set of equations to be solved for p, T, E,U and A. In a r Z single well problem the solution is subject to following boundary and initial conditions t < 0, r > r = PO T = TO Ur, A, = 0 (129) t > 0, r e Qw (130) ~rwe P_ =(130) 3r -2frr Bor = rw Ur =0 (131) Boundary conditions for temperature at the well can be obtained from the continuity of heat flux at the well. 27r Bw~(qr)r C T = 2r BO(qr)r C T - 2rrr BO~ AT w r r f v w w r rpf v w m~r (132) The terms on R.H.S. of (131) are convective and conductive fluxes, respectively. Tw is the temperature at the well. In the case of hot water injection into the aquifer, the condition is QwP C e ~e wfv _T t > 0, r = rw (Tw T ) Br (133) w w2rrB r(13 In the case of hot water withdrawal, T Tw, hence t >0, r = rw (134) t > 0, r = r aT w 0 (134) ar The amount of water injected or pumped out is w = 2r B0~(q ). w r r

46 At infinity, the effect of injection or withdrawal vanishes. Hence r- -a -e -e e ~r -, p IT ~ A, U = 0 (135)

47 SUMMARY AND CONCLUSIONS By averaging three-dimensional coupled equations of mass, energy and equilibrium over the vertical thickness of an aquifer, and introducing thermoelastic stress-strain relations, we obtain a mathematical model composed of a set of partial differential equations in terms of pore water pressure, temperature, vertical subsidence and horizontal displacements, as dependent variables. The development is based on Terzaghi's concept of effective stress, and assumptions of essentially horizontal flow in the aquifer and shear free boundaries. The effects of viscous dissipation and compressible work have also been included. The model developed in this study can be employed to simulate hot water injection into an aquifer for energy storage, as well as for geothermal production by pumping hot water from a reservoir or for cold water flooding a warm reservoir. In the first two cases the mobility ratio (which for a single fluid saturated reservoir is defined as the ratio of viscosities of the displaced fluid to the displacing one) is greater than unity. This produces instabilityl (fingering) at the advancing front. This phenomenon is neglected in this paper. In the last case, the mobility ratio is less than one, and no fingering occurs. So far we have been discussing the case of a single liquid saturating the porous medium. In principle, the methodology presented here is applicable also to multiphase flow (e.g.

48 oil and water in an oil reservoir, or water and steam in a geothermal reservoir), however, certain modifications have to be introduced. The mass conservation will be treated separately for the two fluids. This means one (macroscopic) mass conservation equation will be developed for each fluid. Each fluidwill have its own pressure and its own saturation. The saturations will constitute additional variables (with their sum = 1) and the difference between the pressures of the non-wetting and wetting fluids (two additional variables) will be given by the capillary pressure which is a function of the saturation. The effect of phase changes may also be added. The relative permeability of each of the fluids depends on its saturation. For the energy balance equations, we shall usually assume that the temperature (at the macroscopic level) is the same for both fluids and for the solid, and combine the individual equations into a single equation for the porous medium as a whole with the two fluids regarded as a mixture. The coefficients of heat capacity and conductivity will be functions of the varying saturations. To obtain the equilibrium equations, we shall regard again the two fluids as a mixture. The effective stress in this case will be given by a = a' -X(S )PwI, where.X is some function of the -~ w wz water's saturation [Narisimhan and Witherspoon, 1977; Bear and Corapcioglu, 1980c]. In order to perform the integration of the three-dimensional equations, we need boundary conditions on the top and the bottom

49 boundary surfaces with respect to the flux of each of the fluids and the thermal flux and stresses for the system as a whole. Some of these modifications are presented by Faust and Mercer (1979). Due to nonlinearity of the equations (even after they have been linearized with respect to B) and because of the solution dependent n,p,k, p, etc., superposition is not applic — able. A further attempt to linearize the equations, with parameters taken at some initial or intermediate pressure and temperature values would fail due to product terms. When numerical solution techniques are employed, the continuous changes in these parameters have to be taken into account. In the derivation of the mass conservation equation we have started from the known three-dimensional equation, in which the mass dispersive (<peVf > - DVPf) and diffusive (-DdVPf) fluxes have been neglected. In the next step of averaging over the vertical, we have neglected macrodispersion (=pfq) (Bear, 1979, p. 256) due to both specific discharge and density variations along the vertical. Actually, in order to take macrodispersion into account we should have integrated V'pfq + anpf/Dt=O over the vertical, leading to pfq = pfq + fq, where the last term is the macrodispersive flux. These mass fluxes can easily be added to the equation, if so desired. It is assumed here that they are much smaller than the convective flux. The macroscopic level is also the starting point for developing the equilibrium equations.

50 For the energy balance equation, we have started from the microscopic level and averaged over the vertical, in order to emphasize the various assumptions involved in employing the balance equation at the macroscopic level. In the macroscopic equation, the thermal dispersive flux <(pfC*OT*)(V*)>f has been neglected, assuming that it is much smaller than the convective flux expressed by <p*CvT* > < V > f Similar considertions apply also to the solid matrix. At the next step of averaging, we have neglected the thermal macrodispersive flux (PfCVT)q. Again, the assumptions are that these are relatively small terms. In employing the stress strain relationship (39), we have assumed that although they are strictly valid for a solid continuum, they are also valid, at least as far as their structure is concerned, for the solid porous matrix. In doing so, the two Lame coefficients G and X are no more those of the solid continuum, but new coefficients of the porous matrix to be determined experimentally. Another alternative is to obtain the macroscopic constitutive equations by averaging the microscopic onesover the REV. The same observation is true also with respect to some of the other coefficients. To simplify the discussion on the equilibrium equations, we have assumed that the thermoelastic solid is isotropic (although we maintained the symbols indicating the tensorial nature of k and Am If the thermoelastic body is anisotropic

51 (and it is usually so if its structure is such that permeability is anisotropic, and hence also thermal conductivity), we have to use (35) instead of (36). Instead of two Lame constants, Cijk depends on more basic constants depending on the kind ijkk of symmetry we assume to exist (see any text on the theory of elasticity). Finally we would like to emphasize again that in view of the various assumptions and approximations, the results, which may be obtained by solving the model by numerical techniques, should be viewed as estimates. These estimates, however, should be useful for assessing the effects of various pumping and injection schemes on vertical subsidence and horizontal displacement. Obviously, much more work is still needed in order to determine the values of the various coefficients and parameters- appearing in the equations.

52 REFERENCES Aktan, T., and S. M. Farouq Ali, Finite element analysis of temperature and thermal stresses induced by hot -water injection, Soc. Pet. Eng. J., 457-469, 1978. Bear, J., Dynamics of fluids in porous media, American Elsevier Co., New York, 1972. Bear, J., On the aquifer's integrated balance equations, Adv. Water Resour.,,(1),15-23, 1977. Bear, J., Hydraulics of Groundwater, McGraw Hill, New York, 1979. Bear, J., Technical comment "On the aquifer's integrated balance equation," accepted for publication in Adv. Water Resour., 1980. Bear, J., and M. Y. Corapcioglu, Mathematical model for regional land subsidence due to pumping, Part 1: Integrated aquifer subsidence equations based on vertical displacement only, submitted for publication, 1980a. Bear, J. and M. Y. Corapcioglu, Mathematical model for regional land subsidence due to pumping, Part 2: Integrated aquifer subsidence equations for vertical and horizontal displacements, for publication, 1980b. Bear, J., and M. Y. corapcioglu, Centrifugal filtration in deformable porous media, submitted for publication, 1980c. Bear, J. and G. F. Pinder, Porous medium deformation in multiphase flow, J. Engrg. Mech. Div. Amer. Soc. Civil Eng., 104 (EM4), 881-894, 1978. Bird, R. B., Stewart, W. E., and E. N., Lightfoot, E. N., Transport phenomena, Wiley, New York, 1960. Brownwell, D. H., Garg. S. K. and J. W. Pritchett, Governing equations for geothermal reservoirS Water Resour. Res., 13, 929-934, 1977. Corapcioglu, M. Y., and N. Karahanoglu, Simulation of geothermal production, Proceedings of Second Miami International Conference on Alternative Energy Sources, Miami, Hemisphere Publ. Corp., Wash., D. C., 1979.

53 Derski, W., and S. J. Kowalski, Equations of linear thermoconsolidation, Arch. Mech., 31, 303-316, 1979. Dorsey, N. E., Properties of ordinary water - substance, Hafner Pub. Co., New York, 1968. Ertekin, T., Numerical simulation of the compaction-subsidence phenomena in a reservoir for two-phase non-isothermal flow, Ph.D. Thesis, The Pennsylvania State Univ., Penn, 1978. Faust, C. R., and J. W. Mercer, Geothermal reservoir simulation, 1. Mathematical models for liquid and vapordominated hydrothermal systems, Water Resour. Res., 15, 23-30, 1979. Fernandez, R. T., Natural convection from cylinders buried in porous media, Ph.D. Thesis, Univ. of Calif., Berkeley, Calif., 1972. Garg, S. K., J. W. Pritchett, M. H. Rice, and T. D. Riney, U. S. Gulf Coast geopressured geothermal reservoir simulation, Final Report, Systems, Science and Software, Rep. SSS-R-77-3147, La Jolla, Calif, 1977. Hassanizadeh, M., and W. G. Gray, General conservation equations for multi-phase systems, I. Averaging procedure, Adv. Water Resour., 2, 131-144, 1979. Horne, R. N., and M. J. O'Sullivan, Numerical modelling of a desaturating geothermal reservoir, Numerical Heat Transfer, 1, 203-216, 1978. Huyakorn, P., and G. F. Pinder, A pressure-enthalpy finite element model for simulating hydrothermal reservoirs, Math. Comp. in Simulation, 20, 167-178, 1978. Lippmann, M. J., T. N. Narasimhan, and P. A. Witherspoon, Numerical simulation of reservoir compaction in. liquid dominated geothermal systems, Proceedings of Second International Symposium of Land Subsidence, International Association of Hydrological Sciences, Anaheim, Calif., Pub. no. 121, 179-189, 1976. Mercer, J. W. and C. R. Faust, A review of numerical simulation of hydrothermal systems, Hydrol. Sci. Bull., 24, 335-343, 1979. Narasimhan, T., and P. A. Witherspoon, Numerical model for saturated-unsaturated flow in deformable porous media, Water Resour. Res., 13, 657-664, 1977. Nowacki, W. Dynamic problems of thermoelasticity, Noordhoff Int. Publ., Leyden, The Netherlands, 1975.

54 Pritchett, J. W., S. K. Garg, D. H. Brownwell, L. F. Rice, M. H. Rice, T. D. Riney, and R. R. Hendrickson, Geohydrological environmental effects of geothermal power production, Phase IIa, Systems, Science and Software, Rep. SSS-R-77-2998, La Jolla, Calif., 1976. Sorey, M. L., Numerical modelling of liquid geothermal systems, U. S. G. S. Prof. Paper 1044-in, 1978. Verruijt, A., Elastic storage of aquifers, in Flow Through Porous Media, edited by R. J. M. DeWiest, pp. 331376, Academic Press, New York, 1969. Whitaker, S., Diffusion and dispersion in porous media, J. Amer. Inst. Chem. Eng., 13, 420-427, 1967. Wooding, R. A., Steady state free thermal convection of liquid in a saturated permeable medium, J. Fluid. Mech., 2, 273-285, 1957.

55 LIST OF NOTATIONS A = phase average = U (o) dU 0 0 < A >= intrinsic phase average = U oU )A d blb=2 elevation of bottom and top aquifer bounding surfaces, respectively. B = thickness of aquifer ( = F 2 = b2 - b1 C = fluid's heat capacity (specific heat) [per unit volume] at constant volume C, = solid's heat capacity (specific heat) at constant strain CijkZ = (i k) IT = constant d( ) dt = total derivative of ( ) with respect to the moving water =a ( ) - -- *v V- ( )]'t' w ds( dt = total derivative of ( ) with respect to the moving solid = ( ) +Vs ~*( )] e = (as subscript) denotes excess value f = (as superscript) denotes fluid f = body force acting on saturated porous medium F =O, F2=0 = equations describing top and bottom surfaces, respectively g = gravitational acceleration G = a Lame constant hf = fluid's enthalpy J = conductive flux k = medium's permeability n = porosity o = (as superscript) denotes initial steady state value

56 o = (above a letter) deviation from intrinsic phase average p = pore water pressure q = specific discharge 4q = specific discharge relative to the solid QW = rate of withdrawal or injection r = radial coordinate s = (as subscript) denotes solid t = time T = temperature T = tortuosity u = velocity of microscopic fluid solid interface, velocity of surface bounding an aquifer from above or below U = solid displacement V,V = water and solid velocities, respectively -f -s w = (as subscript) denotes water aa = (a(/aP) IT = constant "T p == (ae/aT) p = constant B = fluid's coefficient of compressibility (at constant temperature) ST = fluid's coefficient of thermal (volumetric) expansion (at constant pressure) Bij = - (aa'e/T) ~ = constant Yij ~1j s =i constant =( - 3aT ) E constant sT

57 A = B - B~ z = Kronocker delta ij ~.. ~ = solid's strain 13 ~ = volumetric strain (dilatation) X = a Lame constant XfXs = heat conductivities of fluid and solid, respectively. XfaX = heat conductivity of fluid filled void space and solid matrix, respectively Am -= thermal conductivity of saturated porous medium zm VP = dynamic viscosity of water p = density (IC)m = heat capacity (per unit volume) of saturated m porous medium a = total stress oa' = effective stress T. = stress in fluid k* = (as superscript) denotes a microscopic value = (over a letter) average over the vertical = (over a letter) deviation from average over the vertical =over a vector, a tensor, or an operator (V',V' ) denotes vector components or operators in a xy plane only

58 LIST OF FIGURES Figure 1 - Nomenclature of Volume Figure 2 - Nomenclature for a single-layered aquifer

REV U0, otTid A~ f t uid

60'if, ^F2= 0 F =0 x lyr ur Y D tur d~~ e. ~1, ~ e Ii ~ t ~ ~ t ]11 ~' -"'' ~r ~ C,r ~~, ~ ~V bl(Xb x~,y t ~~~~~~~~~~~~UNVRSIT`. Y OF ICHIGA 111111111IIII1ll11l1lllllllll1llllllllll1l1 II1IIIII1I1IIIIIIII 3 9015 02224 7350