Numerical models in land subsidence control

COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 5 (1975) 227-237. 0 NORTH-HOLLAND PUBLISHING COMPANY

NUMERICAL MODELS IN LAND SUBSIDENCE CONTROL Giuseppe GAMBOLATI IBM Scientific Center, S. Polo 1364, Venice, Italy

Received 22 April 1974 Revised manuscript received 1 October 1974

One consequence of growing concern over subsurface fluid pumping is the irreversible soil compaction and the related land subsidence. The mechanism of subsidence is an interaction between the stress field and the flow field operating within a deforming and permeable system. The use of high-speed digital computers appears quite promising in its ability to handle large geological systems and to predict the impact on the environment of an intensive exploitation of the underground resources. The models described in this paper can predict land sinking due to fluid withdrawals in natural two or three-dimensional basins and reservoirs. They can be applied at both the reconnaissance and intermediate stage to assess the convenience of alternative pumping schedules in relation to the permanent modifications that land subsidence produces in the environment.

1. Introduction

Regional land subsidence induced by human activities is developing in many areas, and the related consequences are destined to become a matter of growing concern in the next few decades. Among the major causes of land sinking are the extensive withdrawals of fluid from subsurface resources. Man has been using natural underground wealth for a long time but it is only recently that technological progress is allowing for large exploitations which are no longer counterbalanced under natural conditions. This results in a modification of the environment which is often irreversible or reversible at the expense of enormous costs and great waste of human energy. Land sinking is one of these consequences. It is especially dangerous for the areas lying at sea level because of the risk that stretches of coast are permanently submerged. Major subsidence of several oil-gas fields has been reported in the literature [ 1 ] . We can recall Wilmington oil field in the harbor area of Los Angeles and Long Beach, California (more than 8 m from 1937 to 1962), the oil fields on the shore of Lake Maracaibo in Venezuela (3.4 m from 1926 to 1954), Goose Creek oil field in Texas (almost 1 m from 1918 to 1925). Gas fields at Niigata, Japan, and in the PO Delta, Italy, experienced rapid subsidence in the late fifties. Lately subsidence is expected to develop in the Groningen gas field, the Netherlands [ 21. Significant surface lowering due to groundwater pumping has been observed in San Joaquin Valley, California (8 m from 1935 to 1966), in Mexico City (8 m from 1936 to 1968), in Santa Clara Valley, California (4 m from 1920 to 1967), in Tokio and Osaka, Japan (3-4 m from 1928 to 1965). A very recent case of particular interest is. Venice [ 3 I. There the quantities of water extracted have been less than those that caused considerable sinking in the cited areas. The lowering, though

228

G. Gumbolati, Numerical models in kznd subsidence control

modest, is nevertheless a matter of concern because it compounds the effect of the exceptionally high tides in the city during winter. Many other gas-oil fields or hydrologic basins above sea level and not subject to inundation probably have shown subsidence, but with no marked deleterious consequence, therefore they have not been publicized. The mechanism of subsurface compaction leading to the surface lowering is an interaction between the flow field and stress field within the deforming system. The artificial withdrawals of fluid modify the balanced natural conditions by causing a decrease of the pore pressure and therefore an increase of the effective intergranular stress. Under the influence of the increase of the effective stress the formations compact and the land surface subsides. Regardless of the nature of the fluid removed (gas, oil or water), the principles underlying the physical occurence are the same. However the geological setting can be different according to whether a gas-oil reservoir or a hydrologic basin is being analyzed. This accounts for the different models which are applied in either case. If representative values of the soil and formation properties are available, we will be able, to analyze these complex natural systems. One of the most powerful and effective methods of analysis involves the use of the mathematical and numerical prediction models which are formulated from the equation of subsurface flow combined with the equations of the elastic equilibrium and then solved nume~cally with the aid of a digital computer. With such models it is possible to determine the changes in the flow field due to well pumping and the compaction of the various formations in relation to both depth and time. The models can be applied to predict the effects owing to different pumping schedules and to assess the suitability of alternative actions and the consequences of various proposed designs. Models can also provide guidance in the building of monitoring systems and even in the preparation of appropriate legislation and regularoty codes.

2. Numerical models to control land subsidence It has long been understood that land subsidence is best analyzed with reference to the theory of consolidation, wherein it is recognized that consolidation represents the response of a compressible porous medium to changes in the fluid flow field operating within it. Till recently almost all mathematical models of Land subsidence have been direct applications of the Terzaghi consolidation theory [ 11, [4] - [ 61. This simplified theory has been adapted also to the computation of oil field compaction (e.g. [ 71). This approach, in the most sophisticated applications, involves analytical time-dependent solutions to the boundary value problem based on the one-dimensional vertical equation of flow in which the piezometric declines are specified as boundary conditions. An excellent review of this approach as it applies to land subsidence is provided in [ 81. In a more recent review, presented at the Tokyo symposium on land subsidence, Poland [ 91 refers to results based only on the classic on~dimensio~al theory, and he notes the need for some sophisticated mathematical or numerical models. In the entire symposium proceedings only one paper [ 101 presents an analysis based on a model of greater complexity. A complete analysis of land subsidence requires the determination of the three-dimensional deformation field that accompanies the three-dimensional flow field. Schiffman et al. [ 1 I] have

G. Gambolati, Numerical models in land subsidence control

229

classified condolidation theories into: 1) one-dimensional Terzaghi theory, 2) pseudo three-dimensional consolidation theory based on the independent solution of the diffusion equation, 3) complete three-dimensional consolidation theory based on the Biot equations. In the remainder of this section we will develop the second and third approach as they apply to land subsidence and are used by the digital computer. starting with the latter and treating the second as a subset of the third. We shall show that the third approach must often be dismissed on the grounds of impracticality, whereas the second offers at present the best trade-off between data availability and computer limitation. We shall give indications as to how the working mathematical relations are modified according to whether a gas-oil reservoir or a hydrologic system is being simulated. Finally in the next section we shall describe two representative applications recently made at Venice and in the Netherlands to predict land subsidence that was caused respectively by groundwater pumping and gas production. 2. I. General three-dimensional consolidation equations The consolidation equations were first derived by Biot [ 121, [ 131 who stated the interrelations existing in first approximation between the three-dimensional flow field of the fluid and the threedimensional stress field of the medium. Neglecting the compressibility of the single grains (usually much smaller than that of the rock bulk), Biot’s theory relates the components of the rock matrix deformations u,. to the incremental “effective stresses” ulk acting upon the solid skeleton and the incremental pore pressure p. Including the equation of flow, the results take the form of the following set of equations:

(1)

where ‘ik

= $

‘ik,jh

‘jh

.

In eq. (2) Cik,jh are the components of the elasticity tensor (recall that only 21 coefficients are distinct for real anisotropic media). The principal components of eii are linked to the displacements Ui bY aui ‘ii = G

I

’

(3

If we replace the ujk in eq. (1) by their values taken from eq. (2), keeping in mind eq. (3) for the strain components, we obtain a set of four equations in four unknowns: the incremental fluid pressure p and the three components of the displacement vector z+_Solutions of the system will provide the distribution of ui (in particular zd3at the.ground surface is the land subsidence). At the theoretical level the Biot equations offer the most sophisticated approach to the prediction of land subsidence. Their Achilles heel is the large number of formation parameters needed: the 2 1 classical coefficients of anisotropic elasticity Crk,jk, the 6 coefficients of anisotropic hydraulic conductivity kii and porosity n - 28 parameters in all. These parameters are needed for each formation in the case of non-homogeneous media. A second important drawback is the prohibitive amount of storage required to solve Biot equations numerically in natural three-dimensional compacting systems. Even if homogeneity and isotropy for both hydraulic and mechanical properties are assumed throughout the medium, as wan done for instance in [ 141, the limitation still remains. The finite element solution supplied in [ 1O] and I 15 1 does not go far beyond the analytical fo~ulation. The relevant applications are necessarily reduced to small and unr~listi~ally simple configurations, The practical applicability of the Biot system of equations to regional land subsidence problems is at present unfeasible. A worthwhile alternative is provided by the approaches described below.

We can partially remove ourselves from the horns of the previous dilemma by searching for the solution of the eqs. (1) through two related steps. In the first step the equation of flow is solved independently of the remaining equations, thus providing the changes in the flow field. In the second step the outcome from the flow equation is used to compute the displacement distribution and ultimately the surface sinking. To do so, the necessary condition rests on a major a~umption~ namely the volume strain e at a point of the system depends only upon the pore pressure variation p at that point and not upon the pressure changes occuring in the surrounding medium. The worth of this assumption has not yet been tested, but it underlies all developments of the classical flow equation (either in ~oundwater or oil-gas field). The appropriate formula is

where a is related to the elastic properties of the medium. Its value is _(l

--2rr)fl +v) E(1 -V) ’

where E and v are the moduli of Young and Poisson respectively. Making use of eq. (4), the last relationship of the system (1) turns into the well-known parabohc differential equation whose solution is now more easily attainable than that of the complete Biot set of equations. The finite difference scheme is the technique more fr~q~~ently used, but recently the new powerful finite elemenj method has been developing to an increasing degree even in the fluid flow field. The results, expressed in terms of time-dependent pressure distribution within the system, are then used as prescribed disturbing stresses in a three-dimensional calculation of the deformations. To this end

G. Gambolati. Numerical models in land subsidence control

231

the classical equation of elasticity (the first three of eq. (1)) are to be solved. This may be done numerically, but in the case of homogeneous media one can resort to the concept of “strain nuclei”. With such an approach the subsidence is evaluated by integrating the Green function Us?for the displacement field over the entire volume in which changes of p occurred. Geertsma [ 161 provided the vertical displacement U; due to a unit pressure change in a unit volume (tension center) of a semi-infinite homogeneous medium. Gambolati [ 171 extended the formulation so as to allow for a heterogeneous tension center in the homogeneous semi-infinite medium. The value of U; proves to be h - 2 + z + 3h - 4v(z + h) + 6z(z + h)2 P

*3

P

= AFch

*5

z r v)

79,

,

(6)

where the meaning of the symbols is explained in fig. 1 and A is a coefficient depending upon the elastic properties of the nucleus and medium A=

(7)

Y

-

tension

center

cz

Fig. 1. Scheme of the semi-infinite medium with the tension center.

In eq. (7) the subscripts A4and N refer to medium and nucleus, respectively, and cb = 3( 1 - 2v)/E is the compressibility of the rock bulk. Integration of eq. (6) yields the vertical displacement

V

Relation (8) may advantageously be applied to provide the subsidence caused by oil production. In fact the relatively small dimension of the reservoir, its depth and geological setting (usually top and bottom impermeable, consolidated sediments etc.) produce a three-dimensional displacement field, and the effects at the land surface are substantially different from those observable in a hydrologic system where the subsidence is essentially equal to the basin compaction. The main limitation of eq. (8) lies in the requirement of homogeneity for the semi-infinite medium. This approach has been recently used to predict the sinking of the Groningen area, the Netherlands, where a gas field is being developed [ 183. Integral (8) is evaluated numerically over an arbi-

232

G. Gambofati, Mumericai models in land subsidence control

trary-shaped reservoir. It is replaced by a finite sum of the discrete contributions from the small elements AI’ into which the reservoir is subdivided. Thus we obtain at the ground surface

where X, , yn and Z, are the coordinates of the nth reservoir element in which the pore pressure is decreased by p,. It should be recalfed that integral (8) is mathematic~ly correct for any shape of V when the reservoir has the same elastic properties as the surrounding medium. If the tension center is heterogeneous with the rest of the medium, the integration can be performed only if the shape of V approaches a sphere, i.e. when the differences between its dimensions along 3 arbitrary axes are not much pronounced. The appIicab~ity of the above approach faiIs if appreciable contrasts in the parameters of the overlying and underlying rocks occur. One can resort to a finite element analysis, but of course the extent of the system analyzed will be perforce limited.

Natural hydrologic basins usually consist of unconsolidated sedimentary soils. In playered sediments it is customary to discard the horizontal components of displacement in that they are assumed to be small in comparison with the vertical deformation. Flow is still considered to be three-dimensional, Since the consolidation is exclusively vertical, relation (4) is consistent with the approximation involved in the ~ni-dimen~on~ compaction. Coefficient cr in this case is the vertical ~~mpressib~ity of a sample of soil prevented from exp~di~ laterally. The diffusion equation is solved numerically in a multi-aquifer-aquitad system under specified boundary conditions (natural conditions of the system) and imposing a prescribed discharge from the artificial wells. Any approach based on the diffusion equation will necessarily be a two-step procedure involving: a) a hydrologic model to caIcufate the time dependent hydraulic head field, b) a one-dimensional subsidence calculation based on these results. As far as the hydrologic model is concerned, there are at least three possibilities: 1) a full threedimensional analysis, 2) a quasi-three-dimensional analysis [ 191, 3) an analysis in a two-dimensional vertical cross section in Cartesian or radial coordinates. The first possibility must often be dismissed in complex naturaI basins on the grounds that there are serious computer limitations, in terms of both time and storage, to a full three-d~ension~ synthesis at this time. The quasithree-dimensional model is one of a class of models that links together a set of horizontal aquifer models by means of the leakage terms through the intervening aquitards. The fundamental assumption of this approach is that flow is horizontal in the aquifers and vertical in the aquitards. Neuman and Witherspoon f20] have stated that this assumption is satisfactory as long as the p~rmeab~ity contrasts are greater than two orders of ma~itude* The quasi-three-dimensional model of Pinder and Bredehoeft [ 191 belongs to this class of models, with leakage from the aquitards given by a flux calculation based on an analytical solution for flow in the aquitards. The subsidence is the sum of the compaction of the aquitards and its time-dependent value is provided automatically during each time step of the simulation.

G. Gambolati, Numerical models in land subsidence control

233

The primary purpose of a bi-dimensional flow model (either in radial or Cartesian coordinates) is to obtain head drawdown in the aquifers. These values are then inserted as boundary conditions in a refined one-dimensional vertical subsidence model that makes use of the local detailed data by considering the actual multilayer configurations within the aquitards. The model is first calibrated against the past history of available records (observed drawdowns and measured subsidence) by adjusting the regional formation parameters (permeability and compressibility) within their range of variation as provided by the laboratory test analyses. The model is then used as a predictive tool under a variety of possible pumping schedules. At Venice we have opted for an axi symmetrical model: a short description involving theory and results is given in the section below.

3. Two illustrative examples The applications recently made [ 21, [ 2 11 for subsidence prediction in Venice (water pumping) and in the Groningen area (gas production) are representative of the modeling advances to date and of the different approaches that must be used in hydrologic basins and gas-oil fields. At Venice [ 221 we have opted for a two-dimensional vertical cross-section in radial coordinates. It is centered in Marghera (fig. 2a), where the major withdrawals occur, and lies along a line from Mar&era to the Lido through Venice (fig. 2a). The base of the section is located at 300 m depth. The horizontal dimension is 20 km long. Fig. 2b is a schematic diagram of the detailed stratigraphy in the upper 300 m, and fig. 2c shows the nodal mesh used in the simulation performed with a finite element technique. The forcing function is the time-dependent withdrawal rate at Marghera. The calibration of the model was done with a trail and error procedure in which various possible combinations of values (in keeping with the ranges of values suggested by the rather sparse available data) were tested to see which provided the best fit of the past observed records of hydraulic head decline at Venice. In the second step of the procedure the hydraulic head declines that were determined with the calibrated hydrologic model were used to compute the land subsidence (subsidence model). Fig. 3a shows the output from the hydrologic model for the 100 m deep aquifer [ 2 1I. The calibration portion on the left shows the hydraulic head drop at Venice for the period 1930-l 973 as provided by the mathematical model and as calibrated against the available records. The left-hand portion of fig. 3b shows the computed subsidence at Venice. For the period 1930- 1973 a value of 15 cm is indicated. The right-hand portions of the two diagrams of fig. 3 show the predicted changes in hydraulic head and subsidence (or rebound) under two possible pumping schedules. Curves 1 and 2 refer respectively to constant consumption from 1974 onwards (as it has been since 1969) and to shutdown of all wells in 1974. In the first case about 3 cm of further subsidence is to be expected. In the second there would be a modest rebound of perhaps 2 cm in the next 25 years. Artificially recharging the depleted aquifers would increase the rate of rebound but not the ultimate subsidence. Over 85% of the Venice subsidence is irreversible. The model used at Groningen [ 181 is totally different, in accordance with the different geological setting existing there. The reservoir discovered in 1959 began to be exploited in 1963. Its depth of burial is approximately 2900 m and its extent about 900 km2 (fig. 4a). The thickness ranges from 70 m to 240 m. In the initial period of production until 1970 no appreciable increase

G. Gambukzti, iVumeriral models in lami subsidence control

234

bl

a) 0

t5dzu&2

MhRGHERA

VENICE

LID0 0

100 s z zoo

300

Fig. 2. (a) Location map of Venice. (b) Schematic stratigraphic section across the Venetian Lagoon, O-300 m depth. (c) Nodal mesh used in the fmite element model. The location of the aquifers is indicated.

CsMvetion IO,

1930

(a)

50

70

90

2000

1930

50

70

90

2000

U-9

Fig. 3. (a) Calibrated and predicted hydraulic head drawdown at Venice in the aquifer at 100 m depth. (b) Calibrated and predicted subsidence at Venice. The open circles represent the observed records.

z c

G. Gombolati. Numerical models in land subsidence control

235

lction

(a)

01 1975

2000

2025

2050

2075

2100

Cc)

Fig. 4. (a) Location map of the Groningen area. (b) Discretization of the reservoir volume used in the model based on the strain nuclei approach. (c) Maximum predicted subsidence from 1975 to 2100 (after Geertsma and van Opstal [ 181).

of the subsidence with respect to the rate prior to the field development has been measured. The simulation has therefore been used only to predict the future lowering of the area. The authors [ 181 believe that there is sufficient evidence that the underlying and overlying formations have

236

G. Gambolati, Numerical models in land subsidence control

the same elastic properties as the reservoir and apply the strain nuclei approach described earlier. The pressure distribution forecasts have been calculated independently using a specific computer program based on the solution of the diffusion equation under assumed conditions for the field development. The reservoir volume has been discretized in small portions (fig. 4b), the volume and depth of each portion have been estimated and a prediction of the subsidence from 1975to 2 100 has been provided using relation (9). Fig. 4c shows the maximum values of the expected subsidence. When the reservoir is totally depleted and a new equilibrium is reached, a final value of 1 m is indicated. For both cases - Venice and the Groningen area - the subsidence situation should remain under continued evaluation as further data become available. These will allow for a more reliable calibration of the models and therefore adjustments in the prediction.

Acknowledgement The author is indebted to R.A. Freeze from whom he learned the importance of flow analysis in real subsidence problems and with whom he shared the fascinating experience of simulating and predicting the subsidence of Venice.

References [l] J.F. Poland and G.H. Davis, Land subsidence due to withdrawal of fluids, Rev. Eng. Goel. 2 (1969) 187-269. [ 21 Nederlandse Aardolie Maatschappij B.V. (ed.), The analysis of surface subsidence from gas production in the Groningen area, the Netherlands, Verhandelingen Kon. Ned. Geol. Mijn. Gen. 28 (1973) pp. 109. [ 3) M. Caputo, G. Folloni, A. Gubellini, L. Pieri and M. Unguendoli, Survey and geometric analysis of the phenomenon of subsidence in the region of Venice and its hinterland, TR. No. 9, CNR, Laboratorio per lo studio della dinamica delle grandi masse, Venice (1971) pp 20. [4] J.H. Green, The effect of artesian-pressure decline on confined aquifer systems and its relation to land subsidence, Paper 1779-T, U.S. Geol. Survey (1964) pp. 11. [5] B.E. Lofgren and R.L. Klausing, Land subsidence due to groundwater withdrawal, Tulare-Wasco area, California, Paper 437-B, U.S. Geol. Survey (1969) pp. 101. [6] G. Gambolati, Estimate of subsidence in Venice using a one-dimensional model of the subsoil, IBM J. Res. Develop. 16 (1972) 130-137. [7] D.R. Allen, The mechanics of compaction and rebound. Wilmington oil field, Long Beach, California. Tee. Rep. of the Dep. of Oil Properties, City of Long Beach, California (1969) pp. 47. [8] P.A. Domenico and M.D. Mifflin, Water from low-permeability sediments and land subsidence, Water Resour. Res. (1965) 5633576. [9] J.F. Poland, Status of present knowledge and needs for additional research on compaction of aquifer systems, in: Land Subsidence, Proc. IASH 1 (1969)11-21. [lo] R.S. Sandhu and E.L. Wilson, Finite element analysis of land subsidence, in: Land Subsidence, Proc. IASH 2 (1969) 393-400. [ 1 l] R.L. Schiffman, A.T.F. Chen and J.C. Jordan, An analysis of consolidation theories, I. Soil Mech. Found. Div. ASCE 65 (1969) 285-312. [12] M.A. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (1941) 155-164. [13] M.A. Biot, Theory of elasticity and consolidation for a porous anisotropic solid, J. Appl. Phys. 26 (1955) 182-185. [ 141 A. Verruijt, Elastic storage of aquifers, in: Flow through porous media, ed. R.J.M. Dewiest (Academic Press., New York, 1969) 331-376. [15] J. Ghaboussi and E.L. Wilson, Flow of compressible fluid in porous elastic media, Int. J. Numer. Meths. Eng. 5 (1973) 419-442. [ 161 J. Geertsma, Problems of rock mechanics in petroleum production engineering, Proc. 1st Cong. Sot. Rock Mechan., Lisbon (1966) 585-594.

G. Gambolati. Numerical models in land subsidence control [ 171 G. Gambolati, A three-dimensional

231

model to compute land subsidence, U&H Bull. 17 (1972) 219-226.

[ 181 J. Geertsma and H. van Opstal, A numerical technique for predicting subsidence above compacting reservoirs, based on the nucleus of strain concept, see ref. [2], (1973) 43-62.

[ 191 G.F. Pmder and J.D. Bredehoeft, Application of the digital computer for aquifer evaluation, Water Resour. Res. 4 (1968) pp. 1069-1093. [20] S.P. Neuman and P.A. Witherspoon, Applicability of current theories of flow in leaky aquifers, Water Resour. Res. 5 (1969) pp. 817-829. [21] G. Gambo1ati.P. Gatto and R.A. Freeze, Mathematical simulation of the subsidence of Venice, 2. Results, Water Resour. Res. 10 (1974) 563-577. [ 221 G. Gambolati and R.A. Freeze, Mathematical simulation of teh subsidence of Venice, 1. Theory, Water Resour. Res. 9 (1973) pp. 721-733.