Calibration and Reliability in Groundwater Modelling (Proceedings of the ModelCARE 99 Conference held at Zurich, Switzerland, September 1999). IAHS Publ. no. 265, 2000.

Characterization of gypsum aquifers using a
coupled continuum-pipe flow model

STEFFEN BIRK, RUDOLF LIEDL
Applied Geology, Geological Institute, University of Tubingen, Sigwartstrasse 10,
D-72076 Tubingen, Germany

e-mail: steffen.birk@uni-tuebingen.de

MARTIN SAUTER
Institute of Geosciences, University of Jena, Burgweg 11, D-07749 Jena, Germany

Abstract The growth of a single conduit in a gypsum layer located between two confined aquifers is studied using a continuum-pipe flow model (EVE).
The initial diameter of the conduit and the exchange of water between the conduit and the gypsum layer are controlling factors at the initiation stage of conduit development. However they are unimportant in the long-term evolution of the conduit. Controlling factors at the later stages of conduit enlargement are identified using an analytical solution.


INTRODUCTION

Groundwater circulation in gypsiferous horizons and subsequent erosion is frequently the reason for geomechanical problems such as subsidence or even collapse. In order to provide a basis for reliable risk assessment, an adequate hydrogeological characterization of gypsum karst aquifers is required, mainly focusing on the geometric and hydraulic properties of the highly conductive conduit network. However this is not easily achieved by standard investigation techniques such as hydraulic tests. Therefore, Liedl & Sauter (1998) suggested an integrated modelling approach combining both long-term simulations of aquifer genesis and short-term simulations of transport phenomena in karst aquifers.
This paper presents a new modelling tool termed EVE (Evaporite Void Evolution). EVE is designed to support the characterization of gypsum aquifers which requires the numerical forward modelling of conduit genesis. The evolution of a single conduit in a gypsum layer is studied in order to identify the controlling parameters for conduit enlargement.
 

MODELLING APPROACH

A karst groundwater flow system may be conceptualized as a flow system consisting of a so-called fissured system, which represents the mass of the fractured rock, and a conduit system representing the karst pipe network. Groundwater flow in the fissured system is modelled using Boussinesq's equation. Flow in the conduit system is governed by Kirchhoff s law that states that the total inflow and total outflow balance Characterization of gypsum aquifers using a coupled contiuum-pipe flow model 17 for each node of the conduit network. For each pipe the discharge is related to the head difference by the Darcy-Weisbach equation which is adapted to laminar and turbulent flow conditions. Exchange of groundwater between the conduits and the fissured system is assumed to be proportional to the hydraulic head difference between these flow systems. Mass transport in the conduits is described by the one-dimensional (1-D) advection equation with an additional source term accounting for the increase in concentration due to the dissolution reaction.

Clemens et al. (1996) implemented this modelling approach for carbonate aquifers in the CAVE (Carbonate Aquifer Void Evolution) code, which added conduit flow and limestone dissolution to the well-known MODFLOW code. To simulate the evolution of voids in gypsum aquifers for a variety of boundary conditions, EVE combines the updated MODFLOW-96 code (McDonald & Harbaugh, 1996) with the conduit flow module from CAVE and a newly developed evaporite dissolution module.
In order to calculate the rate of conduit growth in a gypsum aquifer the dissolution rate of gypsum is required. Since most studies of gypsum dissolution kinetics agree that the diffusion from the solid gypsum phase into the mobile conduit water is the controlling (i.e. limiting) step of the dissolution process (Klimchouk, 1996a), the dissolution rate can be described by a first-order rate law (James & Lupton, 1978):
 
(1)

where δM/δt is the mass loss from the solid surface with the surface area A, h is the mass transfer coefficient, ceq is the equilibrium concentration, and c is the concentration in the bulk solution.
The mass transfer coefficient is related to the diameter of the conduits d and the diffusion coefficient D via the dimensionless Sherwood number Sh:
 
(2)

