Mathematical Modelling of the Groundwater and Heat Flow in the Complicated Hydrogeology Structures

Mathematical Modelling of the Groundwater and Heat Flow in the Complicated Hydrogeology Structures

Baier J1 Uhlík J2 Datel J.V.3

1Faculty of Environmental Sciences, Czech University of Life Sciences Prague, Kamýcká 129, Praha 6

Suchdol, 16521, Czech Republic baier@centrum.cz

2PROGEO l.t.d., Tiché Údolí 113, Roztoky u Prahy, 25263, Czech Republic

3Department of Hydrogeology, Faculty of Science, Charles University in Prague, Albertov street 6, 12843 Prague 2, Czech Republic

Abstract

The Benesov-Usti aquifer system of the Bohemian Cretaceous Basin constitutes the largest accumulation of moderately thermal groundwater (25 – 35°C) in the Czech Republic. The main goal of the research project was the determination of available yield of groundwater resources in the areas of Decin and Usti and effect of withdrawals on the quantity and temperature of the waters. The area is strongly influenced by Saxonian orogeny. Common approaches to the simulation of groundwater flow using MODFLOW packages BCF and LPF fails. Suitable simulation of the hydraulic function of faults and heat transfer simulations is only possible with discretization of the model layers independent of the location of the aquifers and aquitards. The MODFLOW using HUF package and MT3DMS model were evaluated as satisfied tools for the calculation of groundwater and heat flow. Potential for thermal groundwater abstraction up to 275 L.s-1 in Ústí city and 44 L.s- 1in Dín city accompanied by maximal drop of groundwater temperature 1.7°C were evaluated by predictive simulations.

Key words

HUF MODFLOW, groundwater flow model, heat flow model, geological fault, terrestrial heat flow

1. Introduction

   Long term groundwater extraction has a significant effect on the groundwater levels, flow directions and velocity [14]. In the areas with thermal waters the groundwater extraction significantly influences the temperature distribution

   in the aquifer [1]. The most useful tool for evaluation of groundwater pumping in the complicated hydrogeology systems is calibrated mathematical model [15].

   In the Czech Republic one of largest and most exploited aquifer system with low potential thermal groundwater accumulation is the Benesov-Usti aquifer system (BUAS). It occupies N part of the Bohemian Cretaceous basin (BCB) and covers cca 1580 km2 (Fig. 1).

   Fig. 1 Area under study

   A mathematical model of groundwater flow and heat transport has been designed for the needs of ground water management in the area under study. The main goal of the research project was to determine available quantity of thermal groundwater in the areas of Decin and Usti and the effect of withdrawals on the temperature field. The

   mathematical model allowed for evaluation of temperature changes, as well as identification of infiltration zones.

   The majority of pumping wells is situated in the Ústní nad Labem and Dín cities. Thermal groundwater temperature in the area typically varies between 25 and 35°C. Estimates of the existing heat flow in the area differ quite significantly. In [10] the heat flow estimate varies between 90100 mW.m-2. An average heat flow in the BCB of 73 mW.m-2 and 80 mW.m-2 is assumed in the [8] and [3]. In [19] the heat flux values lower than 50 mWm-2 for the recharge and higher than 70 mWm-2 for discharge areas is determined. Many authors expect increased heat flow in the fault areas.

   In the winter period 50 ls-1 of thermal groundwater is used by heat exchangers to allow central heating in the Dín city. Thermal groundwater become heat source of thermal pools in Usti nad Labem and Decin and for Usti nad Labem ZOO.

   The area under study is strongly affected by tectonic deformation resulting in many faults. The groundwater may overflow between different aquifers at the fault lines [4]. The broadly used groundwater model vertical discretization (using MODFLOW BCF and LPF packages) fails to account for complex groundwater pattern existing in the described area.

