employing distributive modeling approaches
4 Groundwater residence time distributions in heterogeneous karst aquifers
4.2 Methods and approach
4.2.2 Model scenarios and parameterization
Since this study focuses on groundwater residence times within the aquifer, i.e. the saturated zone, the influences of the overburden and unsaturated zone are neglected in the model setup and the aquifer is simulated as confined. The effect this simplification might have for an actual field site application is discussed in Chapter 4.4. Two basic model configurations are used for studying the residence time distributions (Figure 4.1). Both configurations are three-‐dimensional with a lateral extent of 5×5 km2 and an aquifer thickness of 100 m. An observation well is inserted south of the conduit system at a depth of 50 m to simulate sampling in the field (Figure 4.1). The conduit volume is identical for both configurations, but in configuration 1 it is distributed on a 3 km long single conduit, while for configuration 2 a dendritic conduit system with a total length of 16 km is employed. The conduit diameter for configuration 1 is kept constant along the conduit length. In configuration 2 a widening of the karst conduits towards the spring is assumed. With a constant diameter conduit velocities would increase drastically at the intersections, which is not only unrealistic since higher flow velocities enhance karst dissolution processes (e.g. Clemens, 1998), but also numerically difficult to solve in the transport simulation. For the increase in radius, the empirical approach of Oehlmann et al. (2015) was employed (Eq. 4.6).
𝑟!,!=𝑚 𝑠+𝑏 (4.6a) 𝑟!,! =𝑚 𝑠+2 𝑟!,!!!!+𝑟!,!!!!+𝑟!,!!!! (4.6b) where rc [L] is the conduit radius, s [L] is the conduit length and m [-‐] is the slope of linear radius increase along the conduit length. For the smallest conduit branches an initial radius b [L] is defined (Eq. 4.6a). At the conduit intersections, the cross-‐sectional areas of the intersecting conduits are added as initial cross-‐section for the downgradient conduit (Eq. 4.6b). The factor 2 was derived empirically to ensure that flow velocities are as uniform as possible during the simulations. The intersections of three conduit branches at the same points lead to a significant increase in flow velocities, if the factor of 1 (Oehlmann et al., 2015) is used.
For selection of parameter values and ranges for sensitivity analysis, the area of the Gallusquelle spring in south-‐western Germany was used. The Gallusquelle is a medium sized karst spring with an average annual discharge of 0.5 m3 d-‐1. The aquifer is characterized as a mixed system, where both, conduit flow and diffuse matrix flow occur and are of significant importance (Sauter, 1992). Extensive field investigations and model studies provide a good database for model parameters (e.g. Sauter, 1992; Oehlmann et al., 2015). Average groundwater transit times in the saturated zone were determined by Geyer (2008) with the 85Kr method to range between 3 and 4 years. Table 4.1 shows the chosen parameters for the reference simulations and the variation ranges. The relative
parameter sensitivity was calculated by using the Root Square Error (RSE) of the average age and life expectancy values with respect to the reference scenario. For a relative difference below 0.5 years the parameter was considered as insensitive with respect to the objective function, i.e. age or life expectancy. For the spring discharge the same value was set as by Oehlmann et al. (2015) and a difference of 10 L s-‐1 was counted as significant.
For the porous system, all boundaries are zero-‐flux Neumann boundaries. In and out flux for the system is only provided by exchange with the fissured system. The fissured system is bounded by no-‐
flux boundaries everywhere except for the top. The whole upper boundary is defined as a Neumann boundary with the value of groundwater recharge as defined in Table 4.1. For the transport equation of groundwater age (Eq. 4.3), top of the domain is a zero-‐flux Neumann condition as well. By definition, groundwater age cannot enter the aquifer from the outside but is only produced inside of it. For the life expectancy the sign of the recharge is reversed and multiplied with the expectancy value to remove the water that reached the inlet boundary. The upper eastern edge of the domain is set as a Dirichlet boundary condition for groundwater flow and a Neumann condition for transport representing a river (Figure 4.1). The conduit system has the same kind of boundary condition for the karst spring. Recharge to the conduit is provided by exchange with the fissured system and by a source term representing the direct recharge component, i.e. recharge reaching the conduit system directly through vertical shafts, if present. For both, river and spring, Neumann conditions are zero flux for the life expectancy and equal to the groundwater discharge multiplied by the age for groundwater age, analogous to the recharge boundary.
Table 4.1. Parameters for the numerical simulation and variation range for the parameter analysis. The corresponding equations are given in Chapter 4.2. The z-‐axis points upwards.
Parameter description Parameter
name Reference
setup 1 Reference
setup 2 Variation range Porous system
porosity θp (%) 1 1 1 – 10
porous–fissured exchange coefficient β (s-‐1) 3.3×10-‐12 3.3×10-‐12 1×10-‐14 – 1×10-‐8 Fissured system
total recharge r (mm d-‐1) 1.5 1.5 0.5 – 50
porosity θf (%) 1 1 1 – 10
hydraulic conductivity Kf (m s-‐1) 5×10-‐5 5×10-‐5 1×10-‐6 – 1×10-‐3
dispersivity εf (m) 50 50 5 -‐ 100
aquifer thickness maf (m) 100 100 10 – 100
Conduit system
vertical position zc (m) 100 100 0 – 100
direct recharge rd (%) 0 0 0 – 95
conduit-‐matrix exchange coefficient α (m2 s-‐1) 6.1×10-‐4 6.1×10-‐4 1×10-‐8 – 1×10-‐3
cross-‐section A (m2) 12 2.24a 0.6 -‐ 35
roughness n (s m-‐1/3) 3 3 0.01 -‐ 20
dispersivity εc (m) 7 7 2 -‐ 50
initial radius b (m) 1.95 0.1 0.001 – 2
radius increase m (-‐) 0 1×10-‐4 0 – 1×10-‐3
aaverage value
The reference simulations and parameter analysis were performed for steady-‐state conditions.
Steady-‐state conditions are useful for protection zone delineation and required as initial values for transient modelling. In order to assess the influence of discrete groundwater recharge events on the residence time distribution, an additional simulation was performed introducing a hypothetical recharge event with the duration of one week. Contrary to the steady-‐state reference simulation (Table 4.1) a direct recharge component of 10% was assumed for the transient simulation to include the influence of the duality of aquifer recharge. Furthermore, the dispersivity of the conduit system was set to 50 m. It was found to be insensitive during the parameter analysis (Chapter 4.3.2) and a
higher dispersion coefficient leads to a higher numerical stability of the fast transport in the karst conduit system.