For laminar flow conditions, a fully developed diffusion boundary layer is assumed, i.e. Sh = 3.66 in circular pipes with a parabolic velocity profile, whereas the Sherwood number for turbulent flow conditions is given by an empirical correlation (Incropera & DeWitt, 1996).
 

MODEL SCENARIO

It is generally believed that the evolution of gypsum karst starts under artesian conditions (Klimchouk, 1996b), a typical setting consisting of two confined aquifers separated by a gypsum layer. The model domain studied in this paper consists of a two-dimensional (2-D) cross-section with a total thickness of 70 m subdivided into a gypsum layer of 22 m and two aquifers of 24 m thickness (Fig. 1). The general head boundaries represent a recharge area over a distance of 2 km at a level of 150 m, and a discharge area at a distance of 2 km at a level of 70 m. In addition, a river drains the upper aquifer with a leakage factor of 2 x 10-9 m-1 and a water level of 75 m. This setup imposes a hydraulic head gradient between the aquifers and causes an upward directed flow from the lower to the upper aquifer. The hydraulic conductivity is set to 10-6 for the aquifers and l0-9 m s-1 for the gypsum layer. A single conduit of 0.4 mm diameter intersects the gypsum layer. At four nodes, which subdivide the conduit into three pipes, an exchange of water between the conduit and the fissured system is allowed. Since flow is always from the lower aquifer into the gypsum layer, the gypsum concentration of water flowing in the lower aquifer is assumed to be zero. The water flowing in the gypsum layer is assumed to be saturated with respect to gypsum.


Fig. 1 Conceptual setting and model domain with initial hydraulic heads.


SIMULATION RESULTS

Regarding pipe diameters and hydraulic heads, three stages of conduit development can be distinguished (Fig. 2). At the initiation stage the water flowing from the conduit into the upper aquifer is saturated with respect to gypsum. Therefore, the outlet of the conduit is not enlarged and restricts the discharge. Due to the sluggish flow conditions the conduit growth is very slow and the hydraulic head gradient between the aquifers is maintained (Fig. 2(a)). However, the inlet of the conduit is enlarged by the strongly undersaturated water entering the conduit. Thus, the hydraulic head of the lower aquifer propagates upward through the conduit, and the hydraulic head gradient along the upper part of the conduit (i.e. pipe 3) increases accordingly (Fig. 2(b)). As a result the discharge through the conduit increases, and eventually the water emerging at the outlet of the conduit is undersaturated with respect to gypsum. In this situation a positive feedback mechanism triggers a rapid conduit growth (breakthrough), since the Characterization of gypsum aquifers using a coupled contiuum-pipe flow model enlargement of the outlet due to gypsum dissolution results in increasing flow rates and decreasing concentrations and thus higher dissolution rates. Finally the hydraulic head gradient between the aquifers is diminished (Fig. 2(c)).
After the breakthrough has occurred the hydraulic resistance of the conduit is less than those of the aquifers. Therefore, the discharge through the conduit is controlled by the boundary conditions and the hydraulic conductivity of the aquifers, but no longer by the diameter of the conduit. Thus, at the main stage of conduit development the flow rate remains constant (Fig. 3). Whereas the exchange of water between the conduit and the gypsum layer is negligible at this late stage, it is an important factor in conduit genesis at the initiation stage. The flow rates increase in the lower part (pipe 1) of the conduit, although the outflow into the upper aquifer is restricted by the small diameter of the outlet (pipe 3). Thus, the conduit is drained by the gypsum layer and supplied with additional aggressive water from the lower aquifer. In a model run with no exchange between the conduit and the gypsum layer, breakthrough did not occur within 1000 years. A similar result showing the importance of the exchange parameter was obtained by Bauer et al. (1999) when they studied the karstification of carbonate rocks beneath a dam.


Fig. 2 Pipe diameters and hydraulic heads at different stages.
 


Fig. 3 Flow rates in the pipes.