2. Hydrogeology, model concept

   1. Surveying history of the area

      The first wells were drilled from 1888 to 1932 in the areas of Bystrany, Usti nad Labem and Decin (Fig. 1). The hydrological and hydrogeological research of the BCB was subject to an extensive assessment in the 1980s. A comprehensive description of the area was published by [7]. In [11] a detailed hydrological and hydrogeological evaluation of the region including double layer mathematical model is published.

      The European project ISPA (2002-2008, No 2000/CZ/16/P/PE/003 – reconstruction of national groundwater monitoring network) brought a new inputs to the research of the area of interest. Sixteen new wells were drilled in the model area. Extensive drill logging and surveying helped verify and precise data of lithology, basin tectonic structure, hydraulic parameters, piezometric altitudes and vertical profiles of groundwater temperatures. Synthesized information from all available sources is in [5].

   2. Fault zones

      The study area forms part of the territory of intense Saxonian orogeny, the predominant trends

      of faults being NE-SW, NW-SE, and E-W. Fig. 2 shows the outline of fault lines. The entire sunken Stredohori block is very strongly affected by block disintegration along faults of various trends. Due to this fact and the varied lithofacies development of the Upper Cretaceous sedimentation, this unit is one of the most complex structures of the BCB. The most important tectonic features include the Ore Mountains fault zone, the Lusatian fault and the complex of the Stredohori fault (Fig. 2).

      Fig. 2 Outline of geological faults: 1 Ore Mountains fault zone, 2 Stredohori fault zone (2a Straz fault , 2b Ceska Lipa fault zone, 2c Ustek fault, 2d Libochovice fault), 3 Malecov-Okresice fault, 4 Valkerice fault, 5 Radec fault, 6 Kerhartice fault, 7 Ceska Kamenice fault zone 8 Svor fault, 9 Decin fault zone , 10 Lusatian fault, 11 Libouchec fault, 12 Zitenice fault, 13 Skalice fault, 14 Bechlejovice fault, 15 Doubice fault zone

   3. System of aquifers and aquitards

      Lithological distinction of Cretaceous beds is the key factor to identify the aquifers and aquitards of the entire unit [11]. The sunken Stredohori block is the main structural element in the area.

      Three aquifers can be distinguished in the sunken Stredohori block. The base aquifer A is accommodated in the Peruc-Korycany formation (Cenomanian sandstones), aquifer BC in the Bila Hora and Jizera formations (Lower to Middle Turonian sandstones), and aquifer D is constituted by the Brezno and Merboltice formations (Coniacian sandstones). Aquitards a/bc (the base of the Bila Hora formation -marlstones) separates aqifers A and BC. Aquitards bc/d (Jizera and predominantly Teplice formations claystones, marlstones) separates aquifers BC and D. The basal Cretaceous aquifer A occurs throughout the model area. Aquitards a/bc did not develop in the area east of Usti nad Labem, where only joint aquifer AB is present. The thickness of aquifer A and/or AB is typically 40-70 m. The thickness increases to approximately 130 m in the vicinity of the Lusatian fault, where, due to the petering out of the aquitard beds, it connects also with the overlying aquifers to form a single permeable aquifer complex of sandstone beds up to 750 m thick. The lowest part of the Bila Hora formation in the rest of the area is constituted by marlstones. It functions as the upper aquitard for aquifer A. Transmissivity of the aquifer A decreases from NE to SW and W in

      agreement wth the decrease of the sandy component. The mean value of the coefficient of transmissivity is 60 m2/day.

      The main Cretaceous aquifer BC, constituted by the sandstones of the Jizera formation (C) and the upper part of the Bila Hora formation (B) is the thickest and, in terms of water management, the most important aquifer of the entire Bohemian Cretaceous Basin. The aquifer BC extends across most of the area, only missing in its SW part, roughly beyond the Decin-Ustek line (Fig. 2). The contiguous body of Jizera formation sandstones splits WSW of this line to several bodies and peters out. Aquifer BC is connected to the sandstones of aquifer A in the zone along the Lusatian fault (the aquifer ABC). Aquifer BC in the central and SW parts of the area is overlaid by the Teplice and Brezno formations in the facies of calcareous claystones and marlstones, which function as the topping of the aquifer (aquitards bc/d). The thickness of aquifer BC varies from 60 m in the SW to 510 in the NE at the Lusatian fault. The above thickness of aquifer BC is the greatest of the entire BCB.

      The uppermost part of the basin (aquifer D) is constituted by the sandstone facies of the Brezno formation and/or the flyschoid facies of the Merboltice formation. Aquifer D in the central part of the sunken Stredohori block is intruded (and partly covered) by extensive volcanic bodies, which complicate groundwater flow. The upper boundary of aquifer D is constituted by free water table [7].

   4. Vertical discretization of the model area Complex groundwater flow simulation (using MODFLOW) is only possible whit the vertical discretization of the model domain independently of the base and ceiling of existing litostratigraphy units. Given the main aquifer BC thickness of 450 m and approx. 30 m of geothermal gradient, groundwater temperature difference of 15°C may occur in the aquifer. The simulation of heat transport also requires a more detailed discretization of the model space. The scheme of the hydrogeological bodies and faults in the study area, groundwater flow and model layers concept is

      in the Fig. 3.

      Anderman and Hill (2000) developed MODFLOW module HUF, based on a new methodology of the input of hydraulic characteristics independently to the model grid. Hydraulic characteristics of the model cells are calculated from the characteristics described for existing hydrogeological units.

      Fig. 3 Concept of hydrogeological units (schematic cross-section NE-SW)

      Above described discretization concept distributes river boundary condition within several model layers. Uncontrolled drying of model cells during groundwater flow calculation followed. Improved vertical discretization, preserving all the stream cells in the first model layer, is described in Fig. 4. The model space was ultimately discretized into 21 model layers.

      Fig. 4. Schema of the hydrogeological units and model layers

