- Email: [email protected]
The effect of dispersion on the establishment of a paleoclimatic record from groundwater
Journal of Hydrology, 58 (1982) 131--147
131
Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands [2]
THE EFFECT OF DISPERSION ON THE ESTABLISHMENT OF A PALEOCLIMATIC RECORD FROM GROUNDWATER
M.R. DAVIDSON and P.L. AIREY
Australian Atomic Energy Commission Research Establishment, Sutherland, N.S. W. 2232 (Australia) (Received February 25, 1981;revised and accepted August 24, 1981)
ABSTRACT
Davidson, M.R. and Airey, P.L., 1982. The effect of dispersion on the establishment of a paleoclimatic record from groundwater. J. Hydrol., 58: 131--147. A mathematical model which can account for the dispersion of solute in underground water over time scales exceeding 104 yr. is described. In its simplest form, the model postulates a diffusion limited exchange of material between flowing groundwater and the stationary water in the confining aquiclude, leading to the spreading of the solute in the mobile water and a consequent reduction in its concentration variability. The model is used to study the increasing attenuation with "age" of the perturbations in groundwater composition induced by paleoclimatic oscillations in the late Quaternary. Specific reference is made to noble-gas concentrations, stable-isotope ratios and chloride levels. The extreme attenuation of stable-isotope ratios in the Australian Great Artesian Basin cannot yet be adequately explained and may provide the basis of a critical evaluation of different dispersion models in large hydraulic basins.
INTRODUCTION
The possible impact of widespread industrialisation and land-use change on future climatic trends is becoming increasingly widely recognized (Galbally and Freney, 1978). Assuming that the extreme of conditions to which the earth has been exposed in the past will define the limits of change in the foreseeable future, there is substantial practical value in assembling continuous records of paleoclimatic variability. Useful information can be obtained under the following conditions: (1) Correctly interpreted climate sensitive parameters can be measured with adequate precision. (2) A sequence has been preserved which can be related to a satisfactory chronological model. (3) Adequate time resolution is possible. Climatic changes cause perturbations in the chemical and isotopic composition of components of the water cycle; the effects have been monitored in
0022-1694/82/0000--0000/$02.75 © 1982 Elsevier Scientific Publishing Company
132 ice cores (Dansgaard et al., 1969; Delmas and Ascenia, 1980) and in marine sediments {Hays et al., 1976}. Continuous records can, in principle, also be obtained from groundwater. Provided that there are no serious complications due to aquifer mixing or subsurface chemistry, the isotopic composition of water and some aspects of the chemistry could reflect conditions in the recharge area over a period commensurate with the residence time of the groundwater. The parameters which have been correlated with climatic trends include: (1) Stable-isotope ratio variations (Bath et al., 1979; Deak, 1979; Evans et al., 1979; Shampine et al., 1979). (2) Variations in the concentrations of dissolved noble gases (Mazor, 1972; Andrews and Lee, 1979; Bath et al., 1979). (3) The chloride composition of some groundwaters in sandstone aquifers (Airey et al., 1979; Bath et al., 1979). The advantages of the groundwater systems are that the samples can be readily collected, are representative of the area and can be dated using hydraulic or isotopic data. The disadvantages are that detailed interpretation requires an understanding of the effect of climate on the complex processes leading ultimately to precipitation and recharge. In addition, resolution of the data depends not only on the density of the sampling wells but also on the dispersion processes within the host material. The dominant paleoclimatic oscillations have been assessed to have periods of a b o u t 23,000, 42,000 and 100,000 yr. (Hays et al., 1976}. Although the relevant groundwater properties do not, in general, respond to climatic parameters in a simple way, it is assumed that periodic variations will occur which can be analysed into combinations of the same frequencies. Investigation of dispersion phenomena which may attenuate these large-scale oscillations i s important in the interpretation of observed groundwater records. Evidence for such attenuation may come from the study of the Australian Great Artesian Basin, for which groundwater residence times exceeding 3" l 0 s yr. (Habermehl, 1980) have been calculated. There is already an extensive b o d y of literature devoted to the description of transport processes in porous media. The dominating large-scale mechanical dispersion in geologic systems is believed to be produced b y gross changes in local flow directions caused b y non-uniformities in hydraulic conductivities (Smith and Schwartz, 1980). A medium length associated with this nonuniformity is the dispersivity D L / V where D L is the longitudinal mechanical dispersion coefficient and V is the pore-water velocity. {See also Notation for symbols used in this paper.) Typical dispersivity values for aquifer systems have been reported by Lenda and Zuber {1970) to range from a few centimetres in sands to a few metres for coarse gravels or for large dispersion distances. The value chosen in this work (3 m) is in the upper range consistent with the large distances and long time scales being modelled, and the choice of V (3 m yr. -1 ) is representative of groundwater velocities in the Australian Great Artesian Basin (Habermehl, 1980). Note that even on large
133 NOTATION List of symbols a
Ao,A(x) b C l ( X , y , t ) , C1 C2 (x, y, t), C2 D DL DT F(t) K l L t, A t T V W x, ~ x X y, Ay ~1, P2
T
parameter (p2D 1/2 x l p l VL ) amplitude of the sinusoidal oscillation at input and at the position x parameter (lID i n ) solute concentration in aquifer solute concentration in aquiclude thickness averaged concentration in aquifer at position x and time t effective diffusion coefficient (aquiclude) scaled for tortuosity longitudinal mechanical dispersion coefficient (aquifer) transverse mechanical dispersion coefficient (aquifer) thickness averaged solute concentration at input (aquifer) dimensionless parameter (P2 I/Pl L) aquiclude thickness aquifer thickness time, and time increment, respectively period of oscillation associated with the sinusoidal input groundwater pore velocity [(Darcy velocity)/ (effective porosity) ] 2~ times the frequency of oscillation associated with the sinusoidal input longitudinal distance from the recharge zone, and distance increment dimensionless distance parameter ( Dx l Vl 2 ) distance associated with the thickness of the aquifer, and distance increment effective porosity in aquifer and aquiclude, respectively dimensionless time (Dt/l 2 )
(mol 1-1 ) (mol 1-1 ) (mol 1-1 ) (mol 1-1 ) (m 2 yr. -1 ) (m 2 yr. -1 ) (m 2 yr. -1 ) (mol 1-1 ) (m) (m) (yr.) (yr.) (m yr. -1 ) (Ty'l".-I)
(m) (m)
t i m e scales, t h e c h o i c e o f a c o n s t a n t value f o r DL m u s t be v i e w e d w i t h c a u t i o n since t h e m a g n i t u d e o f t h e d i s p e r s i o n m a y increase w i t h o u t limit as t h e geological size o f t h e s y s t e m increases b e c a u s e o f t h e successively larger scale o f n o n - u n i f o r m i t y w h i c h m a y b e e n c o u n t e r e d ( S m i t h a n d S c h w a r t z ,
1980). T h e d i s t a n c e t r a v e r s e d b y a m o l e c u l e o f s o l u t e u n d e r t h e a c t i o n o f longit u d i n a l m e c h a n i c a l d i s p e r s i o n in elapsed t i m e t since g r o u n d w a t e r recharge, is ( D L t ) i n . P r e d i c t e d values o f DL are o r d e r s o f m a g n i t u d e g r e a t e r t h a n t h o s e a n t i c i p a t e d f r o m l a b o r a t o r y studies; nevertheless, l o n g i t u d i n a l m e c h a n ical d i s p e r s i o n e f f e c t s c a n b e i g n o r e d w h e n t h e length scale a l o n g t h e aquifer, o v e r w h i c h c o n c e n t r a t i o n d i f f e r e n c e s are t o b e investigated, is m u c h g r e a t e r t h a n t h e d i s p e r s i o n d i s t a n c e (D L t) 1/2 . F o r e x a m p l e , t h e length scale f o r a
134 cyclic input of period T is VT; longitudinal mechanical dispersion can then be ignored when VT >~(DLt) 1/2 . This condition is met when the period T is sufficiently large. For a 104-yr. climatic cycle, and groundwater residence times t in the range 104--106 yr., VT is one to two orders of magnitude greater than ( D L t) I/2 , based o n D L / V ~- 3 m and V = 3 m yr. -1 . Note that the upper value chosen for t (106 yr.) far exceeds typical groundwater residence times, so that there is never sufficient time available for longitudinal mechanical dispersion to significantly influence the large-scale paleoclimatic oscillations. A different dispersion mechanism which may be more influential at long time scales is discussed below. The corresponding (analytical) model approximates conditions f o u n d in m a n y large basins and can be readily applied to very old groundwater. In its simplest form, it ignores longitudinal mechanical dispersion effects (given by D L ) but postulates a molecular diffusion limited exchange of material between the flowing groundwater and the stationary water in the confining aquiclude, leading to a spreading of the solute in the mobile water and a consequent reduction in the concentration variability from that existing at recharge. A characteristic time for a molecule of solute to diffuse across an aquiclude of thickness I is 12/D where D is the effective diffusion coefficient of the solute therein. Solute diffusion coefficients in water are of the order of 10 -s cm 2 s-1 which gives a value for D of 0.01 m 2 yr. -1 after allowing for a tortuosity factor of 3. With this value o f D , 12/D increases from 104 yr. when I = 10 m to ~ 1 0 s yr. when l = 3 0 m , suggesting that dispersion, which results from the interaction of axial flow in the aquifer and transverse molecular diffusion in the aquiclude, can be important on time scales similar to the periods of the major climatic cycles. An analogous model has been used by Turner {1958) and Skopp and Warrick (1974) to represent a porous medium with an internally distributed stationary phase (e.g., water in the boundary region of pore spaces) for which l is very small. Other theoretical considerations involving stratified media have been made by Shamir and Harleman'(1967), Fried and Combarnous (1971), and A1-Niami and R u s h t o n (1979). In this paper, the model is used to interpret aspects of the variations of stable-isotope ratios, inert-gas concentrations and chloride levels in groundwater. An assessment is made of the factors involved in interpreting the spatial variability of chosen groundwater properties in terms of paleoclimatic oscillations at recharge.
MATHEMATICAL
ASPECTS
Problem definition Consider an idealized homogeneous isotropic semi-infinite aquifer of height L with an adjacent confining aquiclude of thickness l which is also
135
stationary
aquic~ude aquifer
L
mobite
PV
Fig. 1. A schematic representation of a double-layer aquifer--aquiclude system. Water in the aquifer is moving longitudinally with pore velocity V and water in the aquiclude is stationary.
homogeneous and isotropic. This layered system is represented schematically in Fig. 1. Attention is focused on the saturated zone and it is assumed that no water reaches the aquiclude from above; all water is presumed to enter the system at the recharge zone (x -- 0). Flow in the aquifer is taken to be one-dimensional with a constant pore velocity V [= (Darcy velocity)/(effective porosity)] while water contained in the aquiclude is stationary. This section is concerned with the mathematical aspects of longitudinal convection and dispersion in the aquifer, from the recharge zone, of a nonreacting solute which can diffuse into and o u t of the adjacent aquiclude. The effects of radioactive decay of the solute and of adsorption onto the solid matrix of the system are ignored here b u t they can be easily included in the model. Solute concentrations are assumed to be low so velocity components induced b y mass transfer between layers can be ignored. The convection--diffusion equations for the solute in terms of rectangular coordinates ( x , y ) and time t are:
0C 1/Ot -k V(OC 1/Ox)
-- D w (02 C 1 / 0 y 2 )
(1)
and
0C2/0t = D(0 2 c 2 / 0 y 2)
(2)
where axial mechanical dispersion in the aquifer and axial diffusion in the aquiclude have b o t h been ignored (i.e. terms involving 0 2 C1/Ox 2 and 0 2 C2 / 0x 2 are not included). Henceforth, the process governed by eqs. 1 and 2 will be called double-layer dispersion. Solute concentrations in aquifer and aquiclude water are given by C1 (x, y, t) and C2 (x, y, t), respectively. D w is the transverse dispersion coefficient in the aquifer and D is the molecular diffusion coefficient (scaled for tortuosity) in the aquiclude. Boundary conditions describing continuity of solute concentration and flux at the interface y = L are: C1 (x, L, t) = C2 (x, L, t)
(3a)
and 0C1 8C2 plDw -~y ( x , L , t ) = IJ2D- : - ( x , L , t ) Oy
(3b)
136
where Pl and #2 are the porosities of the aquifer and the aquiclude, respectively. Zero-flux conditions ~C1/Oy = 0 and ~C2/~y = 0 are chosen at y = 0 and y = L + l, respectively. Taking the average of eq. 1 over the thickness L of the aquifer gives: DC~ + V -~C1 - - - •2 D -[~C2]/
~t
~x
(4)
Pl L \ ~y ly=L
after applying the boundary conditions at y = 0 and y = L, and defining L
Ca ~- L
C1 dy 0
The term on the right of eq. 4 represents the flux of solute passing between aquifer and aquiclude. Since diffusion in the aquiclude is m u c h slower than transverse dispersion in the aquifer (i.e. D ~ DT ), transverse gradients (OC1/Oy) in the aquifer will be m u c h smaller than those in the aquielude provided that the transverse concentration distribution in the aquifer is initially uniform or that sufficient time has elapsed for the ~C1/Oy to decay (on a time scale corresponding to diffusion in the aquiclude, transverse dispersion in the aquifer is instantaneous). Under these circumstances, C1 may be replaced by C1 in condition (3a) and the problem becomes identical to that described by Skopp and Warrick (1974) for which the Laplace transform solution has t h e form: Lap(C1 ) = Lap(F) exp
( sx V
x #2 (sD)lntanhl(s/D)ln) VL Pl
(5)
when Cj (x, 0) = C2 (x,y, 0) = 0, Ca (0, t) = F(t) and Lap denotes a Laplace transform (e.g., Lap(Ca ) = f g e-"tCl dt).
Step-function input (an unsteady problem) S k o p p and Warrick (1974) inverted eq. 5 for a step-function input [F(t) = 1] to give the solution: C1
2 of 7r
x
[ ~(sinhbX--sinbX)] exp -\cosh bX + cos bX
sin
t--
2
2 \e-osh ~ + -e-oosb-X])
7' and
fort>x/V
(6)
137
C1 = O, f o r t ~ x / V where
a = (p2Dinx)/(pl VL)
and
b =
liD 1/2 when Dt/l:
An independent analysis shows that ~ 1 (i.e. the elapsed time is much smaller than the characteristic transverse diffusion time for the aquiclude), this solution m a y be approximated b y that corresponding to an aquiclude of infinite thickness, viz.:
C1 = ,I o,erfc[½a/(t--x/V)U2]'
forf°r t~x/>x/V yt
(7)
where erfc is the c o m p l e m e n t a r y error function. In principle, it is n o w possible to determine the response to an arbitrary input function b y means of a convolution integral. However, in practice it is more convenient to use a numerical method. The details of this technique are given later (p. 10). Fig. 2 represents a comparison between double-layer (solid curves) and longitudinal mechanical (dashed curves) dispersion in the case of a step input. The solid curves are determined numerically for convenience (after verifying that agreement with the results of eqs. 6 and 7 is better than 1%) and the dashed curves are calculated from eq. 10.6.22 in Bear (1972) for
I.--0
0.1
I
05
8
I-0
~___ I-5
2"0
x=£_~ V.~2
c", LLI N N
\ \ \
0
6'0
I
I
7"0
8"0 Dx x=~-2
9'0
I0'0
Fig. 2. R e s p o n s e o f t h e m o d e l a q u i f e r t o a s t e p - f u n c t i o n i n p u t o f s o l u t e c o n c e n t r a t i o n . The solid curves illustrate the dispersion due solely to a diffusion limited exchange of s o l u t e b e t w e e n f l o w i n g a q u i f e r w a t e r a n d i m m o b i l e a q u i c l u d e w a t e r w h e n K = P2 l / P l L = 0.25. The dashed curves represent only macroscopic dispersion due to the non-uniformities in h y d r a u l i c c o n d u c t i v i t i e s f o r V---- 3 m yr. -1 , D L = 10 m 2 y r . -1 a n d are p l o t t e d a g a i n s t d i m e n s i o n l e s s v a r i a b l e s X = D x / V l 2 a n d T = D t / l 2 a s s u m i n g l = 10 m , D = 0.01 m 2 yr. -1 "
138 V = 3 m yr. -1 and D L = 10 m 2 yr. -1 and plotted against dimensionless distance (X = D x / V l 2 ) and dimensionless time (r = D t / l 2 ) in Fig. 2, assuming t h a t l = 10 m and D = 0.01 m E yr. -1 . Increasing l has the effect of steepening the dashed curves. We first note that double-layer dispersion becomes significant when T is of order one (i.e. t of order 12/D) for K [= (P21/PlL)] = 0.25 and t h a t this critical value of r varies inversely with the exchange capacity coefficient K (as was expected). The next important observation is that it is indeed appropriate to ignore longitudinal mechanical dispersion (regardless of a priori arguments about its importance or otherwise at long time scales) when it is compared with double-layer dispersion for K = 0.25 and r ~> 1, since the dashed curves become much steeper than the solid curves in Fig. 2. The choice of K = 0.25 in this model is based on the following parameter values: P2 = 0.05, Pl -- 0.20 and IlL = 1 (i.e. identical thickness for both aquifer and aquiclude). The choice of P2 assumes an aquiclude such as shale. Alternatively, for a clay aquiclude we might choose P2 = 0.5, Pl = 0.2 and IlL = 0.1 to again obtain K = 0.25.
Sinusoidal input (a steady-state problem) Turner (1958) considered the frequency response of a packed-bed model which included stagnant pockets. He obtained the solution for the response to a periodic input, after a steady cyclic state had been reached, directly from eqs. 2 and 4. When C1 (0, t) = A 0 sin Wt where W = 2~r/T and T is the period of the oscillation, Turner's solution becomes:
C1 = Ao e x p ( - - p x / V ) s i n ( W t - - q x / V )
(8)
for the case w i t h o u t longitudinal mechanical dispersion, where:
KO (sinh bO -- sin bO t P = 2--b-\cosh bO + cos bO] KO ( sink bO + sin bO q = W+-2b \ e o s h bO + cos bO !
and
0 = (2W) in
This result can also be obtained quite readily from eq. 5 by using the expansion theorem given by earslaw and Jaeger (1953, p. 280) for large time solutions. If A(x) is the amplitude of the oscillation as it proceeds along the aquifer, the amplitude attenuation factor is given by:
A/Ao
= exp(--px/V)
and the wave speed by WV/q. The attenuation factor is plotted against longitudinal distance in Fig. 3. Note that, as was expected, the higher the frequency of oscillation the more rapidly it damps out. In Fig. 4, the decay constant p is plotted against aquiclude thickness I for
139
¢: o
-5 IE
5
10
15
20
Dx X = Vt 2
Fig. 3. Relative amplitude attenuation with dimensionless distance X = D x / $ r l 2 of a sinefunction input concentration (period T) after a steady cyclic state has been reached throughout the model aquifer, where K = g2 l/P1L = 0.25 and T* = DT/I 2.
% × (3_
©
I
I
I
J
20
40
60
80
100
(m)
Fig. 4. Amplitude decay constant p for a periodic input concentration (period T) is plotted against aquiclude thickness I for fixed K = tl 2 I / p 1 L = 0.25.
t h e t h r e e d o m i n a n t climatic oscillations w h e n K is f i x e d at 0.25 a n d D at 0.01 m 2 yr. -1 . F o r a given p o r o s i t y ratio this c o r r e s p o n d s t o h o l d i n g IlL fixed. Each curve e x h i b i t s an o p t i m u m value o f I w h i c h achieves a m a x i m u m a m p l i t u d e d e c a y rate. These o p t i m u m values are 10, 14 and 20 m for t h e periods 2 3 , 0 0 0 , 4 2 , 0 0 0 a n d 1 0 0 , 0 0 0 yr., respectively.
The numerical m e t h o d A c o n v e n i e n t and rapid m e t h o d o f solving eqs. 2 and 4 for an a r b i t r a r y i n p u t f u n c t i o n is t h e n u m e r i c a l t e c h n i q u e d e s c r i b e d below. (1) D e f i n e n o d a l p o i n t s (i = 1 , . . . , I ) o f c o n s t a n t spacing Ax along the length o f t h e o n e - d i m e n s i o n a l m o d e l a q u i f e r and, at each such p o i n t , d e f i n e
140
nodes (j = 1, . . . , J ) of constant spacing Ay along the transverse (y-) direction. The objective is to solve eqs. 2 and 4 by replacing them with difference equations at each nodal point and then marching through time t in discrete steps (n = 1 , . . . , N) of length At. Let (C 1 )[' be the value of C1 at position x = izkx and (C2)i~,j be the value of C2 at x -- iAx and y -- L = j A y , all evaluated at time t = n a t . (2) Let A x = V A t and move the concentration distribution C1 forward a distance Ax by setting (C1)* = (C1)['-1 • This is equivalent to solving:
0C 1/0t "~ V(0C 1/Ox)
= 0
(9)
(3) Now solve eq. 2 by replacing it with the following difference equation:
[(C2 ~.n+l ,,,~ - - (C2)i,i]/At "
= D I { ( C 2 ) ¢ ,~+1 , s + l - - 2 ( C 2 ) i , 1 .+1
-[- (C2)n,~ll}/(Ay)2] (10)
where (C2 ~,~1 = (C1)* and (C2 )t,J ~,+1 = (C2 Ji,J-I ~,÷1 are the corresponding difference equations for the b o u n d a r y conditions. This equation is called "implicit" because the space derivative is taken at the new time (n + l ) A t rather than at n A t and has the property of numerical stability (i.e. errors are not amplified during the calculation). Its solution is a standard procedure (see, e.g., Richtmyer and Morton, 1967). (4) Finally derive (C1)7 +1 from the difference approximation to aC1
P2 D [ a C 2 ~
(10)
viz.
( e l ) n• + l - (el)~ At
]../2
D ((C2'n+l--(c2)n+l I) Ji,2
Pl L
(11)
Ay
The above procedure involves splitting eq. 4 into two parts, the first being pure convection and the second the exchange of solute between strata. The choice of Ax = V A t ensures t h a t the solution of eq. 9 at nodal points is exact so no artificial numerical dispersion is introduced. As mentioned earlier, results have been checked against the analytical solutions ( 6 ) a n d (7) for a step-function input. In this case, ~ 1% accuracy or better was obtained by choosing Ay = 0.025l and At at least 0.01 times the elapsed time t for D t / l 2 ~ 10. However, for an arbitrary input function, the errors incurred may be different so the sensitivity of the results to At and Ay should always be checked. DISCUSSION General c o m m e n t s
The model described above is applicable to m a n y large hydraulic basins in which well-defined aquifer layers are bounded by low-permeability aqui-
141
cludes. It has the following desirable features: (1) It is the first example of a model evaluated at long time scales which predicts attenuation in the cycle of groundwater solute concentrations which oscillate with periods exceeding 104 yr. The predicted dispersion is a consequence of the differences in transport properties within the aquifer and the confining aquiclude. The choice of K = 0.25 in the model is discussed in the Section " S t e p function input (an unsteady problem)". (2) Analytical solutions of some simple input functions can be obtained. The unsteady response to a step input may be considered as a model of the response to non-steady deviations of solute concentrations from the underlying periodic trend which may, in turn, be modelled by the steady-state cyclic response to a superposition of sine and cosine inputs. (3) The model can readily be extended to study the effects of dispersion of generalized input functions using a simple numerical solution technique. (4) Subject to symmetry restrictions, the analytical model can be extended to include multiple layers of aquifer and aquiclude. These restrictions can be removed when using the numerical model. (5) Without conceptual change, the model can be extended to account for adsorption onto the host material of either the aquifer or aquiclude, and the effects of radioactive decay. (6) It is possible that the model might predict a significant isotopic difference when the isotope effects associated with adsorption and molecular diffusion in the aquiclude are taken into account.
Applications Dispersion of dissolved noble gases in groundwater. Variations in the concentrations of noble gases in groundwater have recently been used to monitor paleotemperature changes in recharge areas. Rainfall has been shown to saturate with inert gases before infiltration (Mazor, 1972). However, because re-equilibration would normally occur during the lengthy recharge process, the temperature monitored is probably that in the vicinity of the water table. Inert-gas concentrations will only reflect climatic changes over long periods in the absence of significant losses in the saturated zones. Some experimental evidence for this assumption has been discussed by Mazor (1972) in his study of the Jordan Rift Valley, Israel. A groundwater piston-flow model is tacitly assumed in the normal approach, where gaseous concentrations are related to recharge temperatures through laboratory solubility data. The current model illustrates that the simplification of ignoring dispersion effects may not be justified, depending on the aquifer and aquiclude thicknesses. In the Section "Mathematical aspects", we indicated that when K = 0.25, the double-layer dispersion of a step input becomes significant after an elapsed time of 12/D and that the relative amplitude attenuation of the dominant climatic oscillatory modes may be calculated using Fig. 4. The value of D was taken to be 0.01 m 2 yr. -1 .
142 A study has been made of the Ar and Kr levels in the Bunter Sandstone groundwater in the English Midlands (Andrews and Lee, 1979; Bath et al., 1979). The schematic cross-section (Bath et al., 1979, fig. 1) indicates that l and L are "~200 m. In this circumstance, 12/D is 4" 106 yr. and the amplitude of the climatic modes is reduced by only 2--4% at an equivalent distance of 3" 104 yr. from the recharge area. Clearly, the Ar and Kr levels predicted from the dispersion model will be little changed over this distance owing to the large layer thickness chosen. Had the sequence comprised an assembly of layers with thicknesses of order 10 m, the effect of dispersion would have been significant. For a double layer with l = L = 10 m, the amplitude of a 23,000-yr. climatic cycle for example is reduced by 57% at a point 30,000 yr. from the recharge area and the critical time for significant dispersion of a step input is of the order 104 yr. (12/D). Extension of the model to multiple layers is possible as indicated above. There has been a t e n d e n c y in more recent times to recommend monitoring the Kr/Ne or Ar/Ne ratios rather than the absolute concentrations (Mazor, 1972; Phillips and Davis, 1979}. The solubilities of the lower-atomicweight inert gases are far less temperature-dependent. Experimentally the advantages lie in the fact that inert-gas ratios are less sensitive than absolute concentrations to changes accompanying the evolution of soil gas from air, and to any degassing that might occur in the saturated zone. It is therefore important to note t h a t the model can only be applied directly to absolute solute concentrations, because the predicted dispersion depends on the molecular diffusion coefficient in the aquiclude and therefore on the atomic weight of the gas. Different gases disperse to different extents, and simple expressions cannot be developed for their ratios.
Stable-isotope ratio variation in the Australian Great Artesian Basin. A study has recently been reported of aspects of the isotope hydrology of the principal Jurassic aquifer of the Australian Great Artesian Basin (Airey et al., 1979; Hartley, 1979). The " a g e " of the groundwater, which was calculated from the 1900 potentiometric surface and the horizontal hydraulic conductivities increased regularly with distance from the recharge area from less than 2" 104 to more than 3 • 10 s yr. During this portion of the late Quaternary there were two major ice ages and an interglacial period. Altogether, 69 wells were sampled along defined flow lines in the principal study area. Despite the climatic changes the D/H ratios were f o u n d to be remarkably constant. The mean value was --41.8°/0o, and the standard deviation of the population 1.3. Part of the variability is u n d o u b t e d l y due to experimental error which is conservatively assessed to be 0.7~/00 . The 51SO-values were 6.6 + 0.9~/00 . Because of the relatively large experimental error (~ 0.5°/0o ), 40 of the samples were re-run on equipment in which the assessed errors at the l a level varied from 0.03 to 0.080/00 . The average value of this set was - - 6 . 8 7 + 0 . 7 7 % . Since the 5D- and 51SO-values are consistent with the "meteoric line", 5D = 85180 + 10, it would be expected that A(SD)/A{51So) -
-
143 would be of the order of 8. Experimentally the variations in 5D and 6180 were comparable. If confirmed, this result would be very significant and will u n d o u b t e d l y help to distinguish between various transport models. Confirmation is required because the samples were stored (albeit in tightly stoppered glass bottles) for over a year between collection and the more precise measurements. This discussion will therefore be concerned principally with the constancy of the 5D-values. It resolves itself into two parts: (1) the input function; and (2) dispersion within the saturated zone. (1) Input function: There is a good linear correlation on the global scale b e t w e e n the stable-isotopic composition of rainfall samples and the mean annual surface temperature. Dansgaard (1964) has shown that: A(SD)/A(temperature) ~ 5.6
and
A(5 lSo)/A(temperature) ~ 0.69
Care must be taken when assuming that variations in the isotopic composition of groundwater simply reflect temperature changes in the recharge zone. In fact they respond to a range of global and local parameters. Hartley (1981) has gone so far as to suggest that: "there seems little possibility of relating climate to the isotopic composition of groundwater other than by empirical relationships developed from other paleoclimatic evidence" Nevertheless, a simple explanation based on temperature changes and the Dansgaard correlation proved adequate in a number of instances. Examples include recent studies of the isotope hydrology of a number of aquifers in Saudi Arabia (Shampine et al., 1979) and of thermal waters in the Budapest region, Hungary (Deak, 1979). On the other hand, there are many cases in which the paleotemperatures calculated from stable-isotope variations are substantially less than those assessed from other evidence. For instance, in the Bunter Sandstone study (Bath et al., 1979), the 5 lSO-values may be interpreted as indicating a change in cloud formation temperature of at least 1--2°C over the time scales involved (35,000 yr.). Surface temperature changes of from 5 to 6°C were estimated from gas solubilities and other evidence. Again the range in 5180values from the other British aquifer systems discussed b y Evans et al. (1979) is similar. From an assessment of the stable-isotope variability at a number of sites it can be reasonably expected that, over the time scales accessible to 14C, the 51So range would be at least 2~/00, and the 5D-values between a b o u t four and eight times this value. Thus, the range in 5D-values found in the Great Artesian Basin was a b o u t an order of magnitude less than expected. Before discussing the effects of dispersion of 5D variability, t w o possible explanations associated with the nature of the input function need to be considered. As was pointed o u t by Evans et al. (1979), the coefficient of variation of 5D with temperature for m o n t h l y rainfall samples in a particular region may
144 be much less than t h a t predicted from Dansgaard's global correlation. Such arguments cannot be directly used in this study because the 8D-values of the rainfall from the two available stations in the general recharge area (Charleville and Longreach) exhibit partial correlation coefficients with rainfall a m o u n t which are substantially greater than those with temperature (Hartley, 1979, 1981). It is, in principle, possible that various climatic factors could act in opposite directions to effectively cancel the expected 5D variability. This argum e n t is countered by the observation that a group of fourteen wells in South Australia which are recharged from the western extremities of the Great Artesian Basin are also very constant. They exhibit a mean of--48.1°/00 and a standard deviation of 1.3°/00 (Airey et al., 1979; Hartley, 1979). The probability that the combination of climatic factors in two recharge areas separated by 1200 km would affect the variation of D/H ratios in the groundwater in the same way must be very small. The explanation for the observed constancy in the 5D-values must be associated with dispersion over the long residence times. (2) Dispersion effects. The Great Artesian Basin (Habermehl, 1980) extends to distances equivalent to ~ 3 . 1 0 s yr. from the recharge area; it is therefore appropriate to consider the 10 s -yr. climatic cycle (Hays et al., 1976). The amplitude of this cycle will attenuate more slowly than the other more rapid climatic oscillations and is considered in relation to the observed reduction in D/H variability. When l = L = 100 m, for example, calculations based on Fig. 4 indicate an amplitude reduction, for the l 0 s -yr. cycle, of 25% at an equivalent distance of 2" l 0 s yr. from the recharge area. This reduction increases to 46% when l = L = 10 m is chosen. Apparently the model cannot fully explain the gross smoothing of the 5D-values in the Great Artesian Basin. Further work is required, because it appears that the attenuation of stable-isotope ratio variability may provide a critical test of solute dispersion in very large hydraulic basins.
Groundwater chloride. Chloride concentration at recharge depends on a number of factors including: (1) The rate of input in the recharge area by chloride in rainfall and dry fall out; (2) The rate of input from the dissolution of aquifer rocks and sediments; (3) The rate of loss due to surface run off; (4) The net groundwater discharge. Because of the complexity of the overall process, it is not possible to define the principal climate determining parameter. Nevertheless, it has been shown in studies of the Jurassic sandstone aquifer of the Australian Great Artesian Basin (Airey et al., 1979) and in the Bunter Sandstone aquifer, England (Bath et al., 1979) that the plots of chloride concentration with groundwater age fall on well-defined curves with a qualitative similarity. In both
145 cases the levels decrease regularly with age from relatively high levels characterizing the recharge water. In the Great Artesian Basin study there is evidence for a secondary m a x i m u m during the last interglacial, ~ 1 . 1 . l 0 s yr. ago. The current model can be used to assess the extent of chloride ion dispersion over this time scale and, in principle, to provide a better measure of the input function by a systematic trial-and-error procedure, using the numerical model. This could be important when studying the variations in the chloride balance in the recharge area t h r o u g h o u t the late Quaternary. Information on the rainfall variation over the past 1.2" l 0 s yr. is already available (Kershaw, 1978) from a station in north Queensland. At this stage no a t t e m p t has been made to establish these balances.
CONCLUSIONS (1) The dispersion process described here can significantly reduce the variability of solute concentrations within an aquifer to below that present in recharging water over time scales which depend on the aquifer and aquiclude thickness. When K - - 0 . 2 5 , the dispersion of a step input becomes significant after an elapsed time of 12/D (e.g., 104 yr. when l = 10 m and D = 0.01 m 2 yr. -I ) and the relative amplitude attenuation of d o m i n a n t steady-state climatic oscillations m a y be calculated using Fig. 4. Such effects need to be taken into account when interpreting noble-gas concentrations, stable-isotope ratios and chloride levels in terms of paleoclimatic parameters. (2) The extent of dispersion predicted by the present model is not sufficient to explain the observed reduction in variability of stable-isotope ratios in the very old groundwater of the Great Artesian Basin. Further work is required as this observation m a y form the basis of a m e t h o d for the critical comparison of transport models in very large aquifer systems. (3) The numerical solution described here is both fast and reliable enabling the easy determination of the response of an aquifer--aquiclude system to an arbitrary input function from an arbitrary initial state. The technique m a y be used to extend the model to multiple aquifer--aquiclude layers.
ACKNOWLEDGEMENTS The modelling ideas in this work owe much to some stimulating discussions with Peter Hartley of the A.A.E.C. The authors are also indebted to Professor A.T. Wilson of the University of Waikato, New Zea 1~ .d, for the high precision 1so--160 mass spectrometry, and to Mr. D. R o m a , (A.A.E.C.) for preparing the Great Artesian Basin samples.
146 REFERENCES Airey, P.L., Calf, G.E., Campbell, B.L., Habermehl, M.A., Hartley, P.E. and Roman, D., 1979. Aspects of the isotope hydrology of the Great Artesian Basin, Australia. In: Proc. Symp. on Isotope Hydrology, 1978, Vol. I, Neuherberg, Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 205--219. Airey, P.L., Calf, G.E., Hartley, P.E. and Roman, D., 1980. Aspects of the isotope hydrology of two sandstone aquifers in arid Australia. In: I.A.E.A. Adv. Gr. on Application of Isotope Techniques to Arid Zones Hydrology, Vienna, Nov. 1978. Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 93--111. A1-Niami, A.N.S. and Rushton, K.R., 1979. Dispersion in stratified porous media: analytical solutions. Water Resour. Res., 15: 1044--1048. Andrews, J.N. and Lee, D.J., 1979. Inert gases in groundwater from the Bunter Sandstone of England as indicators of age and paleoclimatic trends. J. Hydrol., 41: 233--252. Bath, A.H., Edmunds, W.M. and Andrews, J.N., 1979. Paleoclimatic trends deduced from the hydro-geochemistry of a triassic sandstone aquifer, U.K. In: Proc. Symp. on Isotope Hydrology, 1978, Vol. II, Neuherberg, Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 545--568. Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier, New York, N.Y., 764 pp. Carslaw, H.S. and Jaeger, J.C., 1953. Operational Methods in Applied Mathematics. Oxford University Press, London, 2nd ed., 359 pp. Dansgaard, W., 1964. Stable isotopes in precipitation. Tellus, 16: 436--468. Dansgaard, W., Johnsen, S.J., M~bller, J. and Langway Jr., C.C., 1969. One thousand centuries of climatic record from Cape Century on the Greenland ice sheet. Science, 166: 377--381. Deak, J., 1979. Radiocarbon dating of thermal waters in the Budapest area. In: Int. Congr. on Isotopes in Nature, Leipzig, Pap. AUS-8-1/80. Delmas, R.T. and Ascenia, J.M., 1980. Polar ice evidence that atmospheric CO 2 20,000 year B.P. was 50 per cent of present. Nature (London), 284: 155--156. Evans, G.V., Otlet, R.L., Downing, R.A., Monkhouse, R.A. and Rae, G., 1979. Some problems in the interpretation of isotope measurements in United Kingdom aquifers. In: Proc. Symp. on Isotope Hydrology, 1978, Vol. II, Neuherberg, Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 679--708. Fried, J.J. and Combarnous, M.A., 1971. Dispersion in porous media. Adv. Hydrosci., 7 : 169--282. Galbally, I.E. and Freney, J.R., 1978. Modification of climate. In: A.B. Pittock, L.A. Frakes, D. Jenssen, J.A. Petersen and Z.W. Zillman (Editors), Climate Change and Variability - - A Southern Perspective. Cambridge University Press, Cambridge, pp. 269--282. Habermehl, M.A., 1980. The Great Artesian Basin, Australia. Aust. Bur. Miner. Resour., J. Geol. Geophys., 5: 9--38. Hartley, P.E., 1979. Deuterium/hydrogen ratios in Australian rain and groundwater. M.Sc. Thesis, University of New South Wales, Kensington, N.S.W. Hartley, P.E., 1981. Deuterium/hydrogen ratios in Australian rainfall. J. Hydrol., 50: 217--229. Hays, J.D., Imbrie, J. and Shakleton, N.J., 1976. Variations in the earth's orbit; pacemaker of the ice ages. Science, 194: 1121--1132. Kershaw, A.P., 1978. Record of last interglacial--glacial cycle from northeastern Queensland. Nature (London), 272: 159. Lenda, A. and Zuber, A., 1970. Tracer dispersion in groundwater experiments. In: Proc. Symp. on Isotope Hydrology, 1970. Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 619--637.
147 Mazor, E., 1972. Paleotemperatures and other hydrological parameters deduced from noble gases dissolved in groundwaters, Jordan Rift Valley, Israel. Geochim. Cosmochim. Acta, 36: 1321--1336. Phillips, F.M. and Davis, S.N., 1979. Noble gases in groundwater as paleoclimatic indicators. EOS (Trans. Am. Geophys. Union), 6 0 : 8 3 0 (abstract). Richtmyer, R.D. and Morton, K.W., 1967. Difference Methods for Initial-Value Problems. Interscience, New York, N.Y., 2nd ed., 405 pp. Shamir, U.Y. and Harleman, D.R.F., 1967. Dispersion in layered porous media. J. Hydraul. Div., Proc. Am. Soc. Civ. Eng., 93(HY5): 237--260. Shampine, J.W., Dinner, T. and Noory, M., 1979. An evaluation of isotope concentrations in the groundwater of Saudi Arabia. In: Proc. Symp. on Isotope Hydrology, 1978, Vol. II, Neuherberg, Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 443--463. Skopp, J. and Warrick, A.W., 1974. A two-phase model for the miscible displacement of reactive solutes in soils. Soil Sci. Soc. Am. Proc., 38: 525--550. Smith, L. and Schwartz, F.W., 1980. Mass transport, I. A stochastic analysis of macroscopic dispersion. Water Resour. Res., 16: 303--313. Turner, G.A., 1958. The flow-structure in packed beds -- a theoretical investigation utilizing frequency response. Chem. Eng. Sci., 7: 156--165.
131
Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands [2]
THE EFFECT OF DISPERSION ON THE ESTABLISHMENT OF A PALEOCLIMATIC RECORD FROM GROUNDWATER
M.R. DAVIDSON and P.L. AIREY
Australian Atomic Energy Commission Research Establishment, Sutherland, N.S. W. 2232 (Australia) (Received February 25, 1981;revised and accepted August 24, 1981)
ABSTRACT
Davidson, M.R. and Airey, P.L., 1982. The effect of dispersion on the establishment of a paleoclimatic record from groundwater. J. Hydrol., 58: 131--147. A mathematical model which can account for the dispersion of solute in underground water over time scales exceeding 104 yr. is described. In its simplest form, the model postulates a diffusion limited exchange of material between flowing groundwater and the stationary water in the confining aquiclude, leading to the spreading of the solute in the mobile water and a consequent reduction in its concentration variability. The model is used to study the increasing attenuation with "age" of the perturbations in groundwater composition induced by paleoclimatic oscillations in the late Quaternary. Specific reference is made to noble-gas concentrations, stable-isotope ratios and chloride levels. The extreme attenuation of stable-isotope ratios in the Australian Great Artesian Basin cannot yet be adequately explained and may provide the basis of a critical evaluation of different dispersion models in large hydraulic basins.
INTRODUCTION
The possible impact of widespread industrialisation and land-use change on future climatic trends is becoming increasingly widely recognized (Galbally and Freney, 1978). Assuming that the extreme of conditions to which the earth has been exposed in the past will define the limits of change in the foreseeable future, there is substantial practical value in assembling continuous records of paleoclimatic variability. Useful information can be obtained under the following conditions: (1) Correctly interpreted climate sensitive parameters can be measured with adequate precision. (2) A sequence has been preserved which can be related to a satisfactory chronological model. (3) Adequate time resolution is possible. Climatic changes cause perturbations in the chemical and isotopic composition of components of the water cycle; the effects have been monitored in
0022-1694/82/0000--0000/$02.75 © 1982 Elsevier Scientific Publishing Company
132 ice cores (Dansgaard et al., 1969; Delmas and Ascenia, 1980) and in marine sediments {Hays et al., 1976}. Continuous records can, in principle, also be obtained from groundwater. Provided that there are no serious complications due to aquifer mixing or subsurface chemistry, the isotopic composition of water and some aspects of the chemistry could reflect conditions in the recharge area over a period commensurate with the residence time of the groundwater. The parameters which have been correlated with climatic trends include: (1) Stable-isotope ratio variations (Bath et al., 1979; Deak, 1979; Evans et al., 1979; Shampine et al., 1979). (2) Variations in the concentrations of dissolved noble gases (Mazor, 1972; Andrews and Lee, 1979; Bath et al., 1979). (3) The chloride composition of some groundwaters in sandstone aquifers (Airey et al., 1979; Bath et al., 1979). The advantages of the groundwater systems are that the samples can be readily collected, are representative of the area and can be dated using hydraulic or isotopic data. The disadvantages are that detailed interpretation requires an understanding of the effect of climate on the complex processes leading ultimately to precipitation and recharge. In addition, resolution of the data depends not only on the density of the sampling wells but also on the dispersion processes within the host material. The dominant paleoclimatic oscillations have been assessed to have periods of a b o u t 23,000, 42,000 and 100,000 yr. (Hays et al., 1976}. Although the relevant groundwater properties do not, in general, respond to climatic parameters in a simple way, it is assumed that periodic variations will occur which can be analysed into combinations of the same frequencies. Investigation of dispersion phenomena which may attenuate these large-scale oscillations i s important in the interpretation of observed groundwater records. Evidence for such attenuation may come from the study of the Australian Great Artesian Basin, for which groundwater residence times exceeding 3" l 0 s yr. (Habermehl, 1980) have been calculated. There is already an extensive b o d y of literature devoted to the description of transport processes in porous media. The dominating large-scale mechanical dispersion in geologic systems is believed to be produced b y gross changes in local flow directions caused b y non-uniformities in hydraulic conductivities (Smith and Schwartz, 1980). A medium length associated with this nonuniformity is the dispersivity D L / V where D L is the longitudinal mechanical dispersion coefficient and V is the pore-water velocity. {See also Notation for symbols used in this paper.) Typical dispersivity values for aquifer systems have been reported by Lenda and Zuber {1970) to range from a few centimetres in sands to a few metres for coarse gravels or for large dispersion distances. The value chosen in this work (3 m) is in the upper range consistent with the large distances and long time scales being modelled, and the choice of V (3 m yr. -1 ) is representative of groundwater velocities in the Australian Great Artesian Basin (Habermehl, 1980). Note that even on large
133 NOTATION List of symbols a
Ao,A(x) b C l ( X , y , t ) , C1 C2 (x, y, t), C2 D DL DT F(t) K l L t, A t T V W x, ~ x X y, Ay ~1, P2
T
parameter (p2D 1/2 x l p l VL ) amplitude of the sinusoidal oscillation at input and at the position x parameter (lID i n ) solute concentration in aquifer solute concentration in aquiclude thickness averaged concentration in aquifer at position x and time t effective diffusion coefficient (aquiclude) scaled for tortuosity longitudinal mechanical dispersion coefficient (aquifer) transverse mechanical dispersion coefficient (aquifer) thickness averaged solute concentration at input (aquifer) dimensionless parameter (P2 I/Pl L) aquiclude thickness aquifer thickness time, and time increment, respectively period of oscillation associated with the sinusoidal input groundwater pore velocity [(Darcy velocity)/ (effective porosity) ] 2~ times the frequency of oscillation associated with the sinusoidal input longitudinal distance from the recharge zone, and distance increment dimensionless distance parameter ( Dx l Vl 2 ) distance associated with the thickness of the aquifer, and distance increment effective porosity in aquifer and aquiclude, respectively dimensionless time (Dt/l 2 )
(mol 1-1 ) (mol 1-1 ) (mol 1-1 ) (mol 1-1 ) (m 2 yr. -1 ) (m 2 yr. -1 ) (m 2 yr. -1 ) (mol 1-1 ) (m) (m) (yr.) (yr.) (m yr. -1 ) (Ty'l".-I)
(m) (m)
t i m e scales, t h e c h o i c e o f a c o n s t a n t value f o r DL m u s t be v i e w e d w i t h c a u t i o n since t h e m a g n i t u d e o f t h e d i s p e r s i o n m a y increase w i t h o u t limit as t h e geological size o f t h e s y s t e m increases b e c a u s e o f t h e successively larger scale o f n o n - u n i f o r m i t y w h i c h m a y b e e n c o u n t e r e d ( S m i t h a n d S c h w a r t z ,
1980). T h e d i s t a n c e t r a v e r s e d b y a m o l e c u l e o f s o l u t e u n d e r t h e a c t i o n o f longit u d i n a l m e c h a n i c a l d i s p e r s i o n in elapsed t i m e t since g r o u n d w a t e r recharge, is ( D L t ) i n . P r e d i c t e d values o f DL are o r d e r s o f m a g n i t u d e g r e a t e r t h a n t h o s e a n t i c i p a t e d f r o m l a b o r a t o r y studies; nevertheless, l o n g i t u d i n a l m e c h a n ical d i s p e r s i o n e f f e c t s c a n b e i g n o r e d w h e n t h e length scale a l o n g t h e aquifer, o v e r w h i c h c o n c e n t r a t i o n d i f f e r e n c e s are t o b e investigated, is m u c h g r e a t e r t h a n t h e d i s p e r s i o n d i s t a n c e (D L t) 1/2 . F o r e x a m p l e , t h e length scale f o r a
134 cyclic input of period T is VT; longitudinal mechanical dispersion can then be ignored when VT >~(DLt) 1/2 . This condition is met when the period T is sufficiently large. For a 104-yr. climatic cycle, and groundwater residence times t in the range 104--106 yr., VT is one to two orders of magnitude greater than ( D L t) I/2 , based o n D L / V ~- 3 m and V = 3 m yr. -1 . Note that the upper value chosen for t (106 yr.) far exceeds typical groundwater residence times, so that there is never sufficient time available for longitudinal mechanical dispersion to significantly influence the large-scale paleoclimatic oscillations. A different dispersion mechanism which may be more influential at long time scales is discussed below. The corresponding (analytical) model approximates conditions f o u n d in m a n y large basins and can be readily applied to very old groundwater. In its simplest form, it ignores longitudinal mechanical dispersion effects (given by D L ) but postulates a molecular diffusion limited exchange of material between the flowing groundwater and the stationary water in the confining aquiclude, leading to a spreading of the solute in the mobile water and a consequent reduction in the concentration variability from that existing at recharge. A characteristic time for a molecule of solute to diffuse across an aquiclude of thickness I is 12/D where D is the effective diffusion coefficient of the solute therein. Solute diffusion coefficients in water are of the order of 10 -s cm 2 s-1 which gives a value for D of 0.01 m 2 yr. -1 after allowing for a tortuosity factor of 3. With this value o f D , 12/D increases from 104 yr. when I = 10 m to ~ 1 0 s yr. when l = 3 0 m , suggesting that dispersion, which results from the interaction of axial flow in the aquifer and transverse molecular diffusion in the aquiclude, can be important on time scales similar to the periods of the major climatic cycles. An analogous model has been used by Turner {1958) and Skopp and Warrick (1974) to represent a porous medium with an internally distributed stationary phase (e.g., water in the boundary region of pore spaces) for which l is very small. Other theoretical considerations involving stratified media have been made by Shamir and Harleman'(1967), Fried and Combarnous (1971), and A1-Niami and R u s h t o n (1979). In this paper, the model is used to interpret aspects of the variations of stable-isotope ratios, inert-gas concentrations and chloride levels in groundwater. An assessment is made of the factors involved in interpreting the spatial variability of chosen groundwater properties in terms of paleoclimatic oscillations at recharge.
MATHEMATICAL
ASPECTS
Problem definition Consider an idealized homogeneous isotropic semi-infinite aquifer of height L with an adjacent confining aquiclude of thickness l which is also
135
stationary
aquic~ude aquifer
L
mobite
PV
Fig. 1. A schematic representation of a double-layer aquifer--aquiclude system. Water in the aquifer is moving longitudinally with pore velocity V and water in the aquiclude is stationary.
homogeneous and isotropic. This layered system is represented schematically in Fig. 1. Attention is focused on the saturated zone and it is assumed that no water reaches the aquiclude from above; all water is presumed to enter the system at the recharge zone (x -- 0). Flow in the aquifer is taken to be one-dimensional with a constant pore velocity V [= (Darcy velocity)/(effective porosity)] while water contained in the aquiclude is stationary. This section is concerned with the mathematical aspects of longitudinal convection and dispersion in the aquifer, from the recharge zone, of a nonreacting solute which can diffuse into and o u t of the adjacent aquiclude. The effects of radioactive decay of the solute and of adsorption onto the solid matrix of the system are ignored here b u t they can be easily included in the model. Solute concentrations are assumed to be low so velocity components induced b y mass transfer between layers can be ignored. The convection--diffusion equations for the solute in terms of rectangular coordinates ( x , y ) and time t are:
0C 1/Ot -k V(OC 1/Ox)
-- D w (02 C 1 / 0 y 2 )
(1)
and
0C2/0t = D(0 2 c 2 / 0 y 2)
(2)
where axial mechanical dispersion in the aquifer and axial diffusion in the aquiclude have b o t h been ignored (i.e. terms involving 0 2 C1/Ox 2 and 0 2 C2 / 0x 2 are not included). Henceforth, the process governed by eqs. 1 and 2 will be called double-layer dispersion. Solute concentrations in aquifer and aquiclude water are given by C1 (x, y, t) and C2 (x, y, t), respectively. D w is the transverse dispersion coefficient in the aquifer and D is the molecular diffusion coefficient (scaled for tortuosity) in the aquiclude. Boundary conditions describing continuity of solute concentration and flux at the interface y = L are: C1 (x, L, t) = C2 (x, L, t)
(3a)
and 0C1 8C2 plDw -~y ( x , L , t ) = IJ2D- : - ( x , L , t ) Oy
(3b)
136
where Pl and #2 are the porosities of the aquifer and the aquiclude, respectively. Zero-flux conditions ~C1/Oy = 0 and ~C2/~y = 0 are chosen at y = 0 and y = L + l, respectively. Taking the average of eq. 1 over the thickness L of the aquifer gives: DC~ + V -~C1 - - - •2 D -[~C2]/
~t
~x
(4)
Pl L \ ~y ly=L
after applying the boundary conditions at y = 0 and y = L, and defining L
Ca ~- L
C1 dy 0
The term on the right of eq. 4 represents the flux of solute passing between aquifer and aquiclude. Since diffusion in the aquiclude is m u c h slower than transverse dispersion in the aquifer (i.e. D ~ DT ), transverse gradients (OC1/Oy) in the aquifer will be m u c h smaller than those in the aquielude provided that the transverse concentration distribution in the aquifer is initially uniform or that sufficient time has elapsed for the ~C1/Oy to decay (on a time scale corresponding to diffusion in the aquiclude, transverse dispersion in the aquifer is instantaneous). Under these circumstances, C1 may be replaced by C1 in condition (3a) and the problem becomes identical to that described by Skopp and Warrick (1974) for which the Laplace transform solution has t h e form: Lap(C1 ) = Lap(F) exp
( sx V
x #2 (sD)lntanhl(s/D)ln) VL Pl
(5)
when Cj (x, 0) = C2 (x,y, 0) = 0, Ca (0, t) = F(t) and Lap denotes a Laplace transform (e.g., Lap(Ca ) = f g e-"tCl dt).
Step-function input (an unsteady problem) S k o p p and Warrick (1974) inverted eq. 5 for a step-function input [F(t) = 1] to give the solution: C1
2 of 7r
x
[ ~(sinhbX--sinbX)] exp -\cosh bX + cos bX
sin
t--
2
2 \e-osh ~ + -e-oosb-X])
7' and
fort>x/V
(6)
137
C1 = O, f o r t ~ x / V where
a = (p2Dinx)/(pl VL)
and
b =
liD 1/2 when Dt/l:
An independent analysis shows that ~ 1 (i.e. the elapsed time is much smaller than the characteristic transverse diffusion time for the aquiclude), this solution m a y be approximated b y that corresponding to an aquiclude of infinite thickness, viz.:
C1 = ,I o,erfc[½a/(t--x/V)U2]'
forf°r t~x/>x/V yt
(7)
where erfc is the c o m p l e m e n t a r y error function. In principle, it is n o w possible to determine the response to an arbitrary input function b y means of a convolution integral. However, in practice it is more convenient to use a numerical method. The details of this technique are given later (p. 10). Fig. 2 represents a comparison between double-layer (solid curves) and longitudinal mechanical (dashed curves) dispersion in the case of a step input. The solid curves are determined numerically for convenience (after verifying that agreement with the results of eqs. 6 and 7 is better than 1%) and the dashed curves are calculated from eq. 10.6.22 in Bear (1972) for
I.--0
0.1
I
05
8
I-0
~___ I-5
2"0
x=£_~ V.~2
c", LLI N N
\ \ \
0
6'0
I
I
7"0
8"0 Dx x=~-2
9'0
I0'0
Fig. 2. R e s p o n s e o f t h e m o d e l a q u i f e r t o a s t e p - f u n c t i o n i n p u t o f s o l u t e c o n c e n t r a t i o n . The solid curves illustrate the dispersion due solely to a diffusion limited exchange of s o l u t e b e t w e e n f l o w i n g a q u i f e r w a t e r a n d i m m o b i l e a q u i c l u d e w a t e r w h e n K = P2 l / P l L = 0.25. The dashed curves represent only macroscopic dispersion due to the non-uniformities in h y d r a u l i c c o n d u c t i v i t i e s f o r V---- 3 m yr. -1 , D L = 10 m 2 y r . -1 a n d are p l o t t e d a g a i n s t d i m e n s i o n l e s s v a r i a b l e s X = D x / V l 2 a n d T = D t / l 2 a s s u m i n g l = 10 m , D = 0.01 m 2 yr. -1 "
138 V = 3 m yr. -1 and D L = 10 m 2 yr. -1 and plotted against dimensionless distance (X = D x / V l 2 ) and dimensionless time (r = D t / l 2 ) in Fig. 2, assuming t h a t l = 10 m and D = 0.01 m E yr. -1 . Increasing l has the effect of steepening the dashed curves. We first note that double-layer dispersion becomes significant when T is of order one (i.e. t of order 12/D) for K [= (P21/PlL)] = 0.25 and t h a t this critical value of r varies inversely with the exchange capacity coefficient K (as was expected). The next important observation is that it is indeed appropriate to ignore longitudinal mechanical dispersion (regardless of a priori arguments about its importance or otherwise at long time scales) when it is compared with double-layer dispersion for K = 0.25 and r ~> 1, since the dashed curves become much steeper than the solid curves in Fig. 2. The choice of K = 0.25 in this model is based on the following parameter values: P2 = 0.05, Pl -- 0.20 and IlL = 1 (i.e. identical thickness for both aquifer and aquiclude). The choice of P2 assumes an aquiclude such as shale. Alternatively, for a clay aquiclude we might choose P2 = 0.5, Pl = 0.2 and IlL = 0.1 to again obtain K = 0.25.
Sinusoidal input (a steady-state problem) Turner (1958) considered the frequency response of a packed-bed model which included stagnant pockets. He obtained the solution for the response to a periodic input, after a steady cyclic state had been reached, directly from eqs. 2 and 4. When C1 (0, t) = A 0 sin Wt where W = 2~r/T and T is the period of the oscillation, Turner's solution becomes:
C1 = Ao e x p ( - - p x / V ) s i n ( W t - - q x / V )
(8)
for the case w i t h o u t longitudinal mechanical dispersion, where:
KO (sinh bO -- sin bO t P = 2--b-\cosh bO + cos bO] KO ( sink bO + sin bO q = W+-2b \ e o s h bO + cos bO !
and
0 = (2W) in
This result can also be obtained quite readily from eq. 5 by using the expansion theorem given by earslaw and Jaeger (1953, p. 280) for large time solutions. If A(x) is the amplitude of the oscillation as it proceeds along the aquifer, the amplitude attenuation factor is given by:
A/Ao
= exp(--px/V)
and the wave speed by WV/q. The attenuation factor is plotted against longitudinal distance in Fig. 3. Note that, as was expected, the higher the frequency of oscillation the more rapidly it damps out. In Fig. 4, the decay constant p is plotted against aquiclude thickness I for
139
¢: o
-5 IE
5
10
15
20
Dx X = Vt 2
Fig. 3. Relative amplitude attenuation with dimensionless distance X = D x / $ r l 2 of a sinefunction input concentration (period T) after a steady cyclic state has been reached throughout the model aquifer, where K = g2 l/P1L = 0.25 and T* = DT/I 2.
% × (3_
©
I
I
I
J
20
40
60
80
100
(m)
Fig. 4. Amplitude decay constant p for a periodic input concentration (period T) is plotted against aquiclude thickness I for fixed K = tl 2 I / p 1 L = 0.25.
t h e t h r e e d o m i n a n t climatic oscillations w h e n K is f i x e d at 0.25 a n d D at 0.01 m 2 yr. -1 . F o r a given p o r o s i t y ratio this c o r r e s p o n d s t o h o l d i n g IlL fixed. Each curve e x h i b i t s an o p t i m u m value o f I w h i c h achieves a m a x i m u m a m p l i t u d e d e c a y rate. These o p t i m u m values are 10, 14 and 20 m for t h e periods 2 3 , 0 0 0 , 4 2 , 0 0 0 a n d 1 0 0 , 0 0 0 yr., respectively.
The numerical m e t h o d A c o n v e n i e n t and rapid m e t h o d o f solving eqs. 2 and 4 for an a r b i t r a r y i n p u t f u n c t i o n is t h e n u m e r i c a l t e c h n i q u e d e s c r i b e d below. (1) D e f i n e n o d a l p o i n t s (i = 1 , . . . , I ) o f c o n s t a n t spacing Ax along the length o f t h e o n e - d i m e n s i o n a l m o d e l a q u i f e r and, at each such p o i n t , d e f i n e
140
nodes (j = 1, . . . , J ) of constant spacing Ay along the transverse (y-) direction. The objective is to solve eqs. 2 and 4 by replacing them with difference equations at each nodal point and then marching through time t in discrete steps (n = 1 , . . . , N) of length At. Let (C 1 )[' be the value of C1 at position x = izkx and (C2)i~,j be the value of C2 at x -- iAx and y -- L = j A y , all evaluated at time t = n a t . (2) Let A x = V A t and move the concentration distribution C1 forward a distance Ax by setting (C1)* = (C1)['-1 • This is equivalent to solving:
0C 1/0t "~ V(0C 1/Ox)
= 0
(9)
(3) Now solve eq. 2 by replacing it with the following difference equation:
[(C2 ~.n+l ,,,~ - - (C2)i,i]/At "
= D I { ( C 2 ) ¢ ,~+1 , s + l - - 2 ( C 2 ) i , 1 .+1
-[- (C2)n,~ll}/(Ay)2] (10)
where (C2 ~,~1 = (C1)* and (C2 )t,J ~,+1 = (C2 Ji,J-I ~,÷1 are the corresponding difference equations for the b o u n d a r y conditions. This equation is called "implicit" because the space derivative is taken at the new time (n + l ) A t rather than at n A t and has the property of numerical stability (i.e. errors are not amplified during the calculation). Its solution is a standard procedure (see, e.g., Richtmyer and Morton, 1967). (4) Finally derive (C1)7 +1 from the difference approximation to aC1
P2 D [ a C 2 ~
(10)
viz.
( e l ) n• + l - (el)~ At
]../2
D ((C2'n+l--(c2)n+l I) Ji,2
Pl L
(11)
Ay
The above procedure involves splitting eq. 4 into two parts, the first being pure convection and the second the exchange of solute between strata. The choice of Ax = V A t ensures t h a t the solution of eq. 9 at nodal points is exact so no artificial numerical dispersion is introduced. As mentioned earlier, results have been checked against the analytical solutions ( 6 ) a n d (7) for a step-function input. In this case, ~ 1% accuracy or better was obtained by choosing Ay = 0.025l and At at least 0.01 times the elapsed time t for D t / l 2 ~ 10. However, for an arbitrary input function, the errors incurred may be different so the sensitivity of the results to At and Ay should always be checked. DISCUSSION General c o m m e n t s
The model described above is applicable to m a n y large hydraulic basins in which well-defined aquifer layers are bounded by low-permeability aqui-
141
cludes. It has the following desirable features: (1) It is the first example of a model evaluated at long time scales which predicts attenuation in the cycle of groundwater solute concentrations which oscillate with periods exceeding 104 yr. The predicted dispersion is a consequence of the differences in transport properties within the aquifer and the confining aquiclude. The choice of K = 0.25 in the model is discussed in the Section " S t e p function input (an unsteady problem)". (2) Analytical solutions of some simple input functions can be obtained. The unsteady response to a step input may be considered as a model of the response to non-steady deviations of solute concentrations from the underlying periodic trend which may, in turn, be modelled by the steady-state cyclic response to a superposition of sine and cosine inputs. (3) The model can readily be extended to study the effects of dispersion of generalized input functions using a simple numerical solution technique. (4) Subject to symmetry restrictions, the analytical model can be extended to include multiple layers of aquifer and aquiclude. These restrictions can be removed when using the numerical model. (5) Without conceptual change, the model can be extended to account for adsorption onto the host material of either the aquifer or aquiclude, and the effects of radioactive decay. (6) It is possible that the model might predict a significant isotopic difference when the isotope effects associated with adsorption and molecular diffusion in the aquiclude are taken into account.
Applications Dispersion of dissolved noble gases in groundwater. Variations in the concentrations of noble gases in groundwater have recently been used to monitor paleotemperature changes in recharge areas. Rainfall has been shown to saturate with inert gases before infiltration (Mazor, 1972). However, because re-equilibration would normally occur during the lengthy recharge process, the temperature monitored is probably that in the vicinity of the water table. Inert-gas concentrations will only reflect climatic changes over long periods in the absence of significant losses in the saturated zones. Some experimental evidence for this assumption has been discussed by Mazor (1972) in his study of the Jordan Rift Valley, Israel. A groundwater piston-flow model is tacitly assumed in the normal approach, where gaseous concentrations are related to recharge temperatures through laboratory solubility data. The current model illustrates that the simplification of ignoring dispersion effects may not be justified, depending on the aquifer and aquiclude thicknesses. In the Section "Mathematical aspects", we indicated that when K = 0.25, the double-layer dispersion of a step input becomes significant after an elapsed time of 12/D and that the relative amplitude attenuation of the dominant climatic oscillatory modes may be calculated using Fig. 4. The value of D was taken to be 0.01 m 2 yr. -1 .
142 A study has been made of the Ar and Kr levels in the Bunter Sandstone groundwater in the English Midlands (Andrews and Lee, 1979; Bath et al., 1979). The schematic cross-section (Bath et al., 1979, fig. 1) indicates that l and L are "~200 m. In this circumstance, 12/D is 4" 106 yr. and the amplitude of the climatic modes is reduced by only 2--4% at an equivalent distance of 3" 104 yr. from the recharge area. Clearly, the Ar and Kr levels predicted from the dispersion model will be little changed over this distance owing to the large layer thickness chosen. Had the sequence comprised an assembly of layers with thicknesses of order 10 m, the effect of dispersion would have been significant. For a double layer with l = L = 10 m, the amplitude of a 23,000-yr. climatic cycle for example is reduced by 57% at a point 30,000 yr. from the recharge area and the critical time for significant dispersion of a step input is of the order 104 yr. (12/D). Extension of the model to multiple layers is possible as indicated above. There has been a t e n d e n c y in more recent times to recommend monitoring the Kr/Ne or Ar/Ne ratios rather than the absolute concentrations (Mazor, 1972; Phillips and Davis, 1979}. The solubilities of the lower-atomicweight inert gases are far less temperature-dependent. Experimentally the advantages lie in the fact that inert-gas ratios are less sensitive than absolute concentrations to changes accompanying the evolution of soil gas from air, and to any degassing that might occur in the saturated zone. It is therefore important to note t h a t the model can only be applied directly to absolute solute concentrations, because the predicted dispersion depends on the molecular diffusion coefficient in the aquiclude and therefore on the atomic weight of the gas. Different gases disperse to different extents, and simple expressions cannot be developed for their ratios.
Stable-isotope ratio variation in the Australian Great Artesian Basin. A study has recently been reported of aspects of the isotope hydrology of the principal Jurassic aquifer of the Australian Great Artesian Basin (Airey et al., 1979; Hartley, 1979). The " a g e " of the groundwater, which was calculated from the 1900 potentiometric surface and the horizontal hydraulic conductivities increased regularly with distance from the recharge area from less than 2" 104 to more than 3 • 10 s yr. During this portion of the late Quaternary there were two major ice ages and an interglacial period. Altogether, 69 wells were sampled along defined flow lines in the principal study area. Despite the climatic changes the D/H ratios were f o u n d to be remarkably constant. The mean value was --41.8°/0o, and the standard deviation of the population 1.3. Part of the variability is u n d o u b t e d l y due to experimental error which is conservatively assessed to be 0.7~/00 . The 51SO-values were 6.6 + 0.9~/00 . Because of the relatively large experimental error (~ 0.5°/0o ), 40 of the samples were re-run on equipment in which the assessed errors at the l a level varied from 0.03 to 0.080/00 . The average value of this set was - - 6 . 8 7 + 0 . 7 7 % . Since the 5D- and 51SO-values are consistent with the "meteoric line", 5D = 85180 + 10, it would be expected that A(SD)/A{51So) -
-
143 would be of the order of 8. Experimentally the variations in 5D and 6180 were comparable. If confirmed, this result would be very significant and will u n d o u b t e d l y help to distinguish between various transport models. Confirmation is required because the samples were stored (albeit in tightly stoppered glass bottles) for over a year between collection and the more precise measurements. This discussion will therefore be concerned principally with the constancy of the 5D-values. It resolves itself into two parts: (1) the input function; and (2) dispersion within the saturated zone. (1) Input function: There is a good linear correlation on the global scale b e t w e e n the stable-isotopic composition of rainfall samples and the mean annual surface temperature. Dansgaard (1964) has shown that: A(SD)/A(temperature) ~ 5.6
and
A(5 lSo)/A(temperature) ~ 0.69
Care must be taken when assuming that variations in the isotopic composition of groundwater simply reflect temperature changes in the recharge zone. In fact they respond to a range of global and local parameters. Hartley (1981) has gone so far as to suggest that: "there seems little possibility of relating climate to the isotopic composition of groundwater other than by empirical relationships developed from other paleoclimatic evidence" Nevertheless, a simple explanation based on temperature changes and the Dansgaard correlation proved adequate in a number of instances. Examples include recent studies of the isotope hydrology of a number of aquifers in Saudi Arabia (Shampine et al., 1979) and of thermal waters in the Budapest region, Hungary (Deak, 1979). On the other hand, there are many cases in which the paleotemperatures calculated from stable-isotope variations are substantially less than those assessed from other evidence. For instance, in the Bunter Sandstone study (Bath et al., 1979), the 5 lSO-values may be interpreted as indicating a change in cloud formation temperature of at least 1--2°C over the time scales involved (35,000 yr.). Surface temperature changes of from 5 to 6°C were estimated from gas solubilities and other evidence. Again the range in 5180values from the other British aquifer systems discussed b y Evans et al. (1979) is similar. From an assessment of the stable-isotope variability at a number of sites it can be reasonably expected that, over the time scales accessible to 14C, the 51So range would be at least 2~/00, and the 5D-values between a b o u t four and eight times this value. Thus, the range in 5D-values found in the Great Artesian Basin was a b o u t an order of magnitude less than expected. Before discussing the effects of dispersion of 5D variability, t w o possible explanations associated with the nature of the input function need to be considered. As was pointed o u t by Evans et al. (1979), the coefficient of variation of 5D with temperature for m o n t h l y rainfall samples in a particular region may
144 be much less than t h a t predicted from Dansgaard's global correlation. Such arguments cannot be directly used in this study because the 8D-values of the rainfall from the two available stations in the general recharge area (Charleville and Longreach) exhibit partial correlation coefficients with rainfall a m o u n t which are substantially greater than those with temperature (Hartley, 1979, 1981). It is, in principle, possible that various climatic factors could act in opposite directions to effectively cancel the expected 5D variability. This argum e n t is countered by the observation that a group of fourteen wells in South Australia which are recharged from the western extremities of the Great Artesian Basin are also very constant. They exhibit a mean of--48.1°/00 and a standard deviation of 1.3°/00 (Airey et al., 1979; Hartley, 1979). The probability that the combination of climatic factors in two recharge areas separated by 1200 km would affect the variation of D/H ratios in the groundwater in the same way must be very small. The explanation for the observed constancy in the 5D-values must be associated with dispersion over the long residence times. (2) Dispersion effects. The Great Artesian Basin (Habermehl, 1980) extends to distances equivalent to ~ 3 . 1 0 s yr. from the recharge area; it is therefore appropriate to consider the 10 s -yr. climatic cycle (Hays et al., 1976). The amplitude of this cycle will attenuate more slowly than the other more rapid climatic oscillations and is considered in relation to the observed reduction in D/H variability. When l = L = 100 m, for example, calculations based on Fig. 4 indicate an amplitude reduction, for the l 0 s -yr. cycle, of 25% at an equivalent distance of 2" l 0 s yr. from the recharge area. This reduction increases to 46% when l = L = 10 m is chosen. Apparently the model cannot fully explain the gross smoothing of the 5D-values in the Great Artesian Basin. Further work is required, because it appears that the attenuation of stable-isotope ratio variability may provide a critical test of solute dispersion in very large hydraulic basins.
Groundwater chloride. Chloride concentration at recharge depends on a number of factors including: (1) The rate of input in the recharge area by chloride in rainfall and dry fall out; (2) The rate of input from the dissolution of aquifer rocks and sediments; (3) The rate of loss due to surface run off; (4) The net groundwater discharge. Because of the complexity of the overall process, it is not possible to define the principal climate determining parameter. Nevertheless, it has been shown in studies of the Jurassic sandstone aquifer of the Australian Great Artesian Basin (Airey et al., 1979) and in the Bunter Sandstone aquifer, England (Bath et al., 1979) that the plots of chloride concentration with groundwater age fall on well-defined curves with a qualitative similarity. In both
145 cases the levels decrease regularly with age from relatively high levels characterizing the recharge water. In the Great Artesian Basin study there is evidence for a secondary m a x i m u m during the last interglacial, ~ 1 . 1 . l 0 s yr. ago. The current model can be used to assess the extent of chloride ion dispersion over this time scale and, in principle, to provide a better measure of the input function by a systematic trial-and-error procedure, using the numerical model. This could be important when studying the variations in the chloride balance in the recharge area t h r o u g h o u t the late Quaternary. Information on the rainfall variation over the past 1.2" l 0 s yr. is already available (Kershaw, 1978) from a station in north Queensland. At this stage no a t t e m p t has been made to establish these balances.
CONCLUSIONS (1) The dispersion process described here can significantly reduce the variability of solute concentrations within an aquifer to below that present in recharging water over time scales which depend on the aquifer and aquiclude thickness. When K - - 0 . 2 5 , the dispersion of a step input becomes significant after an elapsed time of 12/D (e.g., 104 yr. when l = 10 m and D = 0.01 m 2 yr. -I ) and the relative amplitude attenuation of d o m i n a n t steady-state climatic oscillations m a y be calculated using Fig. 4. Such effects need to be taken into account when interpreting noble-gas concentrations, stable-isotope ratios and chloride levels in terms of paleoclimatic parameters. (2) The extent of dispersion predicted by the present model is not sufficient to explain the observed reduction in variability of stable-isotope ratios in the very old groundwater of the Great Artesian Basin. Further work is required as this observation m a y form the basis of a m e t h o d for the critical comparison of transport models in very large aquifer systems. (3) The numerical solution described here is both fast and reliable enabling the easy determination of the response of an aquifer--aquiclude system to an arbitrary input function from an arbitrary initial state. The technique m a y be used to extend the model to multiple aquifer--aquiclude layers.
ACKNOWLEDGEMENTS The modelling ideas in this work owe much to some stimulating discussions with Peter Hartley of the A.A.E.C. The authors are also indebted to Professor A.T. Wilson of the University of Waikato, New Zea 1~ .d, for the high precision 1so--160 mass spectrometry, and to Mr. D. R o m a , (A.A.E.C.) for preparing the Great Artesian Basin samples.
146 REFERENCES Airey, P.L., Calf, G.E., Campbell, B.L., Habermehl, M.A., Hartley, P.E. and Roman, D., 1979. Aspects of the isotope hydrology of the Great Artesian Basin, Australia. In: Proc. Symp. on Isotope Hydrology, 1978, Vol. I, Neuherberg, Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 205--219. Airey, P.L., Calf, G.E., Hartley, P.E. and Roman, D., 1980. Aspects of the isotope hydrology of two sandstone aquifers in arid Australia. In: I.A.E.A. Adv. Gr. on Application of Isotope Techniques to Arid Zones Hydrology, Vienna, Nov. 1978. Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 93--111. A1-Niami, A.N.S. and Rushton, K.R., 1979. Dispersion in stratified porous media: analytical solutions. Water Resour. Res., 15: 1044--1048. Andrews, J.N. and Lee, D.J., 1979. Inert gases in groundwater from the Bunter Sandstone of England as indicators of age and paleoclimatic trends. J. Hydrol., 41: 233--252. Bath, A.H., Edmunds, W.M. and Andrews, J.N., 1979. Paleoclimatic trends deduced from the hydro-geochemistry of a triassic sandstone aquifer, U.K. In: Proc. Symp. on Isotope Hydrology, 1978, Vol. II, Neuherberg, Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 545--568. Bear, J., 1972. Dynamics of Fluids in Porous Media. American Elsevier, New York, N.Y., 764 pp. Carslaw, H.S. and Jaeger, J.C., 1953. Operational Methods in Applied Mathematics. Oxford University Press, London, 2nd ed., 359 pp. Dansgaard, W., 1964. Stable isotopes in precipitation. Tellus, 16: 436--468. Dansgaard, W., Johnsen, S.J., M~bller, J. and Langway Jr., C.C., 1969. One thousand centuries of climatic record from Cape Century on the Greenland ice sheet. Science, 166: 377--381. Deak, J., 1979. Radiocarbon dating of thermal waters in the Budapest area. In: Int. Congr. on Isotopes in Nature, Leipzig, Pap. AUS-8-1/80. Delmas, R.T. and Ascenia, J.M., 1980. Polar ice evidence that atmospheric CO 2 20,000 year B.P. was 50 per cent of present. Nature (London), 284: 155--156. Evans, G.V., Otlet, R.L., Downing, R.A., Monkhouse, R.A. and Rae, G., 1979. Some problems in the interpretation of isotope measurements in United Kingdom aquifers. In: Proc. Symp. on Isotope Hydrology, 1978, Vol. II, Neuherberg, Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 679--708. Fried, J.J. and Combarnous, M.A., 1971. Dispersion in porous media. Adv. Hydrosci., 7 : 169--282. Galbally, I.E. and Freney, J.R., 1978. Modification of climate. In: A.B. Pittock, L.A. Frakes, D. Jenssen, J.A. Petersen and Z.W. Zillman (Editors), Climate Change and Variability - - A Southern Perspective. Cambridge University Press, Cambridge, pp. 269--282. Habermehl, M.A., 1980. The Great Artesian Basin, Australia. Aust. Bur. Miner. Resour., J. Geol. Geophys., 5: 9--38. Hartley, P.E., 1979. Deuterium/hydrogen ratios in Australian rain and groundwater. M.Sc. Thesis, University of New South Wales, Kensington, N.S.W. Hartley, P.E., 1981. Deuterium/hydrogen ratios in Australian rainfall. J. Hydrol., 50: 217--229. Hays, J.D., Imbrie, J. and Shakleton, N.J., 1976. Variations in the earth's orbit; pacemaker of the ice ages. Science, 194: 1121--1132. Kershaw, A.P., 1978. Record of last interglacial--glacial cycle from northeastern Queensland. Nature (London), 272: 159. Lenda, A. and Zuber, A., 1970. Tracer dispersion in groundwater experiments. In: Proc. Symp. on Isotope Hydrology, 1970. Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 619--637.
147 Mazor, E., 1972. Paleotemperatures and other hydrological parameters deduced from noble gases dissolved in groundwaters, Jordan Rift Valley, Israel. Geochim. Cosmochim. Acta, 36: 1321--1336. Phillips, F.M. and Davis, S.N., 1979. Noble gases in groundwater as paleoclimatic indicators. EOS (Trans. Am. Geophys. Union), 6 0 : 8 3 0 (abstract). Richtmyer, R.D. and Morton, K.W., 1967. Difference Methods for Initial-Value Problems. Interscience, New York, N.Y., 2nd ed., 405 pp. Shamir, U.Y. and Harleman, D.R.F., 1967. Dispersion in layered porous media. J. Hydraul. Div., Proc. Am. Soc. Civ. Eng., 93(HY5): 237--260. Shampine, J.W., Dinner, T. and Noory, M., 1979. An evaluation of isotope concentrations in the groundwater of Saudi Arabia. In: Proc. Symp. on Isotope Hydrology, 1978, Vol. II, Neuherberg, Int. At. Energy Agency (I.A.E.A.), Vienna, pp. 443--463. Skopp, J. and Warrick, A.W., 1974. A two-phase model for the miscible displacement of reactive solutes in soils. Soil Sci. Soc. Am. Proc., 38: 525--550. Smith, L. and Schwartz, F.W., 1980. Mass transport, I. A stochastic analysis of macroscopic dispersion. Water Resour. Res., 16: 303--313. Turner, G.A., 1958. The flow-structure in packed beds -- a theoretical investigation utilizing frequency response. Chem. Eng. Sci., 7: 156--165.