A further parameter, which is expected to be a controlling factor in conduit genesis, is the initial pipe diameter d0. Therefore, several model runs starting with different initial diameters were performed (Fig. 4). Although the breakthrough time is very sensitive to the initial diameter, the final diameter of the conduit (pipe 1) is nearly the same after 1000 years in each case. Thus, considering the long-term evolution, i.e. a period which is long relative to the breakthrough time, conduit development is insensitive to the initial diameter. Moreover, this long-term evolution can be described analytically by inserting equation (2) in equation (1) and setting δM(t) = ½pAδ(t), where p is the density of gypsum, and d(t) is the pipe diameter at time t. The resulting relationship can be solved for the pipe diameter yielding:
 
(3)

Figure 4 shows the temporal evolution of the diameter of pipe 1 comparing EVE results with the analytical solution given by equation (3). Since for the given model Characterization of gypsum aquifers using a coupled contiuum-pipe flow model 21 scenario the gypsum concentration after the breakthrough is very small compared with the equilibrium concentration, a concentration c = 0 was inserted in equation (3), and the function was plotted for an initial diameter of 1 mm. Smaller initial diameters result in the same straight line as for do = 1 mm, whilst the plot of equation (3) for an initial diameter of 1 cm yields the same result as calculated by EVE.


Fig. 4 Diameters of pipe 1 for various initial diameters calculated by EVE, and comparison with diameters from equation (3) for an initial diameter of 1 mm.

Equation (3) demonstrates that the influence of the initial diameter is negligible for long time periods. Conduit enlargement at the main stage is only controlled by the density of the rock, by the undersaturation of the water with respect to gypsum, and by parameters which describe the diffusion process of gypsum from the solid surface into the bulk solution (Sh and D). Note however that the concentration of gypsum in the water depends on the flow rate.
 

CONCLUSIONS

At the initiation stage, conduit development under artesian conditions is controlled by the initial diameter of the conduit and by the exchange of water between the conduit and the gypsum layer. However, at the late stage conduit growth is insensitive to these parameters. The long-term evolution of conduits with different initial diameters leads to similar diameters, which can be predicted by an analytical solution.
Further investigations will have to address the question as to whether the relative importance of other parameters (e.g. the hydraulic gradient along the conduit), on conduit development is superior to the impact of the parameters studied here.

Acknowledgements This work is part of the ROSES (Risk of Subsidence due to Evaporite Solution) project funded by the European Commission Framework IV Programme (Contract number ENV4-CT97-0603).
 

REFERENCES
 
Bauer, S., Birk, S., Liedl, R. & Sauter, M. (1999) Solutionally enhanced leakage rates of dams in karst regions. In: Karst Modeling (ed. by A. N. Palmer, M. V. Palmer & I. D. Sasowsky), 158-162. Karst Waters Institute Special Publ. no. 5.
Clemens, T., Hùckinghaus, D., Sauter, M., Liedl, R. & Teutsch, G. (1996) A combined continuum and discrete network reactive transport model for the simulation of karst development. In: Calibration and Reliability in Groundwater Modelling (ed. by K. Kovar & P. van der Heijde) (Proc. ModelCARE 96 Conf., Golden, Colorado, September 1996), 309-318. IAHS Publ.no. 237.
Incropera, F. P. & DeWitt, D. P. (1996) Fundamentals of Heat and Mass Transfer. Wiley, New York.
James, A. N. & Lupton, A. R. R. (1978) Gypsum and anhydrite in foundations of hydraulic structures. Geotechnique 28, 249-272.
Klimchouk, A. (1996a) The dissolution and conversion of gypsum and anhydrite. Int. J. Speleol. 25(3-4), 21-36.
Klimchouk, A. (1996b) Speleogenesis in gypsum. Int. J. Speleol. 25(3-4), 61-82.
Liedl, R. & Sauter, M. (1998) Modelling of aquifer genesis and heat transport in karst systems. Bulletin d'Hydrogéologie 16, 185-200.
McDonald, M. G. & Harbaugh, A. W. (1996) Programmer's documentation for M0DFL0W-96, an update to the US Geological Survey modular finite-difference ground-water flow model. USGS Open File Report 96-486.