3. Groundwater flow model

   1. Boundary conditions

      Rainfall recharge is the only source of groundwater in the model area. The infiltrated quantity was found by means of hydrological calculation, whereby rainfall infiltration is equal to the observed base flow from the territory. The quantity of recharge in individual model cells depends on the altitude of the surface and the rock type at outcrop. The quantity infiltrating into the entire model area (1580 km2) is 6,315 L.s-1.km2. Model recharge was calculated in the interval from 0.6 to 8.0 L.s-1.km-2.

      An important part of the infiltration is drained into the rivers. The drainage into river streams was simulated using the boundary condition of the third type.

      Groundwater abstraction was simulated using boundary condition of the second type. It is practically impossible to comply with the requirement that groundwater pressure along the open section of the well (including the case of connection of more aquifers) should be approximately constant in the multiaquifer system using MODFLOW module WELL. Groundwater abstraction was simulated using module MODFLOW MNW1 [6].

   2. Calibration of the groundwater flow model

      The groundwater flow model was calibrated applying the head criterion in three stationary simulations describing different exploitation of the system. Model values of hydraulic levels were compared with known water table data from the given period.

      The first simulation represents flow pattern existing prior to the beginning of groundwater withdrawal (in late 19th century).

      The second simulations involved maximum achieved exploitation of the BUAS system in 1980s (583.6 L.s-1). The model water level was compared with several tens of measurements of water tables in wells, listed by [11] (Fig. 5). Model water levels of artesian aquifers in the drainage areas are usually in agreement with real measurement to 5 m. bigger differences were allowed in the recharge areas (Fig. 5).

      Fig. 5 Comparison of model deviations from the measured groundwater levels

      Model drainage in the individual catchments was compared with the evaluated base flows in the river network (Fig. 6).

      Fig. 6 Comparison of the measured and modelled groundwater flows

      The third option involved simulations of groundwater flows at the existing mean rates of withdrawal (280 L.s-1).

   3. Results of the groundwater flow model

      The model results reveal that, given the tectonic affected structure, the piezometric values of groundwater flow in the basal aquifer A and the main aquifer BC are at similar levels. Vertical shift at the fault lines is usually greater than the thickness of the basal aquifer A. The block structure and the complicated system of groundwater overflow between the main and basal aquifers do not allow the existence of a piezometric heads in aquifer A independent to main aquifer.

      An undisturbed stream of groundwater from the area of Usti nad Labem (in 19th century) trended toward the western boundary of the model territory, into the area of the river Bilina near Teplice. Existing withdrawal in the area of Usti nad Labem (commenced in 1911) created a new site of artificial drainage in the region (Fig. 7), substantially altering the trends of groundwater flow in an area of more than 100 km2. Withdrawals of groundwater in the quantity of 147 L.s-1 and 27 L.s-1 in the areas of Decin and Usti nad Labem in the period of maximum exploitation caused a drop of the pressure head of the artesian aquifers by dozens of metres (Fig. 8). Withdrawal of 238 L.s-1 in the area of Ceska Lipa, where aquifer BC is unconfined, caused a drop of water table up to only 8 m.

      Fig. 7 Model contours of groundwater levels in aquifer BC at maximum withdrawals in 1980s

      Fig. 8 Drop of groundwater levels at maximum withdrawals in 1980s

