- Email: [email protected]
Modelling of helium transport in groundwater along a section in the Pannonian basin
Journal of Hydrology 225 (1999) 185–195 www.elsevier.com/locate/jhydrol
Modelling of helium transport in groundwater along a section in the Pannonian basin L. Cserepes, L. Lenkey* Department of Geophysics, Lora´nd Eo¨tvo¨s University, Ludovika te´r 2, 1083 Budapest, Hungary Received 14 April 1998; accepted 20 September 1999
Abstract Underground water flow in sedimentary basins controls the distribution of dissolved salts and gases, and their concentrations may therefore be used as indicators of the flow direction. Recent measurements of the 4He concentration in deep waters of the Pannonian basin have great importance in this respect. This paper presents an example of the simultaneous computation of water flow and helium distribution along a section crossing the Great Hungarian Plain. The model consists of three permeable layers. The boundaries of the layers are prescribed using geologic sections constrained by ample borehole and seismic data. The results of the finite-difference calculations are fitted to the observed helium concentrations using a least-squares algorithm that varies the model parameters. The significance of the model is that it reconstructs the structure and flux of the groundwater flow and estimates the poorly known hydrogeological parameters of the flow regime such as hydraulic conductivities, conductivity anisotropy and dispersion coefficients. The statistical uncertainty of the estimated parameters is around half an order of magnitude. An estimate of the regional average of the incoming helium flux is also obtained. The total helium flux in the Great Hungarian Plain at the surface is within the range observed in old stable continental areas of the Earth. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Groundwater flow modelling; Helium transport; Hydraulic parameters; Least-squares estimation
1. Introduction Groundwater flow in sedimentary basins controls the distribution of dissolved salts and gases. Among them, helium has received increasing attention as a tracer of groundwater movement (Andrews and Lee, 1979; Torgersen and Clarke, 1985; Torgersen and Ivey, 1985; Stute et al., 1992; Zhao et al., 1998). Helium is an inert gas, therefore its underground distribution is governed by relatively simple laws. During recharge, the helium contained in soil air is dissolved at the water table. Along the flow path * Corresponding author. E-mail address: [email protected] (L. Lenkey)
toward the discharge zone, the concentration of helium increases due to (1) the accumulation of helium produced in situ in the aquifer by radioactive decay of U and Th, and (2) by helium entering the aquifer from greater depth. Measuring the helium concentrations, we can obtain information about the direction and flux of the groundwater flow, age of the groundwater, and hydrogeological parameters of the aquifer. Most of the Earth’s helium is of radiogenic origin. The two main sources, U and Th, which produce 4He by a-particle decay, are highly concentrated in the continental crust. A minor part of the helium is primordial and escapes from the mantle at midocean ridges and through the continental lithosphere
0022-1694/99/$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(99)00158-4
186
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
in tectonically active areas, particularly those undergoing extension or volcanism (O’Nions and Oxburgh, 1988). The isotopic composition of helium is considerably different in the mantle, crust and atmosphere. This allows distinction among the radiogenic crustal, mantle-derived, and atmosphere-derived helium by measuring the ratio of helium isotopes ( 3He/ 4He). More than 200 helium concentrations and isotope ratios have been measured in groundwater samples from Hungary (Martel et al., 1989; Stute et al., 1992; see Fig. 1). The results show that 1–16% of the helium derives from the mantle. The presence of mantle-derived helium in groundwaters of the Pannonian basin is in accordance with the Neogene extensional origin of the basin (Royden and Horva´th, 1988) and with geophysical observations showing elevated upper mantle beneath the basin (Horva´th, 1993). All measured groundwater samples from the Pannonian basin contain 4He in excess of the equilibrium saturation with air. The concentration of dissolved helium varies three orders of magnitude, increasing from the recharge areas toward the discharge areas ( 1–1300 × 1012 4He atoms cm 23 H2O). The samples with the highest concentration come from the central part of the discharge zone. Application of simple one- and two-dimensional groundwater flow models for helium transport allows determination of the groundwater age, and quantification of the helium flux coming from the deeper crust and mantle (Martel et al., 1989; Stute et al., 1992). The goal of this study is to show that by calculating the 4He concentrations from a two-dimensional transport model, and fitting the results to the measured data, it is possible to reconstruct groundwater flow and He migration as well as to estimate the value of poorly known hydrogeological parameters such as
hydraulic conductivity and anisotropy. The groundwater flow and the helium concentrations are modelled along a section crossing the Great Hungarian Plain (GHP) (Fig. 1).
2. The study site The Pannonian basin is located in eastern Central Europe. It is a topographically low region which is about 400 km from N to S and 600 km from W to E and is encircled by the Eastern Alps to the west, the Western Carpathians to the north, the Apuseni Mountains to the east and the Dinarides to the south (Fig. 1). It is not a single basin, but a system of deep subbasins separated by shallow basement highs and low hills. The two largest subbasins are the GHP and the Little Hungarian Plain, where the thickness of NeogeneQuaternary sediments reaches 7–8 km in the deepest troughs. The Neogene-Quaternary sedimentary sequence in the GHP consists of alternating layers of coarsegrained and fine-grained sediments. The basin was filled by a large delta system during the Late Miocene. Based on the lithology, the delta sediments are divided into Lower and Upper Pannonian rocks. The Lower Pannonian sequence comprises clay, marl and shale-rich sediments containing sandstone bodies upward. The Upper Pannonian sequence consists of alternating clay and sandstone. The Pliocene is made of lacustrine sands and clays, the Quaternary consists of poorly consolidated fluvial sand and clayey sand deposits. The stratigraphy of the Neogene-Quaternary sedimentary sequence in the GHP is summarized in Table 1 (after Urbancsek, 1977; Juha´sz, 1991, 1994). There are two regional groundwater flow systems in
Table 1 Stratigraphy of the Neogene-Quaternary sediments in the Great Hungarian Plain (after Urbancsek, 1977; Juha´sz, 1991, 1994) Age
Facies
Lithology
Thickness (m)
Quaternary Pliocene Upper-Miocene Upper Pannonian Lower Pannonian Middle-Miocene
Terrestrial, fluvial Fluvial, lacustrine
Sand and clay Sand and clay
100–600 100–1800
Lacustrine, delta plain Delta-slope, prodelta Prodelta, deep basin Basal
Sandstone, clay, marl Clay, marl, sandstone bodies Laminated clayey marl Conglomerate, clay
800–2000 600–3000 20–800 ,100
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195 Fig. 1. Topography of the Pannonian basin and its surroundings. Thick solid line shows the location of the hydrogeological cross section shown in Fig. 2: (1) 4He concentration measurement; (×) 3He/ 4He measurement.
187
188
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
Fig. 2. Generalized hydrodynamic system of the Great Hungarian Plain after Erde´lyi (1979). For location of the section see Fig. 1. RE, Recharge areas, D, Discharge areas, Q, Quaternary, Pl, Pliocene, LP, Lower Pannonian, UP, Upper Pannonian, B, Basement. Note that the depth scale is logarithmic enhancing the flow system in the Quaternary aquifer.
the GHP. The upper flow system is in the Quaternary sediments. These unconsolidated, thick and interconnected sand layers have good conductivity allowing the circulation of large amounts of groundwater of meteoric origin. Water recharge occurs at topographic highs and discharge occurs in topographic lows (Fig. 2) (Erde´lyi, 1979). The aquifer provides water for industrial and drinking purposes. The lower flow system exists in the Upper Pannonian sediments. The Upper Pannonian sand layers constitute a multilayered aquifer. The water is mostly of meteoric origin, mixed with connate water migrating upwards from deeper layers due to compaction. Due to the high geothermal gradient of the GHP (40–508C km 21), the Upper Pannonian aquifer contains thermal water which is used for balneotherapy and heating. The Pliocene sediments form a semipervious aquitard between the Quaternary and Upper Pannonian flow systems. The Lower Pannonian sequence can be regarded as an aquiclude, because the sandstone bodies are separated by thick marl layers and the communication between the high conductivity rocks is very poor. Our model is restricted to a 168 km long section that connects a main recharge and discharge area: a topographic high in the north-eastern part of the GHP (recharge) and a topographic low in the central part of the plain (discharge). The section follows a principal flow direction identified previously by pressure data (Erde´lyi, 1979). The helium concentration data collected in the close neighbourhood of the section are listed in Table 2. All the samples were taken from the
Quaternary aquifer, in the depth range of 71–405 m, most of them between 110 and 250 m depths.
3. Model of helium transport The hydrological model (Fig. 3) is built of three permeable layers, Quaternary, Pliocene and Upper Pannonian. The Lower Pannonian unit is regarded impermeable (Erde´lyi, 1976; Stute and Dea´k, 1989; Stute et al., 1992). The boundaries of the layers are prescribed using geologic sections constrained by ample borehole and seismic data (Urbancsek, 1977). We assumed the layers to be homogeneous. Each permeable layer consists of alternating sand and clay layers. We account for this effect by assigning anisotropy to the layers. The numerical calculations are carried out in a rectangular model domain. Its bottom is prescribed at a depth d 1400 m measured from the mean level of ground surface. This depth is within the impermeable Lower Pannonian unit. The basic law governing the large scale flow of water in porous-permeable media is Darcy’s law. We assume that the density of water is constant and the rock matrix is incompressible. Then a two-dimensional flow can be derived from a scalar stream function, c (e.g. Bear and Verruijt, 1987, p. 117). The horizontal (u) and vertical (w) Darcy velocities (fluxes) are given as u2
2c ; 2z
w
2c 2x
1
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
189
Depth b.s.l. (km)
0 Quaternary K 1xx
- 200 - 400
Pliocene K 2xx
- 600
Upper Pannonian K 3xx
- 800 - 1000 - 1200
Lower Pannonian
- 1400 0
50
100
150
Distance (km) Fig. 3. Geometry of the hydrogeological units as used in the model. Boundaries of the layers are taken from hydrogeological maps (Urbancsek, 1977) and well data. K1xx, K2xx and K3xx denote the different horizontal hydraulic conductivities of the layers. Vertical exaggeration is 40-fold.
(x and z being the horizontal and vertical axes, respectively), and the mass balance yields the following equation: 2 1 2c 2 1 2c 1 0: 2 2x Kzz 2x 2z Kxx 2z The horizontal Kxx and vertical Kzz hydraulic conductivities are assumed to be constant in each layer. Their ratio e Kxx =Kzz is the anisotropy coefficient. The driving force of the groundwater flow is the gradient of the piezometric head at the surface, i.e.
the slope of the water table. The elevation of the water table in the study area (Ro´nai, 1985) gives the upper boundary condition for the calculations (Fig. 4). The water table follows the ground surface closely over all the GHP, where only mild topographic variations occur. The NE end of the section (at x 168 km) was chosen at the summit of a broad topographic high. The lowest topography is reached at the River Tisza, at the starting point of the section x 0: The sides of the model domain are regarded as no flux boundaries; on the lower boundary, a vanishing vertical flux is assumed.
Table 2 Measured helium concentrations along the section (after Martel et al., 1989; Stute et al., 1992) Distance along the section (km)
Sample locality
Tiszake´cske To¨ro¨kszentmiklo´s Mezo˝tu´r ¨ rme´nyes O
0 32 36 44
1. 2. 3. 4.
57
5. Kisu´jsza´lla´s
86 99 111 114 117 132 157 160 168
6. 7. 8. 9. 10. 11. 12. 13. 14.
Pu¨spo¨klada´ny Kaba Hajdu´szoboszlo´ Nagyhegyes Balmazu´jva´ros Debrecen Gesztere´d Nyı´radony Nyı´rlugos
Depth of sampling (m)
Helium concentrations (10 12 atoms per cm 3 H2O)
211–200 239–281 161–170 71–135 258–286 93–106 239–257 192–215 182–209 120–129 110–158 85–98 138–181 144–170 244–254 110–175
49 57 101 104 160 441 649 282 277 373 85 78 6.3 2.0 2.0 2.1
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
Elevation a.s.l. (m)
190 180 160 140
Piezometric level
120 100 80 0
50
100
150
Distance (km) Fig. 4. Piezometric level at the surface along the section. It was used to calculate the upper boundary condition for the stream function. Vertical exaggeration is 440-fold.
We calculated the concentration of helium C in atoms per cm 3 of pore water. The equation describing the transport of helium is (Bear and Verruijt, 1987, p. 169): 2 nC 1 div uC 2 D grad C A; 2t
3
where t is the time, n the porosity, u {ui } u; 0; w; D the dispersion tensor, and A the in situ production of helium. D is the sum of the isotropic molecular diffusivity (D0) and the anisotropic mechanical dispersion (Marsily, 1981): Dij D0 dij 1 aT uuudij 1 aL 2 aT
ui uj ; uuu
4
where a L and a T are the longitudinal and transverse dispersivity coefficients, respectively. At the surface C is fixed as the concentration of helium in equilibrium with the atmosphere. On the sides there is no horizontal helium flux, and at the bottom we prescribed a constant incoming helium flux across a practically impermeable material. Thus, D reduces to D0 and the basal flux is: Fz 2D0
2C : 2z
5
The basal flux represents the helium that is transported from the deep crust and mantle into the aquifers. The helium is possibly introduced into the basin strata along fracture networks or fault systems (e.g. Torgersen et al., 1989). Numerical modelling (Zhao et al., 1998) shows that the concentration and isotopic composition of helium varies strongly near a fracture zone, where helium enters an aquifer. However, an earlier study of 3He/ 4He ratios and helium concentrations in the Pannonian basin did not reveal any
correlation between these quantities and the location of known basement faults (Martel, 1989). Therefore, we assumed constant basal helium flux in the modelling. We solved the equations numerically by a finitedifference method on a rectangular grid covering the model in the (x,z) plane. A typical grid used in this study consists of 168 horizontal and 70 vertical grid elements. Eq. (2) is solved by a successive overrelaxation method and Eq. (3) is solved by the alternating direction implicit scheme (Peaceman, 1977) using central differences. Only the steady-state distribution of the helium concentration is relevant for our purposes. The hydrological parameters of the studied area are poorly known. One of the aims of this paper is to estimate them. For this we chose the method of fitting the calculated concentrations to the measured ones (Table 2) by varying the model parameters. We use a weighted least-squares criterion to find the model which is closest to the measurements: X Ci 2 C m 2 i minimum; m 2 b C i i
6
where Ci is the calculated concentration at the ith observation point, Ci m the measured value, and b a constant. The weighting factor bCi m reflects the assumption that the relative error of the measurements is the same everywhere; otherwise the value of b is irrelevant to the minimization problem. Eq. (6) is solved by Nelder and Mead’s simplex method which is an iteration strategy without the need of calculating gradients in the parameter space (Press et al., 1992). In each step of the iteration, the full set of equations is to
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
191
Table 3 Parameters of the best fit model Parameter
Type
Symbol
Value
Variance a
Length of the section Depth of the section Horizontal hydraulic conductivity of the layers Quaternary Pliocene Upper Pannonian Anisotropy coefficient Diffusivity of helium Longitudinal dispersivity Longitudinal/transverse dispersivity ratio Helium flux through the bottom Ratio of internal He production
Fixed Fixed
l d
168 km 1400 m
– –
^0.64 ^0.58 Limit reached ^0.53 – ^0.45 – ^0.46 –
a b
Kxx Variable Variable Variable Variable Fixed Variable Fixed
e Kxx =Kzz D0 aL aL =aT
3:32 × 1026 ms21 0:52 × 1026 ms21 3:32 × 1026 ms21 30 1:0 × 1029 m2 s21 60 m 10
Variable Fixed b
Fz Ad= Fz 1 Ad
14:8 × 109 atoms m 22s 21 0.1, 0.2, 0.5
Variance in orders of magnitude. Ratio of internal He production was 0.1, 0.2 or 0.5, and the calculations were made for each case (see text).
be solved by the finite-difference approximation mentioned above. One of the main goals of this study is to estimate the layer conductivities. In order to make the model robust, we had to keep the number of optimization variables low and to impose a priori constraints on parameter search. Therefore, the only spatial heterogeneity assumed in the calculations was the threelayered conductivity structure, while all the other parameters, anisotropy, dispersivities, diffusion coefficient and He generation were kept constant in space throughout the model section. The main reason for this choice is that the information content of the 17 measured He concentrations is not enough to resolve too many free model parameters. The proportion of internally generated He flux (Ad) to the total He flux of the models can be defined as a
Ad : Fz 1 Ad
dispersivities at a value of aL =aT 10 (Ja¨hne et al., 1987). Preliminary experimentation with the model showed that a further constraint on the model variables proved to be useful: the conductivities of the lower layers, especially that of the third one (Upper Pannonian) should not exceed the conductivity of the topmost layer (Quaternary). This is a plausible assumption based on the fact that rock consolidation in the deeper strata normally reduces water conductivity, and this is also supported by in situ permeability measurements taken in boreholes of the study area (Szalay, 1982). Summarizing, the quantities varied by the least-squares algorithm were: horizontal conductivities Kxx of the layers, anisotropy coefficient e , longitudinal dispersivity a L and the incoming helium flux at the bottom, Fz. Table 3 lists all the geometrical and physical parameters of the model.
7
Since the amount of the internally generated helium is unknown we made three sets of calculations with values a 0:1; 0.2 and 0.5. Higher values were not investigated, because earlier studies of helium accumulation show that the main source of helium in the aquifers is the basal helium flux (Heaton, 1984; Torgersen and Clarke, 1985). Another simplification is that we fixed the ratio of longitudinal and transverse
4. Results The solution of equations results in a Darcy velocity distribution and helium concentrations. The parameters resulting in best fit between the modelled and measured helium concentrations are also given in Table 3. The characteristic features of the model are shown in Fig. 5. There is vigorous flow, shown by dense streamlines, beneath the slope of the
192
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
Fig. 5. Results of the best fit model: (a) streamlines, note the vigorous flow beneath the recharge area; (b) helium concentration for the case a 0:1; (c) helium profiles computed for the depth of 150 m (asterisks show measured and corrected data, see text); (d) computed profile of the surface heat flow. The measured data (asterisks) are taken from Do¨ve´nyi (1994). In (a) and (b) the contour values are equispaced, vertical exaggeration is 30-fold. He concentrations in (b) increase downwards with a contour interval of 10 15 atoms cm 23 H2O.
topographic high in the eastern part of the section (Fig. 5(a)). The highest Darcy velocity u 3 × 1029 m s21 approximately 9.5 cm y 21) is obtained in the Quaternary at x 120 km: The streamlines are sparser in the poor conductivity Pliocene unit. In the western part of the section there is weak groundwater flow. The distribution of the helium concentration (Fig. 5(b)) reflects the structure of groundwater flow: in the recharge zone and the zone of fast flow there is almost no helium left, but there are high concentrations in the zone of ascending flow. Fig. 5(c) shows He concentration profiles calculated for different internal He productions a 0:1; 0.2 and 0.5. The profiles of Fig. 5(c) are taken at 150 m depth. (The asterisks represent the measured He
concentrations along the model section, rescaled with a linear depth correction: as if all measurements were made at 150 m depth. This depth correction is applied to allow direct visual comparison of the calculated curves with the measurements. The least-squares fit, according to Eq. (6), was made with the measured data of Table 2.) Over most of the section the model is insensitive to the variation of a. However, in the eastern end of the section (where the concentrations are low) the increase of a results in a significant increase of the calculated concentrations, and this, in turn, causes increasing misfit of the models. So we accept a 0:1 as a best fit parameter. This low value of a gives also a posterior justification for the choice that the He production A is a single constant in the model:
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
its possible variation in correlation with the lithology of the layers would not have significant effects on the final results. Fig. 5(c) shows also that the He concentration is a quantity that is very sensitive to the change of model parameters: minor variation in the parameters results in orders of magnitude variation in the helium concentration. For example, in the western part of the section (0–70 km) the piezometric head undulates a few meters, resulting in an order of magnitude change in the helium concentration. The small increase in the piezometric level at around x 35 km causes downward water flow and in spite of its very low flux this downward flow redistributes the helium concentration. The unknown inhomogeneities in the water conducting layers should also cause significant variations in the concentration. A few data are away from the calculated curve, but still the fit is reasonably good. There is only one data point at x 0; which significantly departs from the theoretical curves. Fig. 5(d) shows heat flow data measured in the close vicinity of the section studied here. The solid line on Fig. 5(d) is a theoretical heat flow curve calculated using the flow pattern of Fig. 5(a). Groundwater flow modifies the underground temperature field by advective and dispersive heat transport. The resulting temperature field–and a consequent surface heat flow distribution–can be calculated from a heat transport equation very similar to the He transport equation (3). This calculation was done as a check of our flow model. Fig. 5(d) shows that there is good agreement between the measured and theoretical heat flow values along the section for 0 , x , 168 km; and also that the temperature field is only weakly affected by the groundwater flow. Heat flow values are a little bit lower in the recharge area at the eastern end of the section, and slightly higher than the average near x 110 km where the upward flow is strong. At the western end of the section, outside the modelled region, there is a spot with extremely high heat fluxes which was earlier identified as an area with prevailing hydrothermal convection in the upper 500 m of loose sandy sediments (Lenkey, 1993). This is the only site presently known in the GHP, where free convection occurs in the sedimentary strata. The discrepancy of the measured and modelled He concentrations at x 0 is easily understandable: vigorous convection near
193
x 0 sweeps out helium from the convection zone very efficiently. An important remark should be added concerning the parameters determined during model optimization. The least-squares procedure automatically resulted in the expected variation in the hydraulic conductivities of the layers, that is between the wellconducting Quaternary and Upper Pannonian strata, there is a relatively poor Pliocene conductor. The conductivity of this layer was found to be approximately 1/6 of the conductivity of the Quaternary unit. The conductivity of the third, Upper Pannonian, layer tended to be quite high, but due to the limit applied to this value, the optimization algorithm gave the same conductivity value for the first and third layer. The overall anisotropy was found to be 30, the longitudinal dispersivity 60 m, and the bottom He flux 14:8 × 109 atoms m 22 s 21. Variances obtained from the least-squares minimization are also found in Table 3; they were calculated using classical formulae for the covariance matrix of the variable parameters (Tarantola, 1987). The variances are given in orders of magnitude that shows the usual result that a multi-parameter model fit normally produces large parameter uncertainties. In our case, the average uncertainty of the minimization variables is around half an order of magnitude.
5. Discussion and conclusions This study gives a quantitative example of how it is possible to set up and calibrate a hydrogeological model using helium concentrations as a groundwater flow tracer. The geometry of the model and the piezometric head profile, the driving force of water flow, were constructed using existing data. The hydrogeological and physical parameters were established by model calculations, by fitting the modelled helium concentrations to the measured values. The procedure resulted in a reasonably good fit. The most important factor governing the structure of groundwater flow and helium concentration is the profile of piezometric head at the surface. Conductivity and anisotropy of the hydrogeologic units, diffusivity and dispersivity of helium are also important factors. The significance of the model is that in addition to the structure and flux
194
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
of groundwater flow, it also gives the poorly known hydrogeological parameters of the flow regime. Another important result of the modelling is the helium flux coming from greater depth. The total average helium flux leaving the system at the surface is 16:3 × 109 atoms m 22s 21. Part of this flux comes from the deep crust and mantle (Fz) and part of it is produced in situ (Ad). Fig. 5(c) shows that the helium concentration is sensitive to the ratio of the two He sources in the recharge area only. Our modelling suggests that 90% of the total flux derives from the helium entering through the lower boundary. This is Fz 14:8 × 109 atoms m 22s 21, shown in Table 3. If the concentration of U and Th in the aquifer is known, then it is possible to give an independent estimation for the rate of helium production (A). Stute and Dea´k (1989) and Stute et al. (1992) calculated this value as A 0:85 atoms cm 23s 21. From this it follows that Ad 1:2 × 109 atoms m 22s 21. If this value is subtracted from the total flux found here, then it results in Fz 15:1 × 109 atoms m 22s 21, which is 93% of the total flux, very close to the 90% suggested here. It was Martel et al. (1989), who first estimated the incoming helium flux in the GHP. They used the rate of the total discharge of the groundwater system of the GHP and the helium concentration of the discharging water. They arrived at Fz 80 × 109 atoms m 22s 21. Stute et al. (1992) used one-dimensional groundwater flow models for the typical descending, ascending and horizontal flows in the GHP and obtained Fz 0:7–4:5 × 109 atoms m 22s 21. Our estimation gives Fz 14:8 × 109 atoms m 22s 21, which is between the two previous estimations, but it is closer to the values of Stute et al. (1992). The helium flux beneath the GHP is higher than beneath the old consolidated basins of the Earth (Southern Africa: 9 × 108 to 3 × 1010 atoms m 22s 21, see Heaton, 1984; Namibia: 8 × 109 atoms m 22s 21, Australia: 9 × 109 atoms m 22s 21, see Torgersen and Clarke, 1985; Paris basin: 4 × 109 atoms m 22s 21, see Marty et al., 1993). The difference, if the uncertainty of the data is taken into account, is not significant. Analysing helium isotope ratios in groundwater samples from the Pannonian basin, Martel et al. (1989) and Stute et al. (1992) showed that 1–16% of the helium flux is originating from the mantle. The high mantle-derived helium flux is likely to be
related to the Middle Miocene extensional formation of the basin and an elevated position of the upper mantle. The decrease of radiogenic crustal helium flux due to crustal thinning is probably compensated by the increase of the mantle derived helium flux due to the elevated position of the asthenosphere beneath the Pannonian basin, resulting in almost normal helium flux. Acknowledgements This research was supported by grants from the Hungarian Research Council (Contract N. OTKA T015966 and T026633). L. Lenkey thanks the financial support given by the “Ja´nos Bolyai Research Grant”. The reviews of T. Torgersen and T. Heaton greatly contributed to the improvement of the quality of the manuscript. References Andrews, J.N., 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. Bear, J., Verruijt, A., 1987. Modeling Groundwater Flow and Pollution, D. Reidel, Dordrecht. Do¨ve´nyi, P., 1994. Geophysical investigations of the lithosphere of the Pannonian basin. PhD thesis, Eo¨tvo¨s University, Budapest. Erde´lyi, M., 1976. Outlines of the hydrodynamics and hydrochemistry of the Pannonian basin. Acta Geol. Acad. Sci. Hung. 20, 287–309. Erde´lyi, M., 1979. Hydrodynamics of the Hungarian basin. Proceedings of the Research Centre for Water Resources Development, vol. 18, pp. 1–82. Heaton, T.H.E., 1984. Rates and sources of 4He accumulation in groundwater. Hydrol. Sci. J. 29, 29–47. Horva´th, F., 1993. Towards a mechanical model for the formation of the Pannonian basin. Tectonophysics 226, 333–357. Ja¨hne, B., Heinz, G., Dietrich, W., 1987. Measurement of the diffusion coefficients of sparingly soluble gases in water. J. Geophys. Res. 92, 10 767–10 776. Juha´sz, Gy., 1991. Lithostratigraphical and sedimentological framework of the Pannonian (s.l.) sedimentary sequence in the Hungarian Plain (Alfo¨ld), Eastern Hungary. Acta Geol. Hung. 34, 53–72. Juha´sz, Gy., 1994. Comparison of the sedimentary sequences in Late Neogene subbasins in the Pannonian Basin, Hungary. Fo¨ldt. Ko¨zl. 124, 341–365. Lenkey, L., 1993. Study of the thermal anomaly at Tiszake´cske by numerical modelling of thermal convection. Magyar Geofizika 34, 30–45 in Hungarian with English abstract. Marsily, G., 1981. Hydroge´ologie Quantitative, Masson, Paris.
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195 Martel, D.J., 1989. Aspects of helium production and transport in the continents. PhD thesis, Cambridge University. Martel, D.J., Dea´k, J., Do¨ve´nyi, P., Horva´th, F., O’nions, R.K., Oxburgh, E.R., Stegena, L., Stute, M., 1989. Leakage of helium from the Pannonian basin. Nature 342, 908–912. Marty, B., Torgersen, T., Meynier, V., O’Nions, K.R., Marsily, G., 1993. Helium isotope fluxes and groundwater ages in the Dogger aquifer, Paris Basin. Water Resour. Res. 29, 1025– 1035. O’Nions, R.K., Oxburgh, E.R., 1988. Helium, volatile fluxus and the development of continental crust. Earth Planet. Sci. Lett. 90, 331–347. Peaceman, D.W., 1977. Fundamentals of Numerical Reservoir Simulation, Elsevier, Amsterdam. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P. (Eds.), 1992. Numerical Recipes. Cambridge University Press, Cambridge. Ro´nai, A., 1985. Magyarorsza´g hidrogeolo´giai te´rke´pe. In: Az Alfo¨ld negyedido˝szaki fo¨ldtana. (Hydrogeological map of Hungary. In: Quaternary of the Great Hungarian Plain.). Geol. Hung. Ser. Geol., 21. Royden, L.H.. Horva´th, F. (Eds.), 1988. The Pannonian Basin, a Study in Basin Evolution. Amer. Assoc. Petr. Geol. Mem., 45. Stute, M., Dea´k, J., 1989. Environmental isotope study ( 14C, 13C, 18 O, D, noble gases) on deep groundwater circulation systems in Hungary with reference to paleoclimate. Radiocarbon 31, 902– 918. Stute, M., Sonntag, C., Dea´k, J., Schlosser, P., 1992. Helium in deep circulating groundwater in the Great Hungarian Plain: flow
195
dynamics and crustal and mantle helium fluxes. Geochimica et Cosmochimica Acta 56, 2051–2067. ´ ., 1982. A rekonstrukcio´s szemle´letu˝ fo¨ldtani kutata´s Szalay, A leheto˝se´gei sze´nhidroge´n-perspektı´va´k elo˝rejelze´se´ben. (Possibilities of the reconstruction of basin evolution in the prediction of hydrocarbon prospects.) PhD thesis, Hungarian Academy of Sciences, Budapest, Kandida´tusi e´rtekeze´s, MTA, Budapest. Tarantola, A., 1987. Inverse Problem Theory, Elsevier, Amsterdam. Torgersen, T., Clarke, W.B., 1985. Helium accumulation in groundwater. I. An evaluation of sources and the continental flux of crustal 4He in the Great Artesian Basin, Australia. Geochimica et Cosmochimica Acta 49, 1211–1218. Torgersen, T., Ivey, G.N., 1985. Helium accumulation in groundwater, II. A model for the accumulation of the crustal 4He degassing flux. Geochimica et Cosmochimica Acta 49, 2445– 2452. Torgersen, T., Kennedy, B.M., Hiyagon, H., Chiou, K.Y., Reynolds, J.H., Clarke, W.B., 1989. Argon accumulation and the crustal degassing flux of 40Ar in the Great Artesian Basin, Australia. Earth Planet. Sci. Lett. 92, 43–56. Urbancsek, J. (Ed.), 1977. Magyarorsza´g me´lyfu´ra´su´ ku´tjainak katasztere, VII. ko¨tet. A panno´niai medence me´lyse´gi vı´zta´rozo´i. (Database of wells of Hungary, tomus VII. Deep water reservoirs of the Pannonian basin.) OVH, Vı´zgazda´lkoda´si Inte´zet. Zhao, X., Fritzel, T.L.B., Quinodoz, H.A.M., Bethke, C.M., Torgersen, T., 1998. Controls on the distribution and isotopic composition of helium in deep ground-water flows. Geology 26, 291– 294.
Modelling of helium transport in groundwater along a section in the Pannonian basin L. Cserepes, L. Lenkey* Department of Geophysics, Lora´nd Eo¨tvo¨s University, Ludovika te´r 2, 1083 Budapest, Hungary Received 14 April 1998; accepted 20 September 1999
Abstract Underground water flow in sedimentary basins controls the distribution of dissolved salts and gases, and their concentrations may therefore be used as indicators of the flow direction. Recent measurements of the 4He concentration in deep waters of the Pannonian basin have great importance in this respect. This paper presents an example of the simultaneous computation of water flow and helium distribution along a section crossing the Great Hungarian Plain. The model consists of three permeable layers. The boundaries of the layers are prescribed using geologic sections constrained by ample borehole and seismic data. The results of the finite-difference calculations are fitted to the observed helium concentrations using a least-squares algorithm that varies the model parameters. The significance of the model is that it reconstructs the structure and flux of the groundwater flow and estimates the poorly known hydrogeological parameters of the flow regime such as hydraulic conductivities, conductivity anisotropy and dispersion coefficients. The statistical uncertainty of the estimated parameters is around half an order of magnitude. An estimate of the regional average of the incoming helium flux is also obtained. The total helium flux in the Great Hungarian Plain at the surface is within the range observed in old stable continental areas of the Earth. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Groundwater flow modelling; Helium transport; Hydraulic parameters; Least-squares estimation
1. Introduction Groundwater flow in sedimentary basins controls the distribution of dissolved salts and gases. Among them, helium has received increasing attention as a tracer of groundwater movement (Andrews and Lee, 1979; Torgersen and Clarke, 1985; Torgersen and Ivey, 1985; Stute et al., 1992; Zhao et al., 1998). Helium is an inert gas, therefore its underground distribution is governed by relatively simple laws. During recharge, the helium contained in soil air is dissolved at the water table. Along the flow path * Corresponding author. E-mail address: [email protected] (L. Lenkey)
toward the discharge zone, the concentration of helium increases due to (1) the accumulation of helium produced in situ in the aquifer by radioactive decay of U and Th, and (2) by helium entering the aquifer from greater depth. Measuring the helium concentrations, we can obtain information about the direction and flux of the groundwater flow, age of the groundwater, and hydrogeological parameters of the aquifer. Most of the Earth’s helium is of radiogenic origin. The two main sources, U and Th, which produce 4He by a-particle decay, are highly concentrated in the continental crust. A minor part of the helium is primordial and escapes from the mantle at midocean ridges and through the continental lithosphere
0022-1694/99/$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(99)00158-4
186
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
in tectonically active areas, particularly those undergoing extension or volcanism (O’Nions and Oxburgh, 1988). The isotopic composition of helium is considerably different in the mantle, crust and atmosphere. This allows distinction among the radiogenic crustal, mantle-derived, and atmosphere-derived helium by measuring the ratio of helium isotopes ( 3He/ 4He). More than 200 helium concentrations and isotope ratios have been measured in groundwater samples from Hungary (Martel et al., 1989; Stute et al., 1992; see Fig. 1). The results show that 1–16% of the helium derives from the mantle. The presence of mantle-derived helium in groundwaters of the Pannonian basin is in accordance with the Neogene extensional origin of the basin (Royden and Horva´th, 1988) and with geophysical observations showing elevated upper mantle beneath the basin (Horva´th, 1993). All measured groundwater samples from the Pannonian basin contain 4He in excess of the equilibrium saturation with air. The concentration of dissolved helium varies three orders of magnitude, increasing from the recharge areas toward the discharge areas ( 1–1300 × 1012 4He atoms cm 23 H2O). The samples with the highest concentration come from the central part of the discharge zone. Application of simple one- and two-dimensional groundwater flow models for helium transport allows determination of the groundwater age, and quantification of the helium flux coming from the deeper crust and mantle (Martel et al., 1989; Stute et al., 1992). The goal of this study is to show that by calculating the 4He concentrations from a two-dimensional transport model, and fitting the results to the measured data, it is possible to reconstruct groundwater flow and He migration as well as to estimate the value of poorly known hydrogeological parameters such as
hydraulic conductivity and anisotropy. The groundwater flow and the helium concentrations are modelled along a section crossing the Great Hungarian Plain (GHP) (Fig. 1).
2. The study site The Pannonian basin is located in eastern Central Europe. It is a topographically low region which is about 400 km from N to S and 600 km from W to E and is encircled by the Eastern Alps to the west, the Western Carpathians to the north, the Apuseni Mountains to the east and the Dinarides to the south (Fig. 1). It is not a single basin, but a system of deep subbasins separated by shallow basement highs and low hills. The two largest subbasins are the GHP and the Little Hungarian Plain, where the thickness of NeogeneQuaternary sediments reaches 7–8 km in the deepest troughs. The Neogene-Quaternary sedimentary sequence in the GHP consists of alternating layers of coarsegrained and fine-grained sediments. The basin was filled by a large delta system during the Late Miocene. Based on the lithology, the delta sediments are divided into Lower and Upper Pannonian rocks. The Lower Pannonian sequence comprises clay, marl and shale-rich sediments containing sandstone bodies upward. The Upper Pannonian sequence consists of alternating clay and sandstone. The Pliocene is made of lacustrine sands and clays, the Quaternary consists of poorly consolidated fluvial sand and clayey sand deposits. The stratigraphy of the Neogene-Quaternary sedimentary sequence in the GHP is summarized in Table 1 (after Urbancsek, 1977; Juha´sz, 1991, 1994). There are two regional groundwater flow systems in
Table 1 Stratigraphy of the Neogene-Quaternary sediments in the Great Hungarian Plain (after Urbancsek, 1977; Juha´sz, 1991, 1994) Age
Facies
Lithology
Thickness (m)
Quaternary Pliocene Upper-Miocene Upper Pannonian Lower Pannonian Middle-Miocene
Terrestrial, fluvial Fluvial, lacustrine
Sand and clay Sand and clay
100–600 100–1800
Lacustrine, delta plain Delta-slope, prodelta Prodelta, deep basin Basal
Sandstone, clay, marl Clay, marl, sandstone bodies Laminated clayey marl Conglomerate, clay
800–2000 600–3000 20–800 ,100
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195 Fig. 1. Topography of the Pannonian basin and its surroundings. Thick solid line shows the location of the hydrogeological cross section shown in Fig. 2: (1) 4He concentration measurement; (×) 3He/ 4He measurement.
187
188
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
Fig. 2. Generalized hydrodynamic system of the Great Hungarian Plain after Erde´lyi (1979). For location of the section see Fig. 1. RE, Recharge areas, D, Discharge areas, Q, Quaternary, Pl, Pliocene, LP, Lower Pannonian, UP, Upper Pannonian, B, Basement. Note that the depth scale is logarithmic enhancing the flow system in the Quaternary aquifer.
the GHP. The upper flow system is in the Quaternary sediments. These unconsolidated, thick and interconnected sand layers have good conductivity allowing the circulation of large amounts of groundwater of meteoric origin. Water recharge occurs at topographic highs and discharge occurs in topographic lows (Fig. 2) (Erde´lyi, 1979). The aquifer provides water for industrial and drinking purposes. The lower flow system exists in the Upper Pannonian sediments. The Upper Pannonian sand layers constitute a multilayered aquifer. The water is mostly of meteoric origin, mixed with connate water migrating upwards from deeper layers due to compaction. Due to the high geothermal gradient of the GHP (40–508C km 21), the Upper Pannonian aquifer contains thermal water which is used for balneotherapy and heating. The Pliocene sediments form a semipervious aquitard between the Quaternary and Upper Pannonian flow systems. The Lower Pannonian sequence can be regarded as an aquiclude, because the sandstone bodies are separated by thick marl layers and the communication between the high conductivity rocks is very poor. Our model is restricted to a 168 km long section that connects a main recharge and discharge area: a topographic high in the north-eastern part of the GHP (recharge) and a topographic low in the central part of the plain (discharge). The section follows a principal flow direction identified previously by pressure data (Erde´lyi, 1979). The helium concentration data collected in the close neighbourhood of the section are listed in Table 2. All the samples were taken from the
Quaternary aquifer, in the depth range of 71–405 m, most of them between 110 and 250 m depths.
3. Model of helium transport The hydrological model (Fig. 3) is built of three permeable layers, Quaternary, Pliocene and Upper Pannonian. The Lower Pannonian unit is regarded impermeable (Erde´lyi, 1976; Stute and Dea´k, 1989; Stute et al., 1992). The boundaries of the layers are prescribed using geologic sections constrained by ample borehole and seismic data (Urbancsek, 1977). We assumed the layers to be homogeneous. Each permeable layer consists of alternating sand and clay layers. We account for this effect by assigning anisotropy to the layers. The numerical calculations are carried out in a rectangular model domain. Its bottom is prescribed at a depth d 1400 m measured from the mean level of ground surface. This depth is within the impermeable Lower Pannonian unit. The basic law governing the large scale flow of water in porous-permeable media is Darcy’s law. We assume that the density of water is constant and the rock matrix is incompressible. Then a two-dimensional flow can be derived from a scalar stream function, c (e.g. Bear and Verruijt, 1987, p. 117). The horizontal (u) and vertical (w) Darcy velocities (fluxes) are given as u2
2c ; 2z
w
2c 2x
1
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
189
Depth b.s.l. (km)
0 Quaternary K 1xx
- 200 - 400
Pliocene K 2xx
- 600
Upper Pannonian K 3xx
- 800 - 1000 - 1200
Lower Pannonian
- 1400 0
50
100
150
Distance (km) Fig. 3. Geometry of the hydrogeological units as used in the model. Boundaries of the layers are taken from hydrogeological maps (Urbancsek, 1977) and well data. K1xx, K2xx and K3xx denote the different horizontal hydraulic conductivities of the layers. Vertical exaggeration is 40-fold.
(x and z being the horizontal and vertical axes, respectively), and the mass balance yields the following equation: 2 1 2c 2 1 2c 1 0: 2 2x Kzz 2x 2z Kxx 2z The horizontal Kxx and vertical Kzz hydraulic conductivities are assumed to be constant in each layer. Their ratio e Kxx =Kzz is the anisotropy coefficient. The driving force of the groundwater flow is the gradient of the piezometric head at the surface, i.e.
the slope of the water table. The elevation of the water table in the study area (Ro´nai, 1985) gives the upper boundary condition for the calculations (Fig. 4). The water table follows the ground surface closely over all the GHP, where only mild topographic variations occur. The NE end of the section (at x 168 km) was chosen at the summit of a broad topographic high. The lowest topography is reached at the River Tisza, at the starting point of the section x 0: The sides of the model domain are regarded as no flux boundaries; on the lower boundary, a vanishing vertical flux is assumed.
Table 2 Measured helium concentrations along the section (after Martel et al., 1989; Stute et al., 1992) Distance along the section (km)
Sample locality
Tiszake´cske To¨ro¨kszentmiklo´s Mezo˝tu´r ¨ rme´nyes O
0 32 36 44
1. 2. 3. 4.
57
5. Kisu´jsza´lla´s
86 99 111 114 117 132 157 160 168
6. 7. 8. 9. 10. 11. 12. 13. 14.
Pu¨spo¨klada´ny Kaba Hajdu´szoboszlo´ Nagyhegyes Balmazu´jva´ros Debrecen Gesztere´d Nyı´radony Nyı´rlugos
Depth of sampling (m)
Helium concentrations (10 12 atoms per cm 3 H2O)
211–200 239–281 161–170 71–135 258–286 93–106 239–257 192–215 182–209 120–129 110–158 85–98 138–181 144–170 244–254 110–175
49 57 101 104 160 441 649 282 277 373 85 78 6.3 2.0 2.0 2.1
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
Elevation a.s.l. (m)
190 180 160 140
Piezometric level
120 100 80 0
50
100
150
Distance (km) Fig. 4. Piezometric level at the surface along the section. It was used to calculate the upper boundary condition for the stream function. Vertical exaggeration is 440-fold.
We calculated the concentration of helium C in atoms per cm 3 of pore water. The equation describing the transport of helium is (Bear and Verruijt, 1987, p. 169): 2 nC 1 div uC 2 D grad C A; 2t
3
where t is the time, n the porosity, u {ui } u; 0; w; D the dispersion tensor, and A the in situ production of helium. D is the sum of the isotropic molecular diffusivity (D0) and the anisotropic mechanical dispersion (Marsily, 1981): Dij D0 dij 1 aT uuudij 1 aL 2 aT
ui uj ; uuu
4
where a L and a T are the longitudinal and transverse dispersivity coefficients, respectively. At the surface C is fixed as the concentration of helium in equilibrium with the atmosphere. On the sides there is no horizontal helium flux, and at the bottom we prescribed a constant incoming helium flux across a practically impermeable material. Thus, D reduces to D0 and the basal flux is: Fz 2D0
2C : 2z
5
The basal flux represents the helium that is transported from the deep crust and mantle into the aquifers. The helium is possibly introduced into the basin strata along fracture networks or fault systems (e.g. Torgersen et al., 1989). Numerical modelling (Zhao et al., 1998) shows that the concentration and isotopic composition of helium varies strongly near a fracture zone, where helium enters an aquifer. However, an earlier study of 3He/ 4He ratios and helium concentrations in the Pannonian basin did not reveal any
correlation between these quantities and the location of known basement faults (Martel, 1989). Therefore, we assumed constant basal helium flux in the modelling. We solved the equations numerically by a finitedifference method on a rectangular grid covering the model in the (x,z) plane. A typical grid used in this study consists of 168 horizontal and 70 vertical grid elements. Eq. (2) is solved by a successive overrelaxation method and Eq. (3) is solved by the alternating direction implicit scheme (Peaceman, 1977) using central differences. Only the steady-state distribution of the helium concentration is relevant for our purposes. The hydrological parameters of the studied area are poorly known. One of the aims of this paper is to estimate them. For this we chose the method of fitting the calculated concentrations to the measured ones (Table 2) by varying the model parameters. We use a weighted least-squares criterion to find the model which is closest to the measurements: X Ci 2 C m 2 i minimum; m 2 b C i i
6
where Ci is the calculated concentration at the ith observation point, Ci m the measured value, and b a constant. The weighting factor bCi m reflects the assumption that the relative error of the measurements is the same everywhere; otherwise the value of b is irrelevant to the minimization problem. Eq. (6) is solved by Nelder and Mead’s simplex method which is an iteration strategy without the need of calculating gradients in the parameter space (Press et al., 1992). In each step of the iteration, the full set of equations is to
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
191
Table 3 Parameters of the best fit model Parameter
Type
Symbol
Value
Variance a
Length of the section Depth of the section Horizontal hydraulic conductivity of the layers Quaternary Pliocene Upper Pannonian Anisotropy coefficient Diffusivity of helium Longitudinal dispersivity Longitudinal/transverse dispersivity ratio Helium flux through the bottom Ratio of internal He production
Fixed Fixed
l d
168 km 1400 m
– –
^0.64 ^0.58 Limit reached ^0.53 – ^0.45 – ^0.46 –
a b
Kxx Variable Variable Variable Variable Fixed Variable Fixed
e Kxx =Kzz D0 aL aL =aT
3:32 × 1026 ms21 0:52 × 1026 ms21 3:32 × 1026 ms21 30 1:0 × 1029 m2 s21 60 m 10
Variable Fixed b
Fz Ad= Fz 1 Ad
14:8 × 109 atoms m 22s 21 0.1, 0.2, 0.5
Variance in orders of magnitude. Ratio of internal He production was 0.1, 0.2 or 0.5, and the calculations were made for each case (see text).
be solved by the finite-difference approximation mentioned above. One of the main goals of this study is to estimate the layer conductivities. In order to make the model robust, we had to keep the number of optimization variables low and to impose a priori constraints on parameter search. Therefore, the only spatial heterogeneity assumed in the calculations was the threelayered conductivity structure, while all the other parameters, anisotropy, dispersivities, diffusion coefficient and He generation were kept constant in space throughout the model section. The main reason for this choice is that the information content of the 17 measured He concentrations is not enough to resolve too many free model parameters. The proportion of internally generated He flux (Ad) to the total He flux of the models can be defined as a
Ad : Fz 1 Ad
dispersivities at a value of aL =aT 10 (Ja¨hne et al., 1987). Preliminary experimentation with the model showed that a further constraint on the model variables proved to be useful: the conductivities of the lower layers, especially that of the third one (Upper Pannonian) should not exceed the conductivity of the topmost layer (Quaternary). This is a plausible assumption based on the fact that rock consolidation in the deeper strata normally reduces water conductivity, and this is also supported by in situ permeability measurements taken in boreholes of the study area (Szalay, 1982). Summarizing, the quantities varied by the least-squares algorithm were: horizontal conductivities Kxx of the layers, anisotropy coefficient e , longitudinal dispersivity a L and the incoming helium flux at the bottom, Fz. Table 3 lists all the geometrical and physical parameters of the model.
7
Since the amount of the internally generated helium is unknown we made three sets of calculations with values a 0:1; 0.2 and 0.5. Higher values were not investigated, because earlier studies of helium accumulation show that the main source of helium in the aquifers is the basal helium flux (Heaton, 1984; Torgersen and Clarke, 1985). Another simplification is that we fixed the ratio of longitudinal and transverse
4. Results The solution of equations results in a Darcy velocity distribution and helium concentrations. The parameters resulting in best fit between the modelled and measured helium concentrations are also given in Table 3. The characteristic features of the model are shown in Fig. 5. There is vigorous flow, shown by dense streamlines, beneath the slope of the
192
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
Fig. 5. Results of the best fit model: (a) streamlines, note the vigorous flow beneath the recharge area; (b) helium concentration for the case a 0:1; (c) helium profiles computed for the depth of 150 m (asterisks show measured and corrected data, see text); (d) computed profile of the surface heat flow. The measured data (asterisks) are taken from Do¨ve´nyi (1994). In (a) and (b) the contour values are equispaced, vertical exaggeration is 30-fold. He concentrations in (b) increase downwards with a contour interval of 10 15 atoms cm 23 H2O.
topographic high in the eastern part of the section (Fig. 5(a)). The highest Darcy velocity u 3 × 1029 m s21 approximately 9.5 cm y 21) is obtained in the Quaternary at x 120 km: The streamlines are sparser in the poor conductivity Pliocene unit. In the western part of the section there is weak groundwater flow. The distribution of the helium concentration (Fig. 5(b)) reflects the structure of groundwater flow: in the recharge zone and the zone of fast flow there is almost no helium left, but there are high concentrations in the zone of ascending flow. Fig. 5(c) shows He concentration profiles calculated for different internal He productions a 0:1; 0.2 and 0.5. The profiles of Fig. 5(c) are taken at 150 m depth. (The asterisks represent the measured He
concentrations along the model section, rescaled with a linear depth correction: as if all measurements were made at 150 m depth. This depth correction is applied to allow direct visual comparison of the calculated curves with the measurements. The least-squares fit, according to Eq. (6), was made with the measured data of Table 2.) Over most of the section the model is insensitive to the variation of a. However, in the eastern end of the section (where the concentrations are low) the increase of a results in a significant increase of the calculated concentrations, and this, in turn, causes increasing misfit of the models. So we accept a 0:1 as a best fit parameter. This low value of a gives also a posterior justification for the choice that the He production A is a single constant in the model:
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
its possible variation in correlation with the lithology of the layers would not have significant effects on the final results. Fig. 5(c) shows also that the He concentration is a quantity that is very sensitive to the change of model parameters: minor variation in the parameters results in orders of magnitude variation in the helium concentration. For example, in the western part of the section (0–70 km) the piezometric head undulates a few meters, resulting in an order of magnitude change in the helium concentration. The small increase in the piezometric level at around x 35 km causes downward water flow and in spite of its very low flux this downward flow redistributes the helium concentration. The unknown inhomogeneities in the water conducting layers should also cause significant variations in the concentration. A few data are away from the calculated curve, but still the fit is reasonably good. There is only one data point at x 0; which significantly departs from the theoretical curves. Fig. 5(d) shows heat flow data measured in the close vicinity of the section studied here. The solid line on Fig. 5(d) is a theoretical heat flow curve calculated using the flow pattern of Fig. 5(a). Groundwater flow modifies the underground temperature field by advective and dispersive heat transport. The resulting temperature field–and a consequent surface heat flow distribution–can be calculated from a heat transport equation very similar to the He transport equation (3). This calculation was done as a check of our flow model. Fig. 5(d) shows that there is good agreement between the measured and theoretical heat flow values along the section for 0 , x , 168 km; and also that the temperature field is only weakly affected by the groundwater flow. Heat flow values are a little bit lower in the recharge area at the eastern end of the section, and slightly higher than the average near x 110 km where the upward flow is strong. At the western end of the section, outside the modelled region, there is a spot with extremely high heat fluxes which was earlier identified as an area with prevailing hydrothermal convection in the upper 500 m of loose sandy sediments (Lenkey, 1993). This is the only site presently known in the GHP, where free convection occurs in the sedimentary strata. The discrepancy of the measured and modelled He concentrations at x 0 is easily understandable: vigorous convection near
193
x 0 sweeps out helium from the convection zone very efficiently. An important remark should be added concerning the parameters determined during model optimization. The least-squares procedure automatically resulted in the expected variation in the hydraulic conductivities of the layers, that is between the wellconducting Quaternary and Upper Pannonian strata, there is a relatively poor Pliocene conductor. The conductivity of this layer was found to be approximately 1/6 of the conductivity of the Quaternary unit. The conductivity of the third, Upper Pannonian, layer tended to be quite high, but due to the limit applied to this value, the optimization algorithm gave the same conductivity value for the first and third layer. The overall anisotropy was found to be 30, the longitudinal dispersivity 60 m, and the bottom He flux 14:8 × 109 atoms m 22 s 21. Variances obtained from the least-squares minimization are also found in Table 3; they were calculated using classical formulae for the covariance matrix of the variable parameters (Tarantola, 1987). The variances are given in orders of magnitude that shows the usual result that a multi-parameter model fit normally produces large parameter uncertainties. In our case, the average uncertainty of the minimization variables is around half an order of magnitude.
5. Discussion and conclusions This study gives a quantitative example of how it is possible to set up and calibrate a hydrogeological model using helium concentrations as a groundwater flow tracer. The geometry of the model and the piezometric head profile, the driving force of water flow, were constructed using existing data. The hydrogeological and physical parameters were established by model calculations, by fitting the modelled helium concentrations to the measured values. The procedure resulted in a reasonably good fit. The most important factor governing the structure of groundwater flow and helium concentration is the profile of piezometric head at the surface. Conductivity and anisotropy of the hydrogeologic units, diffusivity and dispersivity of helium are also important factors. The significance of the model is that in addition to the structure and flux
194
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195
of groundwater flow, it also gives the poorly known hydrogeological parameters of the flow regime. Another important result of the modelling is the helium flux coming from greater depth. The total average helium flux leaving the system at the surface is 16:3 × 109 atoms m 22s 21. Part of this flux comes from the deep crust and mantle (Fz) and part of it is produced in situ (Ad). Fig. 5(c) shows that the helium concentration is sensitive to the ratio of the two He sources in the recharge area only. Our modelling suggests that 90% of the total flux derives from the helium entering through the lower boundary. This is Fz 14:8 × 109 atoms m 22s 21, shown in Table 3. If the concentration of U and Th in the aquifer is known, then it is possible to give an independent estimation for the rate of helium production (A). Stute and Dea´k (1989) and Stute et al. (1992) calculated this value as A 0:85 atoms cm 23s 21. From this it follows that Ad 1:2 × 109 atoms m 22s 21. If this value is subtracted from the total flux found here, then it results in Fz 15:1 × 109 atoms m 22s 21, which is 93% of the total flux, very close to the 90% suggested here. It was Martel et al. (1989), who first estimated the incoming helium flux in the GHP. They used the rate of the total discharge of the groundwater system of the GHP and the helium concentration of the discharging water. They arrived at Fz 80 × 109 atoms m 22s 21. Stute et al. (1992) used one-dimensional groundwater flow models for the typical descending, ascending and horizontal flows in the GHP and obtained Fz 0:7–4:5 × 109 atoms m 22s 21. Our estimation gives Fz 14:8 × 109 atoms m 22s 21, which is between the two previous estimations, but it is closer to the values of Stute et al. (1992). The helium flux beneath the GHP is higher than beneath the old consolidated basins of the Earth (Southern Africa: 9 × 108 to 3 × 1010 atoms m 22s 21, see Heaton, 1984; Namibia: 8 × 109 atoms m 22s 21, Australia: 9 × 109 atoms m 22s 21, see Torgersen and Clarke, 1985; Paris basin: 4 × 109 atoms m 22s 21, see Marty et al., 1993). The difference, if the uncertainty of the data is taken into account, is not significant. Analysing helium isotope ratios in groundwater samples from the Pannonian basin, Martel et al. (1989) and Stute et al. (1992) showed that 1–16% of the helium flux is originating from the mantle. The high mantle-derived helium flux is likely to be
related to the Middle Miocene extensional formation of the basin and an elevated position of the upper mantle. The decrease of radiogenic crustal helium flux due to crustal thinning is probably compensated by the increase of the mantle derived helium flux due to the elevated position of the asthenosphere beneath the Pannonian basin, resulting in almost normal helium flux. Acknowledgements This research was supported by grants from the Hungarian Research Council (Contract N. OTKA T015966 and T026633). L. Lenkey thanks the financial support given by the “Ja´nos Bolyai Research Grant”. The reviews of T. Torgersen and T. Heaton greatly contributed to the improvement of the quality of the manuscript. References Andrews, J.N., 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. Bear, J., Verruijt, A., 1987. Modeling Groundwater Flow and Pollution, D. Reidel, Dordrecht. Do¨ve´nyi, P., 1994. Geophysical investigations of the lithosphere of the Pannonian basin. PhD thesis, Eo¨tvo¨s University, Budapest. Erde´lyi, M., 1976. Outlines of the hydrodynamics and hydrochemistry of the Pannonian basin. Acta Geol. Acad. Sci. Hung. 20, 287–309. Erde´lyi, M., 1979. Hydrodynamics of the Hungarian basin. Proceedings of the Research Centre for Water Resources Development, vol. 18, pp. 1–82. Heaton, T.H.E., 1984. Rates and sources of 4He accumulation in groundwater. Hydrol. Sci. J. 29, 29–47. Horva´th, F., 1993. Towards a mechanical model for the formation of the Pannonian basin. Tectonophysics 226, 333–357. Ja¨hne, B., Heinz, G., Dietrich, W., 1987. Measurement of the diffusion coefficients of sparingly soluble gases in water. J. Geophys. Res. 92, 10 767–10 776. Juha´sz, Gy., 1991. Lithostratigraphical and sedimentological framework of the Pannonian (s.l.) sedimentary sequence in the Hungarian Plain (Alfo¨ld), Eastern Hungary. Acta Geol. Hung. 34, 53–72. Juha´sz, Gy., 1994. Comparison of the sedimentary sequences in Late Neogene subbasins in the Pannonian Basin, Hungary. Fo¨ldt. Ko¨zl. 124, 341–365. Lenkey, L., 1993. Study of the thermal anomaly at Tiszake´cske by numerical modelling of thermal convection. Magyar Geofizika 34, 30–45 in Hungarian with English abstract. Marsily, G., 1981. Hydroge´ologie Quantitative, Masson, Paris.
L. Cserepes, L. Lenkey / Journal of Hydrology 225 (1999) 185–195 Martel, D.J., 1989. Aspects of helium production and transport in the continents. PhD thesis, Cambridge University. Martel, D.J., Dea´k, J., Do¨ve´nyi, P., Horva´th, F., O’nions, R.K., Oxburgh, E.R., Stegena, L., Stute, M., 1989. Leakage of helium from the Pannonian basin. Nature 342, 908–912. Marty, B., Torgersen, T., Meynier, V., O’Nions, K.R., Marsily, G., 1993. Helium isotope fluxes and groundwater ages in the Dogger aquifer, Paris Basin. Water Resour. Res. 29, 1025– 1035. O’Nions, R.K., Oxburgh, E.R., 1988. Helium, volatile fluxus and the development of continental crust. Earth Planet. Sci. Lett. 90, 331–347. Peaceman, D.W., 1977. Fundamentals of Numerical Reservoir Simulation, Elsevier, Amsterdam. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P. (Eds.), 1992. Numerical Recipes. Cambridge University Press, Cambridge. Ro´nai, A., 1985. Magyarorsza´g hidrogeolo´giai te´rke´pe. In: Az Alfo¨ld negyedido˝szaki fo¨ldtana. (Hydrogeological map of Hungary. In: Quaternary of the Great Hungarian Plain.). Geol. Hung. Ser. Geol., 21. Royden, L.H.. Horva´th, F. (Eds.), 1988. The Pannonian Basin, a Study in Basin Evolution. Amer. Assoc. Petr. Geol. Mem., 45. Stute, M., Dea´k, J., 1989. Environmental isotope study ( 14C, 13C, 18 O, D, noble gases) on deep groundwater circulation systems in Hungary with reference to paleoclimate. Radiocarbon 31, 902– 918. Stute, M., Sonntag, C., Dea´k, J., Schlosser, P., 1992. Helium in deep circulating groundwater in the Great Hungarian Plain: flow
195
dynamics and crustal and mantle helium fluxes. Geochimica et Cosmochimica Acta 56, 2051–2067. ´ ., 1982. A rekonstrukcio´s szemle´letu˝ fo¨ldtani kutata´s Szalay, A leheto˝se´gei sze´nhidroge´n-perspektı´va´k elo˝rejelze´se´ben. (Possibilities of the reconstruction of basin evolution in the prediction of hydrocarbon prospects.) PhD thesis, Hungarian Academy of Sciences, Budapest, Kandida´tusi e´rtekeze´s, MTA, Budapest. Tarantola, A., 1987. Inverse Problem Theory, Elsevier, Amsterdam. Torgersen, T., Clarke, W.B., 1985. Helium accumulation in groundwater. I. An evaluation of sources and the continental flux of crustal 4He in the Great Artesian Basin, Australia. Geochimica et Cosmochimica Acta 49, 1211–1218. Torgersen, T., Ivey, G.N., 1985. Helium accumulation in groundwater, II. A model for the accumulation of the crustal 4He degassing flux. Geochimica et Cosmochimica Acta 49, 2445– 2452. Torgersen, T., Kennedy, B.M., Hiyagon, H., Chiou, K.Y., Reynolds, J.H., Clarke, W.B., 1989. Argon accumulation and the crustal degassing flux of 40Ar in the Great Artesian Basin, Australia. Earth Planet. Sci. Lett. 92, 43–56. Urbancsek, J. (Ed.), 1977. Magyarorsza´g me´lyfu´ra´su´ ku´tjainak katasztere, VII. ko¨tet. A panno´niai medence me´lyse´gi vı´zta´rozo´i. (Database of wells of Hungary, tomus VII. Deep water reservoirs of the Pannonian basin.) OVH, Vı´zgazda´lkoda´si Inte´zet. Zhao, X., Fritzel, T.L.B., Quinodoz, H.A.M., Bethke, C.M., Torgersen, T., 1998. Controls on the distribution and isotopic composition of helium in deep ground-water flows. Geology 26, 291– 294.