- Email: [email protected]
Numerical simulations of topography-induced saltwater upconing in the state of Brandenburg, Germany
Phys. Chem. Earth (B), Vol. 26, No. 4, pp. 353-359,200l
Pergamon
0 2001 Elsevier Science Ltd. All rights reserved 14W1909/01/$ - see front matter
PII: S1464-1909(01)00018-1
Numerical Simulations of Topography-Induced Brandenburg, Germany
Saltwater Upconing in the State of
A. Voss and M. Koch Department Germany
of Geohydraulics
and Engineering
Hydrology, University
of Kassel, Kurt-Wolters-Stasse3,
D-34125 Kassel,
Received 7 March 2000; accepted 25 July 2000
Abstract. Saltwater intrusion from saline formation waters in the shallow lowlands of the state of Brandenburg, Germany, is simulated using a 2D density-dependent (dd) flow and transport model. Based on the geological situation and the present-day chemical composition of the groundwaters in that region, migration scenarios with and without inclusion of density effects are modelled. We find that, due to the shallowness of the aquifer system, the surticial topography has a large effect on the flow and migration patterns and, especially, gives rises to upwelling flow underneath the discharge area of the Nuthe river. Comparing models, with and without density effects included, we then investigate possible saltwater upconing due this natural discharge flow pattern. A sensitivity study of the hydrodynamic dispersion and of the anisotropy of the aquifer is at the heart of the investigation. The results of the models show that density effects are diminishing for large values of the dispersivity and high anisotropy ratios. This means that for the management of saltwater intrusion, instead of using a complicated dd flow and transport model, it may be sufticient in most practical situation to use a passive transport model, which would have a much smaller COIIIpUtiItiO~
tiItK.0
2001 Elsevier Science Ltd. All rights reserved
1 Introduction Saltwater intrusion from saline formation waters into groundwater extraction wells has been of concern for some time in the shallow lowlands of the state of Brandenburg, Germany (Hermsdorf and Rechlin, 1997). The present and future demands for increased groundwater use in the wake of recent demographic changes in the vicinity of the new German cap&al Berlin have significantly accentuated this problem. As new groundwater pump wellfields are being planned to satisfy future drinking water demands, the need Correspondence
to: Manfred Koch
for proper pumping management strategies to prevent saltwater intrusion has become more important. Numerical modelling that may have to include sects of density-dependent (dd) transport provides a powem tool for that purpose (Bear et al., 1999). However, many of the commonly used ddcodes still lack in versatility and are, because of the coupling of flow and transport, computationally an order of magnitude more timeconsuming than those models that do not take density effects into the account (cf. Holzbecher, 1998). Therefore the modeller has to choose between an accurate representation of the basic physics, or an approximate representation that can be used to quickly obtain results which may be adequate for management purposes. The fundamental question then arises if and to what extent density effects will play a significant role as to alter the conclusions drawn for the transport from such a simplified non-density-dependent model. The present paper emphasizes some hydrogeological aspects and effects of saltwater upconing in the study region. Based on the geological situation and the present-day chemical composition of the groundwaters in a particular region where a well field is being planned, we simulate, using the 2D flow and transport model SUTRA, the migration of saline groundwaters with and without inclusion of density effects, in order to address the fundamental issue raised above. However, effects of pumping are presently not included. This would require a fully 3D dd flow and transport model. Its application to the study area is a topic of a future study.
2 Hydrogeological
Situation
The investigation area is located south-west of BerlinPotsdam in the county of Teltow-Fl&nig in the State of Brandenburg. The surface geomorphology can be described by remaining peneplains which are ancient plains created by erosion and breakthrough valleys of the river Nuthe and its
354
A. Voss and M. Koch: Numerical
Simulations
of Topography-Induced
Saltwater
Fig. 1. Map and contour lines within the investigation area
tributaries (Kalatz, 1976). The topography varies between +35 m above msl and +45 m above msl in the northeastern area of the Nuthe river and to approximately 93 m msl in the Glauer Berge in the southwest (Fig. 1). There are a total of about 80 boreholes with a maximum depth of 180 m, and an average depth of approximately 50 m. Most of the boreholes are converted to piezometers, with some of them acting as multi-level wells. A typical cross-section of the geological situation is shown in Fig. 2 (Wichmann et al., 1999). In a general sense the Rupel clay is located below the Quaternary at a depth of 100 m below sea level. The Rupel clay separates the natural saltwater of the Zechstein formation from the upper fresh water. Discontinuities in the Rupel clay layer may lead to leaking conduits, so there might be a danger of saltwater intrusion into the upper aquifer. Moreover, there is a saltdome in the vicinity of Blankensee which rises up to a depth of 60 m below msl. The potential contamination of the upper aquifer is due to vertically upward head gradients, especially where the Rupel clay is missing, and due to leaching processes at the surface of the saltdome. The lower confining layer of the main aquifer is formed by glacial till of the Elster ice-age and basin silt, as well as clayey-silty sediments of the Holstein period. The main aquifer is formed by fluvial sands of the postglacial period and of the Salic ice-age. The average thickness of the aqtier
is about 20 m. In these sands local units of glacial till with a thickness of a few meters are imbedded. In the upper part of the geological formation glacial till of the Salic ice-age is widely found which separates the total aquifer into an upper and lower part. In the area of the planned well field Grofi-Beuthen the glacial till is absent so that the thickness of the aquifer increases up to about 65 m. Above the glacial till there are sand terrasses of the Salic ice-age and fluvial sands deposited during ice-melting periods or in more recent times. These sands form the upper phreatic aquifer with a thickness of about 10 m at a maximum. Sieve analysis and pumping tests were carried out to determine the hydraulic conductivities. The experimental results yielded an average value of K = 5, lo4 m/s for the sands and of K=5- lOa m/s for the glacial till and silty sediments. The regional flow direction and the head gradients are extensively monitored through the large number of observation wells. The last measurements were taken in July, 1998. The piezometric contours are shown in Fig. 1. The main groundwater flow is directed towards the Nuthe river which flows in the SE- NE direction. The hydraulic gradient is between 0.6 ‘A in the northeastern section of the Nuthe river and 1 %o in the Nuthe valley. The riverbed of the Nuthe is barely &mated so one can assume a very good
Fig. 2. Typical cross section of the investigation area contact with the aquifer. The area precipitation is between 520 and 540 mm/a. ~ro~dwater recharge takes place and is estimated from both lysimeter data and the theoretical approach of Zieschang as of 4 l/&km’) (125 m~a) for the sand units and of 1 l/(skmz) (30 mm/a) only for the glacial till area (Miller, 1988). 3 Numerical Method 3.1 Governing equations
We use the well-known groundwater flow and transport model SUTRA (Voss, 1984) that has been applied by numerous authors to the modelling of density-dependent transport problems (cf. Koch and Bang, 1998, for’ a particular application). The basic equations solved in SUTRAare the following (cf. Voss, 1984, for details): Darcy ‘s law:
(k1
y=-
-
*OW-Pge,)
&P Fluid mass balance (groundwater flow equation) :
ps+-
at
+ szg
Notations in Eq(s). (1) to (4) are: v groundwater velocity FJU, P hydrostatic pressure PWU”% a porosity; P W? density; p m dynamic viscosity; k &T permeability; g [wr*] gravity acceleration; Qr, w3!l’l flow source/sink term; e, vertical unit vector; S = (l-n) a + r@ bT?IHJ specific storativiv, with a, j3[LTz/h11,the compressibilityof the aquifer matrix and of water, respectively; C concentration; the components of the hydrodynamic dispersiontensorDareDti=D*6ij+qI v / 6ti+(cq.-czJ Vi Vj / ) v 1 &*ITl; D* &*iT] molecular diffusiOl& &, of &I, longitudinal and transverse dispersion lengths, respectively; p. [ML31 reference density of pure water at reference concentation C, = 0. Unlike most gro~d~er models the flow equation is solved for pressure and the results are then converted back to hydraulic heads. The equation of state assumes a linear relation between concentration and density. Because the concentrations in the study area are about 5000 mg/l, i.e. -l/6 of the natural seawater concentration, the linear assumption in Eq. (4) is valid as shown among others by Koch and Bang (1992).
(1) 3.2 Setup of the 2D model To mimic gravity and density effects, a 2D vertical model is set up (see Fig. 3). The cross-section extends from southwest to northwest, perpendicular to groundwater
+ V’(&PV) = Q,
piezometric Mines over a horizontal distance of 7700 m
Solute transport equation:
and a vertical extension of 70 m. The discretisation is ~uidis~t,
Epz
+ ~pv-VC = V+pD.VC]
+ Q,
(3)
with Ax= 20 m and Ay= 2 m mesh sizes in x-
and y-directions, respectively, resulting in a model with 13,896 nodes and 13,475 elements.
Equation of state:
P(c) = PO +
(4)
As for the solution of the groundwater flow equation (Eq, 2), the letI (southwestern)and the right (northeastern) model boundaries are assumed to be constant potential (fixed or
356
35 mNN
A. Voss and M. Koch: Numerical
HyTy W74
HyTy 5l74 HyTy 25l74 HyTy %I?4
Simulations
of Topography-Induced
Saltwater
HyLwi SB7
HyTyw74
HyTyI u74
HyTy I OEt74
HyLwl23Ei67
38 mNN
Fig. 3. Model grid, boundary conditions, and the basic aquifer and aquiclude units shown with a 10 fold vertical exageration
Fig. 4. Steady state pressure distribution constant head), which means that the pressure is increasing linearly with depth. The piezometric heads chosen are 35 m above msl at the left and 38 m at the right boundary, respectively. To better constrain the model, the heads were also set constant and equal to the observed ones at the locations of the multilevel, observation wells - which show upward hydraulic gradients -. Along the upper boundary of the model only the nodes at the location of the Nuthe river were tied and set to an average stage level of 34 m. As for the solution of the solute transport equation (Eq. 3), a Dirichlet (fixed) concentration boundary condition C = C, is assumed over that section of the aquifer bottom where the assumed breach in the Rupel clay formation (see Fig. 2) may allow vertical leakage of brine into the upper aquifer. The width of this leakage section is 200 m and mimics a continuous source of saltwater with a concentration of C, = 5000 mg/l. It should be noted that the proper choice of the concentration boundary condition to simulate leaching of salt from a dome into groundwater is by no means a trivial task, and it may have a significant effect on the resulting concentration field. This was one of the lessons learned from the HYDROCOIN project (Konikow et al., 1997). In particular, the Dirichlet boundary condition C= C, appears to introduce too much salt into the flow field, since it releases salt by both dispersion and advection. We will not discuss this issue further, but we will be consistent with the above choice of the boundary conditions in the following model simulations. To simulate the process of saltwater upconing in a gravitydriven phreatic aquifer, all models are set-up in steady-state mode for the pressure, but in transient mode for the solution of the solute transport equation.
4 Model results 4.1 Isotropic models with small dispersivities At first model simulations are carried out for isotropic conditions, i.e. the ratio of the horizontal and vertical permeability $ / k, = 1. The dispersivities are set to a, = 10 m and a, = 1 m which, for the mesh size used, are the minimum values to satisfy the Peclet-number criterion for numerical stability of the SUTRA code (Pe = Ax/u < 2) (Voss, 1984). The steady-state pressure distribution for this model is shown in Fig. 4. The vertical distribution of the piezometric isolines with a reduced separation of the latter underneath the Nuthe river indicates that the groundwater flows mainly in horizontal direction, with a small upward flow component, towards the river. This shows clearly the dominant and well-known effect of the sloping surf&x topography which drives a hydraulic gradient in the direction of the river bed (Freeze and Cherry, 1979). Figure 5 shows the distribution of the relative salt concentration in the aquifer for the density-independent (tracer) (top panel) and the density-dependent (dd) (bottom panel) model, after a quasi steady-state has been reached (10000 days). For a better display of the results, only the most interesting part of the aquifer between x = 1500 m and 2 100 m above the saltwater leaking breach of the Rupel clay is shown. A visual comparison of the two plumes shows clearly the stabilizing effect of the dd model (Fig.5, bottom panel): The bulk of the salt plume has risen less for the dd model than for the tracer model (Fig.5, top panel), i.e. upconing of the saltwater interface takes place at a slower rate for the former.
A. Voss and M. Koch: Numerical Simulations of Topography-Induced Saltwater Nuthe River I Tracer
---------
,
Nuthe River
I
I
------------
Nuthe River Density-dependent
Nuthe River
1
- _I_ 1500
16W
1700
1800
IWO
2000
-
-
21w
Fig. 5. Relative concentrations of intrusion plume after 10000 days; isotropic case; a, = 10 m, a, = 1 m; top panel: tracer; bottom panel: density-dependent (dd) case
Fig; 6. Relative concentrations of intrusion plume after 10000 days; isotropic case; large dispersivity q = 20 m, a, = 2 m; top panel: tracer; bottom panel: density-dependent (dd) case
Numerical stability has been investigated by varying the size of the time step. For example, using 10 days instead of 0.1 days for the latter, no differences in both the pressure and the concentration solutions were found. However, increasing the initial concentration up to the solubility limit of salt and using a big time step led to numerical instabilities.
stabilizing effect of the density on the saltwater upconing. However, the differences between the tracer- and the ddcase are now much smaller than for the previous two models with reduced dispersivities (Fig. 5). Thus, one can conclude that density effects are not so important for large values of the hydrodynamic dispersion. Breakthrough Curves: HyTy55l74
4.2 Influence of hydrodynamic dispersion To investigate the influence of the hydrodynamic dispersion, in consideration of the spatial scale effect of the latter (e.g. Gelhar et al., 1992), we did simulations with higher dispersivityvalues of a, = 20 m and or = 2 m. Each of the other parameters was kept unchanged. As expected for the higher dispersivity values, Fig. 6 for the tracer- (top panel) and the dd- (bottom panel) flow simulations shows more spreading of the salt plume in horizontal, as well as in vertical direction. This means that higher concentrations in the middle part of the aquifer are obtained which may lead to a higher potential danger if groundwater is extracted. Comparing the two panels of Fig. 6, one notes again a
0.00 0
2mO
40m
6000 mm
8Om
loQ1
WI
Fig. 7. Simulated concentration breakthrough curves at well HyTy W74 for models of Fig(s). 5 and 6.
358
A. Voss and M. Koch: Numerical Simulations of Topography-Induced Saltwater Nuthe River
A mom quantitative picture of the conclusions drawn above emerges from breakthrough curves of the saltwater front at the particular location of the monitor well HyTy 5974 (see Fig. 5 for location). It is located 18 m above the concentration boundary. One observes from Fig. 7 for the two models discussed that (1) steady state is reached after a relative short time of -2000 days and (2) the relative concentrations for the dd- plumes are always below those of the tracer phunes.
4.3 Influence of anisotropy Finally, we investigate the influence of anisotropy which is often present in fluvial aquifers, similar to the one studied here due to the natural deposition of sediments. It is well known that the vertical permeability k, is smaller than the horizontal one k, in such an environment. In the following models, k, has been decreased by a factor of 10, with k, kept invariant. The hydrodynamic dispersivities are set back againto%= lOmandq= 1 m. Figure 8 show the results of the tracer- (top panel) and of the dd- (bottom panel) simulations. For the first case one notes a horizontal spreading of the salt plume, similar to the two previous models with high dispersivity values (Fig. 6). As a result of the lower vertical permeability k, , the horizontal flow will be more accelerated to preserve the water mass-balance. This also leads to higher horizontal and tmnsve& dispersion. A stabilizing density effect can hardly be recognized, unlike the case for the isotropic models of the previous section. Other models with an even higher anisotropy ratio masked the density effects even more, and thus one might argue that under these conditions the passive tracer modelling approach is completely sufficient to properly mimic the physics of saltwater upconing.
5 Conclusions and outlook Based on the geological situation and the present-day chemical composition of the groundwaters in an area of the state of Brandenburg where a well field is being planned, we have simulated, using the 2D flow and transport model SUTRA, the effects of density on the migration of saline groundwaters. It is found that, due to the shallowness of the aquifer system, the surllcial topography has a large effect on the flow and migration patterns and, especially, gives rises to upwelling flow underneath the discharge area of the Nuthe river. Comparing models, with (dd) and without inclusion of density effects (tracer), we then investigated possible saltwater upconing due this natural discharge flow pattern. As for the need of using a computationally time-consuming density-dependent (dd) flow and transport model, instead of a faster tracer model, the following common conclusions can be drawn: 1) Generally, density has a stabilizing effect on the vertical
Nuthe Dencity-
---------+ I
River
pendent
- I-------
- -
-‘-
-
-
-
I HyTyY74 - - _I_ _ _
’
it
Fig. 8. Relative concentrations of intrusion plume after 10000 days; anisotropic case with k, / k, = 10/l; a, = 10 m, a, = 1 m; top panel: tracer; bottom panel: densitydependent (dd) case uprising of saline formation waters, i.e. the migration rate is smaller for the dd-models. 2) When increasing the hydrodynamical dispersivities, the differences between the tracer- and the dd models are smaller, i.e. density effects appear to be less important in such situations. 3) When including anisotropy into the models, density effects are masked even more and essentially no differences between the two modelling approaches are observed. From conclusions 2) and 3) one may infer that the application of a complicated dd- flow and transport model, instead of a passive (tracer) transport model, may not be warranted in practice. However, the situation might be somewhat different when the upconing due to pumping from the planned well field is taken into account. This requires the use of a three-dimensional density-dependent groundwater flow and transport model. Models for that purpose are currently being implemented and tested by the authors.
Acknowfedgement. This work bas been fun&d by tbe Landesumweltamt Brandenburg, Potsdam.
A. Voss and M. Koch: Numerical Simulations of Topography-Induced Saltwater
Koch, M. and 0. Zhang, Numerical modelling and v
References Hear, J., Cheng, AH.-D., Sorek, S., Ouazar D., and Herrera, I., Seawater intrusion in coastal aquifers - Concepts, Methods and Practices, Kluwer Academic F’ubl.,Dordrecbt, The Netherlands, 640 pp., 1999. Freeze, Rk
and Cherry, LA, Groundwater,
Prentice Hall, Engelwocxl Cliffs,
N.J., 1979 Gelbar, L.W., Welty, C., and Rehfeldt, K.R., A critical review of data on field-scale dispersion in aquifers, WaterResour. Rex, 28,1955-1974,1992. Hermsdo& A and R&in, B., Hydrogeologische Verhilltnisse. In: S&tier, J.H.: Geologis von Berlin und Brandenburg, Nr. 4: Potsdam und Umgebung, Selbstverlag Geowiss. von Berlin und Brandenburg e.V., Berlin, pp. 138-142, 1997. Holzbecher, E.O., Modelling densipdrivenflow in porous media, Springer Verlag, Berlin, Heidelberg, 286 pp., 1998. Kalatq R, Hydrogeologischer Ergebnislxricht mit Vorratsnachweis Tbyrow 1974175, XQ3 Hydrogeologie Nordhausen, unverijffentlicht, 1976. Koch, M. and Zhang G., Numerical simulation of the effects of variable density in a contaminant plume, Groundwater, 5.73 l-742, 1992.
359
of saltwater seepage f?om coastal brackish canals in southeast Florida, In: Envtronmental Coastal Regions, Brebbia, C.A (ed.), WIT Press, So&ton, UK, pp. 395-404,1998. Konikow, L.F., Sanford, W.E., and Campbell, P.J., Comtant-conce&ation boundary condition: Lessons from the HYDROCOIN variable-density groundwater ben&ma& problem, WoterResour. Res., 33,2253-2261.1997. Hydrogeologischer Ergebnisbericht mit Miiller, A., Grundwasservorratsberechnung - Detailerkundung Thyrow, VEB Hydrogeologie Nordhausen, BT Berlin, 1988. Voss, CL, A ftite-elemeni simulation model for salumt&um&naU, fluiddensity- dependent ground-water flow with energy transport a chemicallyreactive single-species solute transport U.S. Geological Survey WaterResources Investigations Report 84-4369,#9p., Reston, VI, 1984. Wichmann, K., Koch, M., Voss, A et al., Die Dynamik der Salz/SiU3wassergrenze im Grundwasser als Kriterium der langfiistigen Sicherheit der der Trinkw asserressourcen im Land Brandenburg. Abschhu&xicht Arbeitsgruppe Grundwasserversalzung G&be&m, Londesumweltamt Brandenburg, Potsdam, 1999.
Pergamon
0 2001 Elsevier Science Ltd. All rights reserved 14W1909/01/$ - see front matter
PII: S1464-1909(01)00018-1
Numerical Simulations of Topography-Induced Brandenburg, Germany
Saltwater Upconing in the State of
A. Voss and M. Koch Department Germany
of Geohydraulics
and Engineering
Hydrology, University
of Kassel, Kurt-Wolters-Stasse3,
D-34125 Kassel,
Received 7 March 2000; accepted 25 July 2000
Abstract. Saltwater intrusion from saline formation waters in the shallow lowlands of the state of Brandenburg, Germany, is simulated using a 2D density-dependent (dd) flow and transport model. Based on the geological situation and the present-day chemical composition of the groundwaters in that region, migration scenarios with and without inclusion of density effects are modelled. We find that, due to the shallowness of the aquifer system, the surticial topography has a large effect on the flow and migration patterns and, especially, gives rises to upwelling flow underneath the discharge area of the Nuthe river. Comparing models, with and without density effects included, we then investigate possible saltwater upconing due this natural discharge flow pattern. A sensitivity study of the hydrodynamic dispersion and of the anisotropy of the aquifer is at the heart of the investigation. The results of the models show that density effects are diminishing for large values of the dispersivity and high anisotropy ratios. This means that for the management of saltwater intrusion, instead of using a complicated dd flow and transport model, it may be sufticient in most practical situation to use a passive transport model, which would have a much smaller COIIIpUtiItiO~
tiItK.0
2001 Elsevier Science Ltd. All rights reserved
1 Introduction Saltwater intrusion from saline formation waters into groundwater extraction wells has been of concern for some time in the shallow lowlands of the state of Brandenburg, Germany (Hermsdorf and Rechlin, 1997). The present and future demands for increased groundwater use in the wake of recent demographic changes in the vicinity of the new German cap&al Berlin have significantly accentuated this problem. As new groundwater pump wellfields are being planned to satisfy future drinking water demands, the need Correspondence
to: Manfred Koch
for proper pumping management strategies to prevent saltwater intrusion has become more important. Numerical modelling that may have to include sects of density-dependent (dd) transport provides a powem tool for that purpose (Bear et al., 1999). However, many of the commonly used ddcodes still lack in versatility and are, because of the coupling of flow and transport, computationally an order of magnitude more timeconsuming than those models that do not take density effects into the account (cf. Holzbecher, 1998). Therefore the modeller has to choose between an accurate representation of the basic physics, or an approximate representation that can be used to quickly obtain results which may be adequate for management purposes. The fundamental question then arises if and to what extent density effects will play a significant role as to alter the conclusions drawn for the transport from such a simplified non-density-dependent model. The present paper emphasizes some hydrogeological aspects and effects of saltwater upconing in the study region. Based on the geological situation and the present-day chemical composition of the groundwaters in a particular region where a well field is being planned, we simulate, using the 2D flow and transport model SUTRA, the migration of saline groundwaters with and without inclusion of density effects, in order to address the fundamental issue raised above. However, effects of pumping are presently not included. This would require a fully 3D dd flow and transport model. Its application to the study area is a topic of a future study.
2 Hydrogeological
Situation
The investigation area is located south-west of BerlinPotsdam in the county of Teltow-Fl&nig in the State of Brandenburg. The surface geomorphology can be described by remaining peneplains which are ancient plains created by erosion and breakthrough valleys of the river Nuthe and its
354
A. Voss and M. Koch: Numerical
Simulations
of Topography-Induced
Saltwater
Fig. 1. Map and contour lines within the investigation area
tributaries (Kalatz, 1976). The topography varies between +35 m above msl and +45 m above msl in the northeastern area of the Nuthe river and to approximately 93 m msl in the Glauer Berge in the southwest (Fig. 1). There are a total of about 80 boreholes with a maximum depth of 180 m, and an average depth of approximately 50 m. Most of the boreholes are converted to piezometers, with some of them acting as multi-level wells. A typical cross-section of the geological situation is shown in Fig. 2 (Wichmann et al., 1999). In a general sense the Rupel clay is located below the Quaternary at a depth of 100 m below sea level. The Rupel clay separates the natural saltwater of the Zechstein formation from the upper fresh water. Discontinuities in the Rupel clay layer may lead to leaking conduits, so there might be a danger of saltwater intrusion into the upper aquifer. Moreover, there is a saltdome in the vicinity of Blankensee which rises up to a depth of 60 m below msl. The potential contamination of the upper aquifer is due to vertically upward head gradients, especially where the Rupel clay is missing, and due to leaching processes at the surface of the saltdome. The lower confining layer of the main aquifer is formed by glacial till of the Elster ice-age and basin silt, as well as clayey-silty sediments of the Holstein period. The main aquifer is formed by fluvial sands of the postglacial period and of the Salic ice-age. The average thickness of the aqtier
is about 20 m. In these sands local units of glacial till with a thickness of a few meters are imbedded. In the upper part of the geological formation glacial till of the Salic ice-age is widely found which separates the total aquifer into an upper and lower part. In the area of the planned well field Grofi-Beuthen the glacial till is absent so that the thickness of the aquifer increases up to about 65 m. Above the glacial till there are sand terrasses of the Salic ice-age and fluvial sands deposited during ice-melting periods or in more recent times. These sands form the upper phreatic aquifer with a thickness of about 10 m at a maximum. Sieve analysis and pumping tests were carried out to determine the hydraulic conductivities. The experimental results yielded an average value of K = 5, lo4 m/s for the sands and of K=5- lOa m/s for the glacial till and silty sediments. The regional flow direction and the head gradients are extensively monitored through the large number of observation wells. The last measurements were taken in July, 1998. The piezometric contours are shown in Fig. 1. The main groundwater flow is directed towards the Nuthe river which flows in the SE- NE direction. The hydraulic gradient is between 0.6 ‘A in the northeastern section of the Nuthe river and 1 %o in the Nuthe valley. The riverbed of the Nuthe is barely &mated so one can assume a very good
Fig. 2. Typical cross section of the investigation area contact with the aquifer. The area precipitation is between 520 and 540 mm/a. ~ro~dwater recharge takes place and is estimated from both lysimeter data and the theoretical approach of Zieschang as of 4 l/&km’) (125 m~a) for the sand units and of 1 l/(skmz) (30 mm/a) only for the glacial till area (Miller, 1988). 3 Numerical Method 3.1 Governing equations
We use the well-known groundwater flow and transport model SUTRA (Voss, 1984) that has been applied by numerous authors to the modelling of density-dependent transport problems (cf. Koch and Bang, 1998, for’ a particular application). The basic equations solved in SUTRAare the following (cf. Voss, 1984, for details): Darcy ‘s law:
(k1
y=-
-
*OW-Pge,)
&P Fluid mass balance (groundwater flow equation) :
ps+-
at
+ szg
Notations in Eq(s). (1) to (4) are: v groundwater velocity FJU, P hydrostatic pressure PWU”% a porosity; P W? density; p m dynamic viscosity; k &T permeability; g [wr*] gravity acceleration; Qr, w3!l’l flow source/sink term; e, vertical unit vector; S = (l-n) a + r@ bT?IHJ specific storativiv, with a, j3[LTz/h11,the compressibilityof the aquifer matrix and of water, respectively; C concentration; the components of the hydrodynamic dispersiontensorDareDti=D*6ij+qI v / 6ti+(cq.-czJ Vi Vj / ) v 1 &*ITl; D* &*iT] molecular diffusiOl& &, of &I, longitudinal and transverse dispersion lengths, respectively; p. [ML31 reference density of pure water at reference concentation C, = 0. Unlike most gro~d~er models the flow equation is solved for pressure and the results are then converted back to hydraulic heads. The equation of state assumes a linear relation between concentration and density. Because the concentrations in the study area are about 5000 mg/l, i.e. -l/6 of the natural seawater concentration, the linear assumption in Eq. (4) is valid as shown among others by Koch and Bang (1992).
(1) 3.2 Setup of the 2D model To mimic gravity and density effects, a 2D vertical model is set up (see Fig. 3). The cross-section extends from southwest to northwest, perpendicular to groundwater
+ V’(&PV) = Q,
piezometric Mines over a horizontal distance of 7700 m
Solute transport equation:
and a vertical extension of 70 m. The discretisation is ~uidis~t,
Epz
+ ~pv-VC = V+pD.VC]
+ Q,
(3)
with Ax= 20 m and Ay= 2 m mesh sizes in x-
and y-directions, respectively, resulting in a model with 13,896 nodes and 13,475 elements.
Equation of state:
P(c) = PO +
(4)
As for the solution of the groundwater flow equation (Eq, 2), the letI (southwestern)and the right (northeastern) model boundaries are assumed to be constant potential (fixed or
356
35 mNN
A. Voss and M. Koch: Numerical
HyTy W74
HyTy 5l74 HyTy 25l74 HyTy %I?4
Simulations
of Topography-Induced
Saltwater
HyLwi SB7
HyTyw74
HyTyI u74
HyTy I OEt74
HyLwl23Ei67
38 mNN
Fig. 3. Model grid, boundary conditions, and the basic aquifer and aquiclude units shown with a 10 fold vertical exageration
Fig. 4. Steady state pressure distribution constant head), which means that the pressure is increasing linearly with depth. The piezometric heads chosen are 35 m above msl at the left and 38 m at the right boundary, respectively. To better constrain the model, the heads were also set constant and equal to the observed ones at the locations of the multilevel, observation wells - which show upward hydraulic gradients -. Along the upper boundary of the model only the nodes at the location of the Nuthe river were tied and set to an average stage level of 34 m. As for the solution of the solute transport equation (Eq. 3), a Dirichlet (fixed) concentration boundary condition C = C, is assumed over that section of the aquifer bottom where the assumed breach in the Rupel clay formation (see Fig. 2) may allow vertical leakage of brine into the upper aquifer. The width of this leakage section is 200 m and mimics a continuous source of saltwater with a concentration of C, = 5000 mg/l. It should be noted that the proper choice of the concentration boundary condition to simulate leaching of salt from a dome into groundwater is by no means a trivial task, and it may have a significant effect on the resulting concentration field. This was one of the lessons learned from the HYDROCOIN project (Konikow et al., 1997). In particular, the Dirichlet boundary condition C= C, appears to introduce too much salt into the flow field, since it releases salt by both dispersion and advection. We will not discuss this issue further, but we will be consistent with the above choice of the boundary conditions in the following model simulations. To simulate the process of saltwater upconing in a gravitydriven phreatic aquifer, all models are set-up in steady-state mode for the pressure, but in transient mode for the solution of the solute transport equation.
4 Model results 4.1 Isotropic models with small dispersivities At first model simulations are carried out for isotropic conditions, i.e. the ratio of the horizontal and vertical permeability $ / k, = 1. The dispersivities are set to a, = 10 m and a, = 1 m which, for the mesh size used, are the minimum values to satisfy the Peclet-number criterion for numerical stability of the SUTRA code (Pe = Ax/u < 2) (Voss, 1984). The steady-state pressure distribution for this model is shown in Fig. 4. The vertical distribution of the piezometric isolines with a reduced separation of the latter underneath the Nuthe river indicates that the groundwater flows mainly in horizontal direction, with a small upward flow component, towards the river. This shows clearly the dominant and well-known effect of the sloping surf&x topography which drives a hydraulic gradient in the direction of the river bed (Freeze and Cherry, 1979). Figure 5 shows the distribution of the relative salt concentration in the aquifer for the density-independent (tracer) (top panel) and the density-dependent (dd) (bottom panel) model, after a quasi steady-state has been reached (10000 days). For a better display of the results, only the most interesting part of the aquifer between x = 1500 m and 2 100 m above the saltwater leaking breach of the Rupel clay is shown. A visual comparison of the two plumes shows clearly the stabilizing effect of the dd model (Fig.5, bottom panel): The bulk of the salt plume has risen less for the dd model than for the tracer model (Fig.5, top panel), i.e. upconing of the saltwater interface takes place at a slower rate for the former.
A. Voss and M. Koch: Numerical Simulations of Topography-Induced Saltwater Nuthe River I Tracer
---------
,
Nuthe River
I
I
------------
Nuthe River Density-dependent
Nuthe River
1
- _I_ 1500
16W
1700
1800
IWO
2000
-
-
21w
Fig. 5. Relative concentrations of intrusion plume after 10000 days; isotropic case; a, = 10 m, a, = 1 m; top panel: tracer; bottom panel: density-dependent (dd) case
Fig; 6. Relative concentrations of intrusion plume after 10000 days; isotropic case; large dispersivity q = 20 m, a, = 2 m; top panel: tracer; bottom panel: density-dependent (dd) case
Numerical stability has been investigated by varying the size of the time step. For example, using 10 days instead of 0.1 days for the latter, no differences in both the pressure and the concentration solutions were found. However, increasing the initial concentration up to the solubility limit of salt and using a big time step led to numerical instabilities.
stabilizing effect of the density on the saltwater upconing. However, the differences between the tracer- and the ddcase are now much smaller than for the previous two models with reduced dispersivities (Fig. 5). Thus, one can conclude that density effects are not so important for large values of the hydrodynamic dispersion. Breakthrough Curves: HyTy55l74
4.2 Influence of hydrodynamic dispersion To investigate the influence of the hydrodynamic dispersion, in consideration of the spatial scale effect of the latter (e.g. Gelhar et al., 1992), we did simulations with higher dispersivityvalues of a, = 20 m and or = 2 m. Each of the other parameters was kept unchanged. As expected for the higher dispersivity values, Fig. 6 for the tracer- (top panel) and the dd- (bottom panel) flow simulations shows more spreading of the salt plume in horizontal, as well as in vertical direction. This means that higher concentrations in the middle part of the aquifer are obtained which may lead to a higher potential danger if groundwater is extracted. Comparing the two panels of Fig. 6, one notes again a
0.00 0
2mO
40m
6000 mm
8Om
loQ1
WI
Fig. 7. Simulated concentration breakthrough curves at well HyTy W74 for models of Fig(s). 5 and 6.
358
A. Voss and M. Koch: Numerical Simulations of Topography-Induced Saltwater Nuthe River
A mom quantitative picture of the conclusions drawn above emerges from breakthrough curves of the saltwater front at the particular location of the monitor well HyTy 5974 (see Fig. 5 for location). It is located 18 m above the concentration boundary. One observes from Fig. 7 for the two models discussed that (1) steady state is reached after a relative short time of -2000 days and (2) the relative concentrations for the dd- plumes are always below those of the tracer phunes.
4.3 Influence of anisotropy Finally, we investigate the influence of anisotropy which is often present in fluvial aquifers, similar to the one studied here due to the natural deposition of sediments. It is well known that the vertical permeability k, is smaller than the horizontal one k, in such an environment. In the following models, k, has been decreased by a factor of 10, with k, kept invariant. The hydrodynamic dispersivities are set back againto%= lOmandq= 1 m. Figure 8 show the results of the tracer- (top panel) and of the dd- (bottom panel) simulations. For the first case one notes a horizontal spreading of the salt plume, similar to the two previous models with high dispersivity values (Fig. 6). As a result of the lower vertical permeability k, , the horizontal flow will be more accelerated to preserve the water mass-balance. This also leads to higher horizontal and tmnsve& dispersion. A stabilizing density effect can hardly be recognized, unlike the case for the isotropic models of the previous section. Other models with an even higher anisotropy ratio masked the density effects even more, and thus one might argue that under these conditions the passive tracer modelling approach is completely sufficient to properly mimic the physics of saltwater upconing.
5 Conclusions and outlook Based on the geological situation and the present-day chemical composition of the groundwaters in an area of the state of Brandenburg where a well field is being planned, we have simulated, using the 2D flow and transport model SUTRA, the effects of density on the migration of saline groundwaters. It is found that, due to the shallowness of the aquifer system, the surllcial topography has a large effect on the flow and migration patterns and, especially, gives rises to upwelling flow underneath the discharge area of the Nuthe river. Comparing models, with (dd) and without inclusion of density effects (tracer), we then investigated possible saltwater upconing due this natural discharge flow pattern. As for the need of using a computationally time-consuming density-dependent (dd) flow and transport model, instead of a faster tracer model, the following common conclusions can be drawn: 1) Generally, density has a stabilizing effect on the vertical
Nuthe Dencity-
---------+ I
River
pendent
- I-------
- -
-‘-
-
-
-
I HyTyY74 - - _I_ _ _
’
it
Fig. 8. Relative concentrations of intrusion plume after 10000 days; anisotropic case with k, / k, = 10/l; a, = 10 m, a, = 1 m; top panel: tracer; bottom panel: densitydependent (dd) case uprising of saline formation waters, i.e. the migration rate is smaller for the dd-models. 2) When increasing the hydrodynamical dispersivities, the differences between the tracer- and the dd models are smaller, i.e. density effects appear to be less important in such situations. 3) When including anisotropy into the models, density effects are masked even more and essentially no differences between the two modelling approaches are observed. From conclusions 2) and 3) one may infer that the application of a complicated dd- flow and transport model, instead of a passive (tracer) transport model, may not be warranted in practice. However, the situation might be somewhat different when the upconing due to pumping from the planned well field is taken into account. This requires the use of a three-dimensional density-dependent groundwater flow and transport model. Models for that purpose are currently being implemented and tested by the authors.
Acknowfedgement. This work bas been fun&d by tbe Landesumweltamt Brandenburg, Potsdam.
A. Voss and M. Koch: Numerical Simulations of Topography-Induced Saltwater
Koch, M. and 0. Zhang, Numerical modelling and v
References Hear, J., Cheng, AH.-D., Sorek, S., Ouazar D., and Herrera, I., Seawater intrusion in coastal aquifers - Concepts, Methods and Practices, Kluwer Academic F’ubl.,Dordrecbt, The Netherlands, 640 pp., 1999. Freeze, Rk
and Cherry, LA, Groundwater,
Prentice Hall, Engelwocxl Cliffs,
N.J., 1979 Gelbar, L.W., Welty, C., and Rehfeldt, K.R., A critical review of data on field-scale dispersion in aquifers, WaterResour. Rex, 28,1955-1974,1992. Hermsdo& A and R&in, B., Hydrogeologische Verhilltnisse. In: S&tier, J.H.: Geologis von Berlin und Brandenburg, Nr. 4: Potsdam und Umgebung, Selbstverlag Geowiss. von Berlin und Brandenburg e.V., Berlin, pp. 138-142, 1997. Holzbecher, E.O., Modelling densipdrivenflow in porous media, Springer Verlag, Berlin, Heidelberg, 286 pp., 1998. Kalatq R, Hydrogeologischer Ergebnislxricht mit Vorratsnachweis Tbyrow 1974175, XQ3 Hydrogeologie Nordhausen, unverijffentlicht, 1976. Koch, M. and Zhang G., Numerical simulation of the effects of variable density in a contaminant plume, Groundwater, 5.73 l-742, 1992.
359
of saltwater seepage f?om coastal brackish canals in southeast Florida, In: Envtronmental Coastal Regions, Brebbia, C.A (ed.), WIT Press, So&ton, UK, pp. 395-404,1998. Konikow, L.F., Sanford, W.E., and Campbell, P.J., Comtant-conce&ation boundary condition: Lessons from the HYDROCOIN variable-density groundwater ben&ma& problem, WoterResour. Res., 33,2253-2261.1997. Hydrogeologischer Ergebnisbericht mit Miiller, A., Grundwasservorratsberechnung - Detailerkundung Thyrow, VEB Hydrogeologie Nordhausen, BT Berlin, 1988. Voss, CL, A ftite-elemeni simulation model for salumt&um&naU, fluiddensity- dependent ground-water flow with energy transport a chemicallyreactive single-species solute transport U.S. Geological Survey WaterResources Investigations Report 84-4369,#9p., Reston, VI, 1984. Wichmann, K., Koch, M., Voss, A et al., Die Dynamik der Salz/SiU3wassergrenze im Grundwasser als Kriterium der langfiistigen Sicherheit der der Trinkw asserressourcen im Land Brandenburg. Abschhu&xicht Arbeitsgruppe Grundwasserversalzung G&be&m, Londesumweltamt Brandenburg, Potsdam, 1999.