4. Heat flow model

   1. Methodology and goals

      Goals of the model were: a/find and test suitable methodology for the heat flow modelling in complicated hydrogeological structures, b/ evaluate effect of groundwater flow on temperature field and c/ evaluate the long term temperature changes due to the groundwater pumping. Thermal characteristics of the model were set an accordance with methodology publicised by [12].

      The influence of temperature changes on the hydraulic and thermal characteristics is omitted.

      The temperature changes in the structure are small. According to [12] heat calculation was performed out using MT3DMS software [16].

      The computed temperature field for unaffected condition by groundwater abstraction was set as initial condition for all others simulations. The effect of groundwater flow and groundwater pumping on temperature field was determined from the difference between the temperature field computed with and without groundwater flow. The temperature field in the saturated rock formation without groundwater flow was computed by setting thousand times smaller values of hydraulic conductivity. and rain infiltration into the groundwater and heat flow model.

      Fig. 9 Model contours of average temperatures in the collector BC – conditions without pumping

      Fig. 10 Groundwater temperature change in collector A caused by groundwater flow

   2. Boundary conditions and heat transport parameters

      The basal heat inflow is simulated using second type boundary condition between 70 – 80 mW.m-2 in the lowest model layer (bedrock of the cretaceous sediments).

      The heat outflow from the model is simulated using first type boundary condition. The constant temperature is set to the first model layer. Part of the heat is drainage with groundwater to the streams, springs and pumping wells.

      There is lack of information describing thermal properties of the saturated rock.

      The thermal parameters of saturated rock environment was set up using methodology publicised in [12]. The thermal conductivity 1.66 W.m-1.°C-1 and thermal capacity 920 J.kg-1.°C were set up in the model.

   3. Heat flow model results

      Model calculated temperatures were compared with temperatures measured usually at the borehole bottom (fig. 11).

      Good agreement of measured and modelled temperature (the difference smaller than 3°C) was reached in the area of Dín and Ústí nad Labem cities. Modelled temperatures lower than 10°C in compare to measurement were computed in the central part of the model area. Simulated downward groundwater flow from the aquifer BC to the base

      aquifer A is probably of higher intensity compared to existing conditions.

      45

      Simulated temperature (°C)

      Simulated temperature (°C)

      40 Temperature

      Best fits

      35

      30

      25

      20

      15

      10

      5

      5 10 15 20 25 30 35 40 45

      Observed temperature (°C)

      Fig. 11 Comparison of the measured and modelled temperatures

      The areal distribution of the average vertical temperatures for the collector BC is shown in the Fig. 9. Lowest model temperatures are in the shallow parts of the hydrogeology structure and in the infiltration areas especially. Highest temperatures are measured and modelled in the areas with slow groundwater flow (Ústí nad Labem area) and in the deepest parts of the structure (close to Beneov nad Plounicí city).

      The groundwater flow has a fundamental effect on the temperature distribution in the model area (Fig. 10). In the area between Ústí and Teplice cities the temperature differences between affected and unaffected temperatures via groundwater flow are relatively small. But in the infiltration areas the groundwater temperature in the basal collector A would be more than 20 °C higher without groundwater flow.

      Fig. 12 Drop of the groundwater temperature in pumped wells in the Dín City area

      The long term development of the model temperature in the pumped wells in the Dín city is in the fig. 12. The temperature decrease is very

      small and its progress much beyond the human lifespan. This model results corresponds to reality. Its known that in the most exploited area the temperature decrease was not observed (after almost hundred years of groundwater pumping).

5. Conclusions

Fault lines fragment the model area of the Benesov-Usti system into 30 separate blocks with different altitudes of the base of Cretaceous sediments (Fig. 2). The vertical shift at individual faults causes discontinuities of the basal aquifer A. Horizontal overflows of groundwater between various aquifers occurs across fault lines.

A system of multiple aquifers in a block structure as occurs in the model area does not allow vertical discretization of the model domain in accordance with the geometry of existing hydrogeological units. Such approach would totally neglect the effect of faults on the groundwater flow pattern. Discretization of the model layers independently to the geological setup seems to be adequate approach. The application of the MODFLOW module HUF for the calculation of groundwater flow allowed automatic input and effective calibration of the hydraulic characteristics of the model space and should be used for groundwater flow simulations in complicated hydrogeology structures.

Potential to increase pumping of thermal groundwater in the areas of Usti nad Labem and Decin from 80.8 L.s-1 and 25.3 L.s-1 to 275 L.s-1 and

44 L.s-1 was evaluated by predictive model simulations. Retention of the piezometric surface of artesian aquifers above the ground surface was the applied criterion in the specification of maximum groundwater withdrawal. The water temperature stays at the same level for at least several generations, as subsequently confirmed the predictive heat transport model.

The computed temperature field is in reasonable agreement with measured temperatures. Model results were achieved by setting almost uniform heat flow values to the basal model layer. Model results suggest that heat flow in the area under study could be quite homogenous. According to model results the areas with higher temperatures are more influenced by upward flow of heated groundwater than by higher natural heat flow from the bedrock.

Results of the groundwater and heat flow model could be supposed new input to structural geologists discussing the function of faults in the area. It is known that the vertical shift of individual blocks is not attained in a single fault line but typically in a zone of several parallel faults. The

tightness of aquitards separating the aquifers can thus be significantly disturbed.

Demands have been voiced for an updated and detailed reassessment of the tectonic structure of the entire area of interest. Temperature field as well as geochemistry should be taken into account into hydrogeological analysis of the fault zones.

Acknowledgements

This research has been supported by the Czech Science Foundation GACR/205/07/0691 and by the Institutional research project of the Czech Ministry for Environment MZP/00020711010021620855 Research and protection of hydrosphere.

References

1. Andrews, Ch.B. (1978). The impact of the use of heat pumps on ground-water temperatures.

   Ground water, Vol. 16, No.6, 437 – 443

2. Anderman, E.R.; Hill, M.C., (2000). MODFLOW-2000, the U.S. Geological Survey modular ground-water model, Documentation of the Hydrogeologic-Unit Flow (HUF) Package: U.S. Geological Survey Open-File Report 00- 342,

   Reston, Virginia

3. Cermak, V. (1977). Geothermal measurements in Paleogene, Cretacecous and Permo-Carboniferous sediments in Northern Bohemia. Geophys. Roy.Astron.Soc., 48:537-541,

   Oxford

4. Datel, J. V.; Kobr, M.; Prochazka, M., (2008). Well logging methods in groundwater's surveys of complicated aquifer systems: Bohemian Cretaceous Basin. Environ.Geol., 57:1021-1034, doi: 10.1007/s00254-008-1388-8

5. Datel, J. V., (2008). Limits for use of thermal waters on the example of the conceptual model of the Benesov-Usti aquifer system of the Bohemian Cretaceous Basin. Dissertations, Charles University in Prague

6. Halford, K. J.; Hanson, R.T., (2002). User guide for the drawdown-limited, multi-node well (MNW) package for the U.S. Geological Survey's modular three-dimensional finite-difference ground-water flow model, versions MODFLOW- 96 and MODFLOW-2000: U.S. Geological Survey Open-File Report 02-293, Reston, Virginia

7. Hercik, F.; Herrmann, Z.; Valecka, J., (1999). Hydrogeology of the Bohemian Cretaceous Basin. Czech Geological Survey, Prague

8. Ibrmajer, J.; Suk, M.; et. al., (1989). Geophysical view of the Czech Republic. Czech geological survey, Prague

9. Jiráková, H.; Procházka, M.; Ddeek, P.; et al. (2011). Geothermal assessment of the deep aquifers of the north-western part of the Bohemian

   Cretaceous basin, Czech Republic. Geothermics 40, 112-124.

10. Myslil, V.; Kukal, Z.; Posmourny, K.; Frydrych, V., (2007). Geothermal energy ecological energy from depths of the Earth. J.

    Planeta, 15:1-32, Praue

11. Nakladal, V. et al., (1987). Groundwater resources and hydrogeological conditions of the Bohemian Cretaceous Basin. Hydrological unit 3.

    Stavebni Geologie, Prague

12. Thorne, D.; Langevin, Ch.; Sukop, C. M., (2006). Addition of simultaneous heat and solute transport and variable fluid viscosity to SEAWAT, Computers & Geosciences 32,1758 1768,

1. Tularam, G. A.; Krishna, M., (2009). Long term consequences of groundwater pumping in Australia: A review of impacts around the globe. Journal of Applied Sciences in Environmental Sanitation, Vol 4, Num. 2: 151-166.

2. Williams, D. D.; et al., (1999). Analysis of Convective Heat Transfer in Deformed and stratified aquifers associated with frasch thermal mining, Groundwater 37, 517 – 522

3. Zheng, C.; Wang, P. P., (1999). MT3DMS, a modular three-dimensional multi-species transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; documentation and users guide SERDP-99-1. U.S. Army Engineer Research and Development Centre, Vicksburg.