Combining implicit geological modeling, field surveys, and hydrogeological modeling to describe groundwater flow in a karst aquifer

Abstract

In three-dimensional (3-D) implicit geological modeling, the bounding surfaces between geological units are automatically constructed from lithological contact data (position and orientation) and the location and orientation of potential faults. This approach was applied to conceptualize a karst aquifer in the Middle Triassic Muschelkalk Formation in southwest Germany, using digital elevation data, geological maps, borehole logs, and geological interpretation. Dip and strike measurements as well as soil-gas surveys of mantel-borne CO₂ were conducted to verify the existence of an unmapped fault. Implicit geological modeling allowed the straightforward assessment of the geological framework and rapid updates with incoming data. Simultaneous 3-D visualizations of the sedimentary units, tectonic features, hydraulic heads, and tracer tests provided insights into the karst-system hydraulics and helped guide the formulation of the conceptual hydrogeological model. The 3-D geological model was automatically translated into a numerical single-continuum steady-state groundwater model that was calibrated to match measured hydraulic heads, spring discharge rates, and flow directions observed in tracer tests. This was possible only by introducing discrete karst conduits, which were implemented as high-conductivity features in the numerical model. The numerical groundwater flow model was applied to initially assess the risk from limestone quarrying to local water supply wells with the help of particle tracking.

Résumé

Dans un modèle géologique tridimensionnel (3-D) implicite, les surfaces limitant les unités géologiques sont construites automatiquement à partir des données de contact lithologique (position et orientation) et de la localisation et l’orientation de failles potentielles. Cette approche a été appliquée pour conceptualiser un aquifère karstique dans la formation du Muschelkalk (Trias moyen) dans le sud-ouest de l’Allemagne, en utilisant un modèle numérique de terrain, des cartes géologiques, des logs de forage et une interprétation géologique. Des mesures de pendages ainsi que des suivis de gaz du sol du CO₂ d’origine mantellique ont été réalisés pour vérifier l’existence d’une faille non cartographiée. La modélisation géologique implicite a permis une évaluation directe du cadre géologique et des mises à jour rapides avec les nouvelles données. Les visualisations simultanées en 3-D des unités sédimentaires, des structures tectoniques, des charges hydrauliques et des essais de traçage ont fourni des informations détaillées sur l’hydraulique du système karstique et ont aidé à guider l’élaboration du modèle hydrogéologique conceptuel. Le modèle géologique 3-D a été automatiquement transféré dans un modèle numérique d’écoulement à continuum unique en régime permanent qui a été calibré pour reproduire les charges hydrauliques mesurées, les débits des sources et les directions d’écoulements indiquées par les essais de traçage. La calibration a été rendue possible seulement après l’introduction de conduits karstiques distincts, qui ont été implémentés comme des éléments de conductivité hydraulique élevée dans le modèle numérique. Le modèle numérique d’écoulements a été utilisé avec la méthode « suivi des particules » pour une évaluation initiale des risques liés à l’exploitation de calcaire en carrières sur les forages d’approvisionnement en eau potable à l’échelle locale.

Resumen

En la modelización geológica implícita tridimensional (3-D), las superficies que delimitan las unidades geológicas se construyen automáticamente a partir de los datos de contacto litológico (posición y orientación) y la ubicación y orientación de las posibles fallas. Este enfoque se aplicó para conceptualizar un acuífero kárstico en la formación Muschelkalk del Triásico Medio en el sudoeste de Alemania, utilizando datos digitales de elevación, mapas geológicos, registros de perforaciones e interpretación geológica. Se llevaron a cabo mediciones de buzamiento y rumbo, así como estudios del gas del suelo del CO₂ transportado por el manto para verificar la existencia de una falla no cartografiada. La modelización geológica implícita permitió la evaluación directa del marco geológico y una rápida actualización con los datos de entrada. Las visualizaciones tridimensionales simultáneas de las unidades sedimentarias, las características tectónicas, las cargas hidráulicas y las pruebas de trazadores permitieron comprender la hidráulica del sistema kárstico y ayudaron a orientar la formulación del modelo hidrogeológico conceptual. El modelo geológico tridimensional se tradujo automáticamente en un modelo numérico de aguas subterráneas en estado estacionario de un solo continuo que se calibró para que coincidiera con las cargas hidráulicas medidas, las tasas de descarga de los manantiales y las direcciones de flujo observadas en las pruebas de trazadores. Esto sólo fue posible mediante la introducción de conductos kársticos discretos, que se implementaron como características de alta conductividad en el modelo numérico. El modelo numérico de flujo de aguas subterráneas se aplicó para evaluar inicialmente el riesgo de la explotación de canteras de caliza en los pozos de abastecimiento de agua potable con la ayuda del trazado de partículas.

摘要

在三维(3-D)隐式地质建模中, 地质体之间的边界面是根据岩层接触数据(位置和方向)以及潜在断层的位置和方向自动构建的。采用数字高程数据, 地质图, 钻孔记录和地质解释方法, 将地质建模方法用于概化德国西南部中三叠世Muschelkalk地层的岩溶含水层。进行了倾角和走向测量以及壁炉架传播CO₂的土壤气调查, 以验证是否存在未发现的断层。隐式地质建模可以对地质结构进行直接评估,并可以快速更新输入数据。沉积单元,构造特征, 水头和示踪试验的3-D可视化同时提供了对岩溶系统水力学的认识, 并有助于建立概念性水文地质模型。将3-D地质模型自动转换为单连续体的稳定地下水数值模型, 并对其进行校准以匹配在示踪剂试验中观测水头, 泉流量和流向。仅通过引入离散的岩溶管道是可能的, 该管流在数值模型中可视为高渗透系数。借助粒子追踪技术, 地下水流数值模型被用于初步评估从灰岩采石场到当地供水井的风险。

Resumo

Na modelagem geológica tridimensional (3-D) implícita, as superfícies delimitadoras entre as unidades geológicas são construídas automaticamente a partir de dados de contato litológico (posição e orientação) e a localização e orientação de potenciais falhas. Essa abordagem foi aplicada para conceitualizar um aquífero cársico na formação do Médio Triássico Muschelkalk no sudoeste da Alemanha, usando dados digitais de elevação, mapas geológicos, relatórios de sondagem e interpretação geológica. Foram realizadas medições de mergulho e direção, bem como investigações no solo de gases CO₂ derivados do manto para verificar a existência de uma falha não mapeada. A modelagem geológica implícita permitiu a avaliação direta da geometria das unidades geológicas e atualizações rápidas com o recebimento de novos dados. As visualizações simultâneas 3-D das unidades sedimentares, contexto tectônico, cargas hidráulicas e ensaios com traçadores forneceram informações sobre a hidráulica do sistema cárstico e ajudaram a orientar a elaboração do modelo hidrogeológico conceitual. O modelo geológico 3-D foi automaticamente convertido em um modelo numérico de fluxo de meio contínuo equivalente em estado estacionário, que foi calibrado para corresponder às cargas hidráulicas observadas, descargas das nascentes e às direções de fluxo observadas nos ensaios com traçadores. Isso foi possível apenas com a introdução de feições discretas de condutos cársticos, que foram implementadas como feições de alta condutividade no modelo numérico. O modelo numérico de fluxo de águas subterrâneas foi aplicado para avaliar inicialmente o risco da extração de calcário para os poços locais de abastecimento de água com a ajuda do rastreamento de partículas.

Introduction

Understanding the geology is the backbone of any groundwater management project, since hydraulic heads and velocities, as well as flow trajectories and residence times, are strongly controlled by the properties and arrangement of lithological entities. The hydraulic connection between aquifers depends on their shape, arrangement, and internal structure (Raiber et al. 2015). Regional flow systems can be segmented by faults, which may form either barriers or preferential flow paths (Caine et al. 1996; Moya et al. 2014; Raiber et al. 2015). Groundwater models that are merely developed on hydraulic interpretations, while disregarding geological insights, are likely to be erroneous or untrustworthy (Fogg 1986).

Shaped by dissolution, internal erosion and subsidence processes, karst systems are commonly characterized by sinkholes, caves, springs, and conduits forming underground drainage systems (Ghasemizadeh et al. 2012). Although karst can develop from various evaporites, only limestone and dolostone karsts are relevant for water management purposes (Hartmann et al. 2014). Describing karst aquifers for the purpose of quantifying groundwater flow poses unique challenges due to the complex and highly uncertain arrangement of karst conduits where water flow often occurs under turbulent conditions at very high velocities (e.g., Jeannin 2001). Karst conduits comprise only a minor fraction of the total aquifer porosity but significantly control the aquifer hydraulics (e.g., Ghasemizadeh et al. 2012; Kuniansky 2016). Understanding the structure and hydraulic significance of fractures and conduits is therefore fundamental for the description of groundwater flow in karst aquifers. While the exact locations and properties of karst conduits are difficult to predict, their general orientation is strongly dominated by basic geological features such as the orientations of bedding planes, joints and faults, as well as the type and extent of the overburden, which strongly controls the karstification process (Ford 2003). To obtain an appropriate hydrogeological conceptual model of a karst aquifer, information on structural geology necessarily needs to be complemented by field observations and surveys such as tracer tests (e.g., Perrin and Luetscher 2008). This conceptualization is a crucial step, as numerical groundwater flow models built on flawed conceptual models inevitably produce inappropriate results (Bredehoeft 2005; Robins et al. 2005). Finally, the conceptual hydrogeological model can be evaluated and improved by a numerical flow model. If the calibrated numerical flow model cannot satisfactorily reproduce measured hydraulic heads and discharge rates, the conceptual geological model (definition of geological units, inclusion of faults, conduits geometry, etc.) needs to be revised.

Three-dimensional geological modeling consists of inferring a realistic spatial representation of the lithological bodies of the studied domain and other relevant features, such as faults, from sparse data and interpretations (Calcagno et al. 2008; Hassen et al. 2016; Pakyuz-Charrier et al. 2017; de la Varga et al. 2019). Several sources of data are usually employed, including borehole logs, geological maps, and geophysical data (Wellmann et al. 2010). A 3-D geological model is composed of various surfaces (e.g., strata contacts and faults) defining discrete volumes of geological units, that can subsequently be populated with petrophysical properties (Perrin et al. 2005) such as permeability and porosity. Inaccuracies are commonly related to data density and quality, geological complexity, interpretations, conceptual uncertainty, simplification requirements, and intrinsic randomness (Wellmann et al. 2010; Moya et al. 2014, Enemark et al. 2019).

Three-dimensional geological models are frequently employed to support decision-making, playing an important role in mining (Collon et al. 2015; Guo et al. 2016), oil and gas (Perrin et al. 2005), and geothermal exploration (Milicich et al. 2014; Collon et al. 2015), environmental management and hydrogeological surveys (Ross et al. 2005; Cox et al. 2013; Raiber et al. 2015), design of geotechnical constructions and nuclear underground-storage facilities (Živec and Žibert 2017), and evaluation of underground resources (Cherpeau et al. 2010), among other quantitative geoscientific applications (Hillier et al. 2016; Laurent et al. 2016; Pakyuz-Charrier et al. 2018). The modeling objectives can vary from spatial depiction of geological features to the generation of input data for physical-process simulations (Calcagno et al. 2008). In groundwater studies, 3-D geological modeling has been suitably used to support the conceptualization of flow regimes and to characterize groundwater recharge (e.g., Cox et al. 2013; Moya et al. 2014; Raiber et al. 2015; Hassen et al. 2016; Martinez et al. 2017; Raiber et al. 2019), and in the setup of numerical models (e.g., Fogg 1986; Martin and Frind 1998; Ross et al. 2005; Borghi et al. 2015). Modeling the fault network is of particular importance throughout 3-D geological modeling of bedrock aquifers, because faults can act as sealing barriers or drains, compartmentalizing and considerably influencing the subsurface flow (Cherpeau et al. 2010).

Over the last decades, various methods have been developed to construct 3-D models of geological objects (e.g. Gjoystdal et al. 1985; Lajaunie et al. 1997; de Kemp 1999; Cowan et al. 2003; Calcagno et al. 2008; Wellmann et al. 2010; Zou et al. 2012; Hillier et al. 2016; Burs et al. 2016; Guo et al. 2016; Martin and Boisvert 2017; Gonçalves et al. 2017; Riesner et al. 2017, de la Varga et al. 2019). Traditional 3-D geological modeling schemes are based on manual digitization of geological units within cross-sections. Simple to complex geometries can be generated by manually connecting two-dimensional (2-D) polylines across multiple sections. However, the intensive manual-drawing demand of this technique, known as explicit modeling, makes model updates a tedious and laborious task (Cowan et al. 2003; Vollgger et al. 2015; Guo et al. 2016). This is a substantial constraint, given that geological interpretations inevitably evolve as more data become available (Cowan et al. 2003). Consequently, the number of hypotheses that can be explored by explicit geological modeling while assessing and interpreting data is limited. Alternatively, geostatistical techniques can be employed, producing reasonable results for simple geometries. Complex shapes can be explicitly interpolated using triangulated irregular grids, but closely spaced sections and a high data density are required (Wellmann et al. 2010). EarthVision, GOCAD, Surfer, and Petrel are examples of software that apply an explicit geometric method, in which surfaces are represented by direct triangulation or interpolation (Burs et al. 2016).

Progress in 3-D interpolation techniques has favored the emergence of a practical alternative to explicit modeling, in which geological surfaces are quickly created by translating geological data into numerical values and subsequently applying computer algorithms (Cowan et al. 2003; Caumon et al. 2013; Vollgger et al. 2015; Burs et al. 2016). This scheme considers geological surfaces as isosurfaces of given functions, that are defined everywhere in space and built from the data location and orientation, besides some supplementary external constraints (Wellmann et al. 2010) such as the stratigraphy, fault geometry, hierarchy of stratal units, and contact types (deposits, erosions, intrusions, veins). This modeling approach is known as implicit modeling, because the explicit definition or digitization of surfaces is no longer required (Cowan et al. 2003). Nevertheless, polylines can still be defined to infer boundary positions, which are thereafter interpolated by the algorithm (Cowan et al. 2004). This way, complex geometries with arbitrary topology can rapidly be created by combining field data and geological interpretation (Alcaraz et al. 2011; Hillier et al. 2014). Three-dimensional implicit geological modeling has been gaining popularity over the last decade (Gonçalves et al. 2017). The major benefit of implicit modeling relies on its speed, which allows faster assessments of new information and uncertainty simulations, accordingly reducing the intrinsic risks of geological modeling (Cowan et al. 2003; Wellmann et al. 2010; Birch 2014; Vollgger et al. 2015). By comparing implicit and explicit geological modeling, Birch (2014) concluded that implicit modeling performs better than the traditional method of manual digitization. Additional benefits over explicit models comprise reproducibility, automation, and reduced user-based bias (Vollgger et al. 2015; Gonçalves et al. 2017). Moreover, the implicit function allows a fast and internally consistent 3-D model, which straightforwardly integrates multiple types of data such as lithological contacts and orientation, and structural constraints (Philippon et al. 2015; Vollgger et al. 2015). GeoModeller, Leapfrog, and Vulcan are examples of commercial software with implicit modeling capabilities (Burs et al. 2016).

Validating geological models constructed on restricted subsurface data is difficult; therefore, the consideration of regional insights facilitates the integration of surface observations and subsurface interpretations (Agar and Geiger 2015). Karst aquifers commonly present a hierarchically arranged network of linked fractures and conduits that may be drained by only one main spring (Hartmann et al. 2014). Thus, artificial tracer experiments with multiple injection points can provide valuable information on the structure of the karst-conduit network, especially combined with data on structural geology, spring hydrology, and speleological observations (Perrin and Luetscher 2008). Emanation of soil gases such as Rn, CO₂, and He, are indicative of active faults (Lombardi and Voltattorni 2010). Although natural releases of CO₂ can also be related to metamorphism, decomposition of organic material, and biological activity, elevated CO₂ degassing seems to be associated with mantle emissions, active and ancient volcanism, and decarbonation processes (Lombardi and Voltattorni 2010). In areas with such characteristics, shallow soil CO₂ degassing surveys may be used to delineate fault zones (e.g., Schütze et al. 2012; Lombardi and Voltattorni 2010), which may locally control flow-path directions (Perrin and Luetscher 2008).

Due to their specific heterogeneity, karst aquifers represent a special challenge for quantitative modeling (Hartmann et al. 2014). A multitude of approaches were proposed, from very simplistic black-box or lumped-parameter models in order to simulate rather phenomenologically the hydrologic dynamics of catchments and spring discharge (Butscher and Huggenberger 2008; Jukić and Denić-Jukić 2009; Tritz et al. 2011; Željković and Kadić 2015; Hosseini et al. 2017; and many others) to spatially distributed numerical models, which explicitly describe the groundwater flow using an equivalent porous medium approach (e.g., Ghasemizadeh et al. 2015; Borghi et al. 2016), a dual-continuum approach (e.g., Sauter 1993), or a discrete conduit-continuum approach (e.g., Kiraly 1979; Liedl et al. 2003; Shoemaker et al. 2008; de Rooij et al. 2013; Chen and Goldscheider 2014; Kuniansky 2016). Reviews on available model concepts are presented, e.g., by Kovács and Sauter (2007), Ghasemizadeh et al. (2012), and Kalhor et al. (2019).

This paper demonstrates how a spatially distributed model can be conceptualized, set up, and calibrated to describe groundwater flow in a karst aquifer by combining field surveys, 3-D implicit geological modeling, and an iterative model calibration method to reasonably consider hydraulically relevant karst conduits in the model. The paper describes the principal approach and its application to a Triassic limestone aquifer in southwestern Germany. The discussion sets a focus on the value of combining, integrating, and processing information and data from various sources to develop a geological and subsequently a hydrogeological model. This is done for the example of a preliminary assessment of possible risks from limestone quarrying for which a steady-state equivalent-porous-medium model was set up and calibrated to inform a particle tracker in order to estimate groundwater flow directions from existing quarries.

Principal approach

The proposed workflow combines the compilation of data, implicit geological modeling, field work, and the development of both a conceptual hydrogeological model and a numerical groundwater flow model, as visualized in Fig. 1. The first phase comprises the acquisition, compilation, and potential digitalization of topographical data (digital elevation models, DEMs), geological data (borehole logs, geological maps, etc.), and hydrogeological data (well construction profiles, hydraulic heads, location of springs and discharge rates, etc.). For the construction of a 3-D geological model, the surfaces to be represented in the model must be chosen in accordance with the main objectives of modeling and the given scale of the study area (Perrin et al. 2005). In order to define the site-specific hydrostratigraphy, hydrogeological relevant units are identified and, if appropriate, units with similar hydraulic properties are merged. If significant uncertainties arise or new interpretations are proposed, supplementary field work can be conducted to test the correctness of the geological model. Once the arrangement of geological objects in the model respects regional and local contexts, the history of possible geological events, and the field observations, i.e. once geological realism is considered reasonable based on the available data and knowledge, the model can be considered validated and is ready to be used to describe the primary spatial distribution of the properties (e.g., hydraulic conductivity, porosity, storage coefficient) of the numerical model. Additionally, the spatial representation of lithological entities can be used to assist the visualization of filter screens and hydrogeological data in three dimensions, e.g., potentiometric surfaces and tracer-test results. On this basis, a conceptual hydrogeological model can be developed that provides the necessary information for defining the boundary conditions of the numerical model. The selection of an appropriate numerical model approach should be part of this conceptualization. During the calibration process, further refinements of the numerical model’s spatial parametrization can be conducted in an iterative way, i.e. using the geological model for visualization of interim simulations, including conceptualizations of the karst-conduit network. Eventually, the results of numerical modeling lead to additional geological insights. In early stages, a common conclusion drawn from the model results is that the existing numerical model cannot be calibrated. This may call for a review of the suitability of the chosen model approach or a revision of the conceptual hydrogeological model, or it requires a modification of the geological model. Such modifications may also be necessary when new data become available that are in contradiction to the current geological model. Depending on the data type and the magnitude of the misfit, new data may require only a re-calibration of the numerical model, a revision of the hydrogeological conceptualization, or a revision of the underlying geological model.

Fig. 1

Proposed workflow. Solid lines: mandatory steps; dashed lines: optional steps

Study area

Overview and previous work

The outlined workflow was applied to a karst aquifer in south-west Germany. The study area is located between the mountain ranges of the Black Forest and the Swabian Alb (see Fig. 2). It covers the catchment of the Ammer River, a tributary of the Neckar River, and relevant neighboring catchments (total area approx. 750 km²). The geology comprises a sequence of Triassic and Jurassic strata with several prominent escarpments. The layers of the Upper Muschelkalk (mo), having a total average thickness of approx. 80 m, together with the underlying 10–12-m-thick porous-to-cavernous uppermost dolomites of the Middle Muschelkalk (mm), form the main aquifer, which is used for the supply of drinking water by several municipalities. Shallow aquifers in fluvial deposits exist only locally.

Fig. 2

a Location and geology of the study area, b cross-sectional view of geological units in the study area, and c location of particular features within the study area; coordinates: ETRS89-UTM32N; source: LGRB (1998)

The work presented in the following builds on the achievements of previous research addressing the geological and hydrogeological situation in the catchment of the Ammer River and its surroundings. Harreß (1973) mapped karst features and springs in the study area and carried out a large number of tracer tests to evaluate the risk from potential waste-water infiltration into the Muschelkalk aquifer. Reuther (1973) examined the tectonics and stratigraphy, and mapped the major faults in the study area. Villinger (1982) analyzed the tracer-test data of Harreß (1973) and further sources (unpublished), in combination with discharge measurement at the main springs and estimates of groundwater recharge, to consistently balance the groundwater budget of the catchment. Plümacher (1999) and Plümacher and Ufrecht (2000) described the regional groundwater flow in the Muschelkalk aquifer for a much larger area (ca. 4,500 km²). Selle et al. (2013) presented a numerical groundwater-flow model for the upstream part of the catchment, without considering any faults. Pavlovskiy and Selle (2015) analyzed environmental-tracer data at selected groundwater-abstraction wells to infer well-catchment areas.

Geological setting

In the study area, Triassic and Jurassic strata gently dip toward SE (circa 1°), due to the Tertiary subsidence of the Alpine foreland Molasse basin and the uplift of the Black Forest crystalline basement (LGRB 2005). The elevation predominantly ranges from 300 to 600 m above sea level. Only in the southeastern corner of the study area, where the Middle Jurassic crops out, does the topography reach elevations of up to 800 m.

Deposited in a cratonic basin, the Triassic sediments are subdivided into (1) Lower Triassic fluvial Buntsandstein, (2) Middle Triassic marine Muschelkalk, and (3) Upper Triassic partly marine, partly terrestrial Keuper (Hornung and Aigner 1999). The deep marine Jurassic sediments occur only at the eastern parts of the study area (Fig. 2), where the Lower and Middle Jurassic units are exposed. Pleistocene glacial unconformities overlaying the Mesozoic sediments frequently occur on the hills (loess) and hill slopes (talus), as well as in the valleys (drift sheet). The latter is mainly covered by Holocene alluvial deposits (GLA 1966; LGRB 1998). Figure 3 provides a short description of the main geological units outcropping in the study area.

Fig. 3

Stratigraphy of the study area. a Geological map (1:1,000,000): LGRB (1998). b Geological maps (1:25,000): GLA (1966, 1986, 1989, 1992, 1994a, b, c), LGRB (1996, 2005).

Two major fault systems characterize the central Europe tectonics: (1) the pre-Variscan (characterized by the Erzgebirgisch SW–NE and the Swabian WSW–ENE direction) and (2) the post-Variscan (containing the Hercynian SE–NW and Rhenish SSW–NNE trends)(Petrovic 2016). The study area is specifically influenced by conjugate fault patterns formed by two large synthetic WSW–ENE fault zones, the Swabian lineament in the south and the Neckar-Jagst lineament in the north (Fig. 4), including antithetic-fault confluences and intersections in between (Schwarz and Kilfitt 2008). Extensional, sinistral and high-angle SE–NW faults branch off to the north of the destral Swabian lineament, forming the Filder-Graben system, with three hanging walls. The footwall comprises the major portion of the study area (Fig. 4). Close to Tübingen, the SE–NW fault zone separating the footwall from the upper hanging wall forms a confluence bow with branches while merging with the Swabian lineament, revealing rotational effects (Schwarz and Kilfitt 2008). In this linkage zone the Swabian lineament forms a graben, also known as Bebenhausen graben, with offsets varying from 80 to 100 m. Commonly, the offsets along the Swabian lineament ranges from 20 to 30 m (Ufrecht 2006). The faults of this region are additionally characterized by a preferential degassing of carbon dioxide, probably emerging from deeper parts of the earth’s crust (Harreß 1973; Schütze et al. 2012).

Fig. 4

Fault lines in Baden-Württemberg with details of the study area (source: LGRB 1998)

A karst system initiates and develops through a combination of processes such as dissolution (subrosion), erosion (corrosion) and incasion (collapse). Fractures are widened by dissolution and at a certain stage a critical point is achieved, with turbulent flow occurring and accelerating the enlargement of fractures. Eventually some fractures start to dominate the dissolution processes and large conduits and caves may evolve (Benson and Yuhr 2016). The area of potential karstification due to carbonate dissolution by infiltrated water is limited to the western part of the study area (Fig. 5), where the Upper Muschelkalk and uppermost dolomites of the Middle Muschelkalk are exposed or only weakly covered. It forms the major recharge area and is characterized by sinkhole fields, low drainage density, influent streams, and karst springs. Further towards ESE, i.e. following the strata dip, the Keuper overburden increases and the Muschelkalk aquifer gets increasingly confined (Villinger 1982; Ufrecht 2006), preventing the Upper Muschelkalk karstification (Pavlovskiy and Selle 2015). Karstification may also be locally favored by the degassing of mantel-borne CO₂ along faults (Harreß 1973), providing an additional structural control to the development of the karst-conduit network. Karstification also occurs in the gypsum-bearing Grabfeld formation (kmGr) due to dissolution of gypsum and anhydrite, but the Grabfeld formation and the Upper Muschelkalk are separated by the Erfurt formation (kuE), which acts as an aquitard so that karstification of the Muschelkalk limestones and the Grabfeld gypsum is decoupled.

Fig. 5

Area of karstification and location of karst features (sinkholes, karst depressions, etc.) within the study area—as provided by the State Geological Survey (LGRB 2014); colors indicate data origin of karst features: brown = digital elevation model (DEM), green = soil map

Sinkholes are by far the most frequent karst feature in the study area. The biggest one (the Herrgottsscheuer) is 30 m wide and 10 m deep. The open shaft is 14 m deep and SE–NW oriented. Only two caves are accessible. The biggest one (the Pommerlessloch) is at least 50 m long and SSW–NNE oriented (Harreß 1973). Although the described shaft and cave clearly present tectonic control, i.e. pre-Variscan and post-Variscan orientations, respectively, the existing knowledge about the karst system does not allow one to discuss the role of individual joints and faults on the karstification process. Moreover, considering the classification proposed by Waltham and Fookes (2003), the karst system can be classified as youthful, with drainage control and small sinkholes and caves.

Hydrological and hydrogeological settings

The annual mean precipitation in the study area ranges between 700 and 850 mm—Fig. S1 of the electronic supplementary material (ESM). Maximum values typically occur in summer. Storm events occur only rarely (average annual number of days with precipitation >10 mm/>20 mm/>30 mm is approx. 20/4/1). The annual average air temperature is approximately 8 °C.

The Upper Muschelkalk and the 10–12-m-thick permeable uppermost dolomites of the Middle Muschelkalk (Diemel Formation, mmD) form the Muschelkalk aquifer (Villinger 1982; Ufrecht 2006), the main groundwater reservoir in the study area. Internally, a 6–8-m-thick continuous limestone-mudstone sequence occurring at the bottommost part of the Upper Muschelkalk, known as Haßmersheimer layers (moH), hydraulically split the Muschelkalk aquifer into a basal and an upper part (Villinger 1982; Ufrecht 2006). Groundwater in the Muschelkalk aquifer mainly recharges in the outcropping areas in the western part of the study area, and roughly flows in a SE direction, following the strata dipping (Villinger 1982). In areas where the Muschelkalk is not covered by low-permeability strata, the aquifer recharge rate is approximately 220 mm/year.

Hydraulic conductivity values vary over several orders of magnitude for all formations considered. In the Upper Muschelkalk, hydraulic conductivity ranges from approx. 10⁻⁶ to 10⁻³ m/s (Villinger 1982). The mean hydraulic conductivity in the Ammer area may be estimated to 4 × 10⁻⁵ m/s according to the transmissivities reported in Villinger (1982) (see also Plümacher 1999). The hydraulic conductivity of the Grabfeld formation (kmGr) depends on its weathering state and how strong gypsum has been leached. kmGr aquifers may have a very complex internal structure where, in extreme cases, impermeable, low-permeability and high-permeability zones can be found contiguously (Ufrecht 2017). Hydraulic conductivity values in leached and nonleached kmGr units range from 10⁻⁹ to 10⁻⁴ m/s and from 10⁻¹³ to 10⁻⁵ m/s, respectively (Schlosser et al. 2007), with mean values between 10⁻⁶ and 2 × 10⁻⁵ m/s (Ufrecht 2017). The hydraulic conductivity of the Quaternary sediments within the Ammer and the Neckar valleys is on the order of 10⁻³ m/s (e.g., Lessoff et al. 2010). Although sand and dolomite facies can build thin and localized aquifers within the Keuper units, the km2, km3, km4, km5, km5, ko, as well as Lower Jura (ju) units play only a minor hydrogeological role in the study area, as they predominantly occur in forest areas at the Schönbuch plateau, where groundwater recharge is assumed to be low (Selle et al. 2013). Note that the drainage pattern clearly follows a WSW–ENE and SE–NW pattern (Fig. 6). Further eastwards, however, dendritic drainage patterns progressively develop, given the low permeability of the Upper Keuper and Lower Jurassic units.

Fig. 6

Hydrological regime of the study area according to spring discharge data, river run-off data and tracer-test results (distinction between minor and major transport directions on a qualitative basis, see Harreß 1973). The subsurface catchment does not coincide with the topographic catchment of the Ammer River; areas of losses of water runoff to neighboring catchments, as well as of runoff gains from neighboring catchments, have been estimated based on tracer test results; numbers indicate IDs of zones used in the modeled water balance

The evaporites of the Middle Muschelkalk characterize the lower boundary of the Muschelkalk aquifer. The Middle Muschelkalk is an aquitard that, together with the Lower Muschelkalk, regionally separates the Muschelkalk aquifer from the underlying Buntsandstein aquifer (Ufrecht 2006).

Over the last decades, a total of 35 tracer tests have been conducted in sinkholes, losing rivers, and dry valleys of the study area, in order to better characterize the Muschelkalk aquifer (Harreß 1973; Villinger 1982; GLA 1997). These tests reveal high flow velocities (~100 m/h on average, >200 m/h maximum) and clearly show that the Ammer surface watershed mismatches the groundwater divides. Tracer compounds injected northward of the topographic divide were recovered in the Ammer catchment while tracers injected in the western part of the Ammer basin were found at the Bronnbach spring in the Neckar Valley south of the topographic Ammer catchment (Fig. 6).

The Muschelkalk aquifer naturally discharges to the Neckar and Ammer rivers, including their associated tributaries (Villinger 1982; Selle et al. 2013). Major perennial springs of the Muschelkalk aquifer are the Bronnbach spring (mean discharge: 530 L/s) at the southern edge of the study area, the Ammer spring, the Herrenberg spring (each approx. 15 L/s), and the Reusten spring (~10 L/s). These springs are examples of focused groundwater discharge through karst features. The Muschelkalk aquifer is used for drinking water supply. Groundwater abstraction wells are situated both in the Ammer and Neckar valleys with average pumping rates of 95 and 55 L/s, respectively (Ammertal-Schönbuchgruppe (ASG), unpublished data, 2012).

Low sulfate concentrations in these production wells suggest that the hydraulic connection between the Grabfeld formation and the Muschelkalk aquifer is weak, meaning that the lower-Keuper Erfurt formation (kuE) act as an aquitard prompting groundwater discharge from the overlaying Grabfeld formation aquifer (Pavlovskiy and Selle 2015). Using environmental tracers (tritium, SF6, and temperature), Pavlovskiy and Selle (2015) estimated mean water transit times to the groundwater abstraction wells ranging between 6 and more than 50 years.

3-D geological modeling

Software and modeling procedure

The geological modeling was carried out using the 3-D implicit modeling software Leapfrog (ARANZ Geo Limited 2016), which can be operated on conventional personal computers (Vollgger et al. 2015). Leapfrog uses the data and parameters supplied such as lithological codes, DEMs and drawn polylines, to implicitly construct surfaces based on spatial interpolation of the point attributes using radial basis functions (RBF) interpolation (Cowan et al. 2003; Alcaraz et al. 2011; Hillier et al. 2014; ARANZ Geo Limited 2016), which is mathematically identical to the function-estimate form of kriging. Multiple models conditioned on data can be built (Cowan et al. 2003). Leapfrog’s hydrogeology module transfers the implicit geological model to input files of MODFLOW and FEFLOW.

Data research and compilation

Data on surface geology, including lithological contacts and tectonic structures, were gathered from nine geological maps at the scale of 1:25,000, in raster (GLA 1966, 1986, 1989, 1992, 1994a, b, c; LGRB 1996, 2005) and vector (GIS) format at the scale of 1:50,000 (LGRB 2012). For the subsurface geology, 463 boreholes were examined (LGRB 2020), in addition to known heights of the Quaternary unconformity at the Neckar Valley, between Rottenburg and Tübingen (LGRB 1999). The terrain’s surface was represented by a DEM with 10-m resolution (LGL 2012).

The quality and relevance of the data were systematically evaluated. From the 463 boreholes inspected, 78 boreholes were subjected to stratigraphic interpretation and classification, while 53 were disregarded due to missing data or ambiguity. Finally, 410 irregularly distributed points were considered in the geological modeling (Fig. 2c), with depths varying from 1 to 392 m. Although the average depth is 27 m, approximately 77% of the boreholes are shallower, 50% are 20–30 m deep and only 4% are deeper than 100 m.

Definition of the model domain

The geological model domain, i.e. the study area (as outlined in Fig. 2c), was defined such that given physical boundaries of the Muschelkalk aquifer are encompassed such as the bottom of the Muschelkalk aquifer in the western and northern parts of the study area, and the footwall border in the eastern part (Figs. 2 and 3), considering the offsets greater than 80 m, at the northern portion of the SE–NW fault zone and along the Bebenhausen graben (Ufrecht 2006). To make sure that the geological model will fully cover the anticipated domain of the hydrogeological model, the southern boundary was set south of the Neckar River.

Setting the site-specific hydrostratigraphy

The geological model incorporates the entire local stratigraphy of the Triassic (Buntsandstein, Muschelkalk, and Keuper), and the Lower and the Middle Jurassic, as well as the major Quaternary alluvial and talus bodies. Loess deposits, which are typically thin in the study area, were not individually considered in the model but are accounted for in the quantification of recharge. Because the Lower and Middle Muschelkalk mostly act as aquitards, they were unified. Similarly, the km2, km3, km4, km5, km5, and ko units were also grouped, as well as the Jurassic units.

3-D implicit geological modeling

Given that geological and geophysical data are inherently uncertain and often scattered, no geological representation can be inferred from the data without a preceding geological interpretation (Perrin et al. 2005). In order to gather an initial understanding of the site geology, as well as to identify inconsistencies, all imported data were promptly visualized and evaluated in three dimensions. Unlike the explicit modeling method, which typically requires the construction of cross sections and the explicit digitization of surfaces, the surfaces are automatically generated from the spatial distribution of lithological contacts and orientation in 3-D implicit geological modeling. As the boreholes are rather irregularly distributed in space (Fig. 2c), the modeling process required significant interpretation to assist the model interpolation in data-scarce areas, which were manually added by drawing points and polylines inferring geological contacts in the subsurface. Hence, the development of the model was driven by both data and geological knowledge.

The modeling process started with a very simple setting, i.e. discontinuities (faults) were successively added as necessary to achieve a parsimonious model, as simple as possible and as complex as necessary. Faults can be defined in Leapfrog using polylines or GIS vectors to create vertical walls or surfaces, interacting with each other based on the faults’ chronology and interactions. Faults divide the geological model into individual blocks. Each block has its own and independent surface chronology, attributes, and spatial interpolation, meaning that hard data (borehole logs and lithological contacts extracted from the geological maps) and inferences of one block do not influence the surfacing of neighboring blocks.

Faults were manually incorporated in the model by firstly digitizing polylines representing their strikes on the DEM. Due to the high-angle characteristics of the faults in the study area, all faults in the model were represented by vertical walls (dip = 90°) for the sake of simplicity. In Leapfrog, a fault must extend to another fault or to the model boundaries. Fault extensions are modeled based on the interpreted fault chronology (e.g., older faults are not allowed to cross younger faults). Where modeled faults exceeded their known fault extension, care was taken to avoid offsets between neighboring blocks beyond the known faults’ length. In data-scarce areas, subsurface lithological boundaries were inferred below mapped contacts based on the estimated average thickness of the strata.

Supplementary field work, model refinement

Field campaigns were conducted to investigate areas where the available data presented ambiguity or high uncertainty as well as to verify raised new hypotheses. For instance, in the central part of the study area, discrepancies between the 1:25,000 geological maps in raster at the scale of 1:25,000 (GLA 1966) and vector formats at the scale of 1:50,000 (LGRB 2012) were verified, with respect to the location of the top of the lower Keuper (kuE) formation on the Ammer Valley slopes. A gamma-ray well logging, conducted in a 60-m-deep artesian well, supposedly cutting the lower Keuper (kuE) and Muschelkalk aquifer (results not shown), indicated that the raster map (GLA 1966) indeed better depicts the local geology, which was then used as a reference to control the geological model outputs.

In another area, the log description of two neighboring boreholes, the surface geology and the local geomorphology (local creek alignment and direction of the escarpment erosion front) indicate the existence of an unmapped S–N oriented fault (Fig. 7). This postulated fault, however, was not in accordance with the site synthetic WSW–ENE and antithetic SE–NW fault trends. In order to assess the plausibility of the suggested S–N oriented fault, dip and strikes measurements were carried out at a former gypsum quarry close to Käsbach Creek. The identified high-angle SW–NE to SSE–NNW fractures, in agreement with synthetic and the postulated fault trend, respectively, provided substantial indications that the existence of a N–S fault is reasonable. The postulated S–N fault was finally validated by two subsequent soil-CO₂-mapping profiles, considering that site faults commonly emanate mantel-borne CO₂ (Harreß 1973; Schütze et al. 2012) and the fact that volumetric CO₂ concentration peaks (~ 13%) in the soil gas were measured a few meters apart from the interpreted fault location (Fig. 7c,d). This fault may represent an S–N element, connecting two synthetic faults, which originated at a late stage of internal rotation of the footwall block in a transpressional environment. A detailed analysis of this tectonic feature is beyond the scope of this paper.

Fig. 7

a Postulated, verified and known faults and tectonic elements (raised and down-dropped blocks) in the center of the area of interest (see Fig. 2c for a better orientation), b stereographic projections of measured fractures, and c northern and d southern soil CO₂ profiles crossing Käsbach Creek

Further unmapped faults were inferred, considering the surface geology and borehole data, where available (Fig. 7a). These faults were identified by unrealistic shapes, e.g. ripples, on the modeled surfaces.

A pair of long and curved faults (ESE–WNW) was interpreted forming a down-dropped block with modest offsets (<20 m) just southward of the Mötzingen quarry (see western part of Fig. 7a). This interpretation is in agreement with Reuther (1973) who also inferred the existence of such faults using soil CO₂ mapping. Junginger (2019) recently measured radon concentrations at the Ammer River in Poltringen (see central part of top-left map in Fig. 7), where a fault was postulated, providing additional evidences that one or more active faults intersect the Ammer River and that the interpretations upon which the model was constructed are reasonable.

The postulated antithetic SE–NW faults are roughly parallel to the direction of maximum compression of the active stress field (about 140°), from the Alpes towards its northern foreland, according to Illies et al. (1981). The interpreted synthetic WSW–ENE curved faults are parallel and associated to the Swabian lineament.

In Wolfenhausen, the postulated synthetic faults indicate the existence of a local down-dropped block with offsets around 80 m, as was also inferred by Villinger (1982). Close to Wurmlingen, however, the postulated synthetic faults form a raised block. Maybe these slip-strike features represent negative and positive flower transpressional structures, respectively.

Once the main lithological surfaces and model features were set, the relevant hydrofacies, i.e. the upper dolomite of the Middle Muschelkalk (mmD) and the Haßmersheimer layers (moH) of the Upper Muschelkalk, were modeled, using the top of the Middle Muschelkalk surface as reference and applying downward and upward constant offsets, respectively. A constant thickness of 12 m was assumed for mmDol. Although moH consist of an upper and lower marl layer, separated by limestone, they were merged to a single 6-m-thick layer in the model, occurring 7 m above the bottom of Upper Muschelkalk (mo). Figure 8 shows the final geological model, which served as the basis for the subsequent hydrogeological modeling.

Fig. 8

Three-dimensional representation of the geological model

Hydrogeological modeling

Conceptual model

In order to develop the conceptual hydrogeological model of the study area, all relevant data and information were compiled, processed, and visualized within Leapfrog. Vadose, phreatic, and confined zones were delineated by intersecting the hydraulic-head surface with the elevations of the aquifer base and top (as represented in terms of geological contacts in the 3-D geological model). Figure 9 shows the overlay of these features, providing a clear picture of the situation.

Fig. 9

Hydrogeological conditions in the Upper Muschelkalk aquifer (3-D representation within Leapfrog). The dashed blue line indicates the transition between vadose and phreatic conditions, while the solid blue line indicates the transition between phreatic and confined conditions for average hydrologic conditions; note the two indentations corresponding to horst structures and the isolated confined area, due to a down-dropped block near Wolfenhausen with expressive faults offset (see Fig. 7). Black and white lines represent tracer test results with major and minor recovery rate, respectively (Harreß 1973, Villinger 1982). Faults appear as translucent gray walls. Beyond the Filder graben system in the northeast, the entire overlaying model stratigraphy is shown

With the Triassic layers progressively dipping towards the southeast, hydraulic conditions in the Muschelkalk aquifer change in the same direction, from vadose (i.e. unsaturated) to phreatic, and finally confined conditions. Along the Neckar Valley in the south of the model domain, northern and southern groundwater flow regimes converge. Given the offsets of the extensional antithetic SE–NW faults, separating the footwall from the upper hanging wall, as well as the Bebenhausen graben along the Swabian Lineament, it can be assumed that groundwater outflow in the study area occurs under confined conditions between the Swabian Lineament and the Neckar Valley.

In the most western part of the study area, where the Upper Muschelkalk is exposed or its overburden is thin, the Upper Muschelkalk is a vadose zone. This region is the major recharge area of the modeled aquifer. Here sinkholes emerge due to carbonate dissolution and underlying subrosion of evaporites in the Middle Muschelkalk (Fig. 5). Following Hartmann et al. (2014), such sinkholes, with their associated subsidence fractures, probably lead infiltrated water directly to the basal part of the Muschelkalk aquifer, thus recharging groundwater below the Haßmersheimer layers. The level of the water table is in the upper dolomite of the Middle Muschelkalk, and the saturated thickness and hydraulic transmissivity of the aquifer are rather low. In dry periods, the groundwater body might even be absent in some areas. In the southeasterly direction, the thickness of the saturated zone continuously increases, and flow in the Upper Muschelkalk is first phreatic below, and then above the Haßmersheimer layers (where it starts to be confined below them). Further southeast, the Muschelkalk aquifer becomes fully confined, with the two horsts being exceptions.

While a large part of groundwater is recharged in the basal part of the Muschelkalk aquifer, it discharges predominantly in the upper part (springs, wells and rivers are all located above the Haßmersheimer layers). This applies also to the tracer tests, which were predominantly injected in sinkholes and most likely reached the basal aquifer part, but were recovered in receptors located above the Haßmersheimer layers at high recovery rates (up to 80%). The 3–D geological model helps to determine likely tracer propagation paths. To illustrate this by an example, a tracer test is considered, that showed a propagation over the distance of about 13 km within ~5 days from the point of injection in the western part of the study area to the Bronnbach spring (labelled as “Test 2017” in Fig. 9). The injected tracer possibly migrated through the basal part of the Muschelkalk aquifer until reaching an extensional fault that hydraulically connects basal and upper parts of the Muschelkalk aquifer (Fig. 10). Further downgradient, the tracer most likely bypassed the Wolfenhausen graben (see Fig. 9), migrating eastward of it and reaching the Bronnbach spring through the postulated antithetic fault. Please note that similar conceptualizations may be done for other tracer tests conducted in the western part of the study area.

Fig. 10

Estimated path of tracer propagation during the test in 2017 from injection to observation point (Bronnbach spring) with help of the 3-D geological model. The dashed line indicates the circumvention of the Wolfenhausen graben. (For the color code of the geological units see Fig. 8)

The flow lines reflected by the tracer experiments are oriented mostly in the SE–NW direction (Fig. 9). Although in shear-zone domains, synthetic faults tend to be partially opened and antithetic faults compressed (Schwarz and Kilfitt 2008), the combined alignment effect of the SE bedding dip and SE–NW antithetic faults, as well as the kinetically related fracture sets, probably played an important role in the competitive formation of karst conduits. Bedding planes may play a more significant role for karstification than joints and faults, as their greatest apertures tend to be more continuous than those of fractures (Ford 2003). Moreover, the hydraulic properties of a fault core zone may change over time, initially acting as a conduit during deformation, and subsequently as a barrier, if mineral precipitation, clay or even breccia clogs the fault plain (e.g., Caine et al. 1996, Bense et al. 2013).

Though the majority of the tracer tests shows similar results with respect to the observed directions of propagation, they differ with respect to propagation distance and velocity as well as recovery rates. The individual analysis of these data in the context of the specific setting of the tracer test (i.e. aquifer conditions at the location of tracer injection and observation) provides essential information for the parameterization of the numerical model. High values of both velocity and recovery rate, for example, are a clear indication for a hydraulically very effective karst feature linking the injection and observation points. These features need to be explicitly accounted for in the numerical groundwater-flow model if tracer tests are to be reproduced by subsequent particle-tracking simulations.

In the analysis of individual tracer tests, the hydraulic setting along the presumed tracer paths is also important. Some tracer tests showed that the tracers may have migrated also through the confined parts of the Muschelkalk aquifer to the Reusten spring (Fig. 9). Karstification under confined settings generates different morphologies than under unconfined settings. Specific hydrogeologic characteristics under confined conditions (restricted inputs/outputs) reduce positive flow-dissolution feedbacks and, consequently, the karstic competition in fracture networks. As a result, more pervasive channeling is likely to be developed under confined settings, i.e. 2-D or 3-D densely packed maze conduits, in contrast to conduit-like morphologies that tend to be broadly spaced and dendritic due to highly competing developments under unconfined conditions. Hence, the values of porosity and specific storage in confined karst systems may be up to one order of magnitude larger than in unconfined settings (e.g., Klimchouk 2006), which should be reflected through larger tracer transit times. Observed tracer migration velocities (based on tracer arrival at the Reusten spring, see GLA 1997), however, suggest an alternative tracer transport path along SE–NW fractures and WSW–ENE synthetic faults (which are connected to WSW–ENE faults), this way bypassing the confined area. This option of explaining the flow regime is also supported by the fact that the estimated apparent water age in the Poltringen wells is relatively young (Pavlovskiy and Selle 2015). However, a valid clarification is not possible—as in almost all issues concerning the role of fractures, faults, and discrete karst conduits for the flow regime. This particularly applies to local effects, which can hardly be resolved in a catchment-scale model. In contrast, rather regional effects of discrete hydraulically effective features, acting as an internal drainage system that controls the overall water budget, can be better identified and described in a quantitative way as part of the model calibration (see the following).

Numerical groundwater-flow model

The numerical groundwater-flow model is set up and operated as a steady-state model to provide average groundwater flow data for subsequent particle-tracking simulations. The model follows the single-continuum equivalent porous medium approach where high-permeability features such as karst conduits are simply represented by a series of model cells with extra-high conductivity values. Such models are also referred to as smeared conduit models as they do not attempt to model the detailed geometry of relatively small high-permeability features but try instead to capture the effects of a conduit on a much larger computational grid block by assigning effective properties to that grid block (Green et al. 2006).

The model encompasses the entire subsurface catchment of the Ammer River and covers an area of ~470 km² (Fig. 11). The model boundaries are set as follows: (A) the western and northern boundary, as well as parts of the southern boundary, coincide with the line where the bottom of mmD crops out, (B) the northeastern boundary follows the Hildrizhausen fault, Bebenhausen graben, and Swabian lineament where the huge offset of layer elevations blocks groundwater flow, (C) the southern boundary runs along the Neckar River, (D) short sections of the model boundary in the southwest and east follow the expected direction of groundwater flow in the Muschelkalk aquifer, (E) in the north the boundary runs along the Würm’, and (F) in the east the model boundary connects boundaries B and C. The outcrop boundary (see A) is represented by a head-dependent flux boundary (i.e. ‘drain’) at the mmD bottom. The northeastern boundary (B) and the boundaries that follow the assumed groundwater flow direction (D) are defined as no-flow boundaries for both the shallow and the Muschelkalk aquifer. The boundaries along the Neckar and Würm rivers (C and E) are represented as head-dependent flux boundary (i.e. ‘river’) for the shallow aquifer and as no-flow boundary for the Muschelkalk aquifer. In the east, a fixed head boundary describes groundwater flow in the Muschelkalk aquifer across the model boundary (F) with the hydraulic head set equal to the level of the Neckar River. Note that the Ammer as well as its tributaries are represented as ‘drains’ to account for the fact that many stretches of the rivers and creeks are intermittent. With this, any infiltration of surface waters was deliberately neglected, which is, although temporarily observed in the study area, presumably of minor relevance for the steady-state description under average conditions. In the absence of detailed information, the conductance of both drain and river cells were set as high as needed to simulate the assumed good hydraulic contact between groundwater and surface water such that water exchange occurs under small hydraulic potential differences.

Fig. 11

Model domain, recharge distribution and lateral boundary conditions (HDFB = head-dependent flow boundary)

Groundwater recharge is set according to a 500 m × 500 m raster of average rate values (within period 2001 to 2015), provided by LUBW (2016), the state authority for the environment, measurement and conservation in Baden Württemberg that runs the regional soil-water balance model GWN-BWN (see Gudera and Morhard 2015). The recharge rate varies spatially between 20 and 280 mm/year (Fig. 11). A simplification has been introduced in the northwestern part of the model where groundwater flow presumably takes place only in the shallow, near-surface parts of the subsurface and runs off in the Goldersbach River. In this area, the groundwater flow model considers only the Muschelkalk layers, and the groundwater recharge is set to zero due to the thick overburden.

MODFLOW-NWT (Niswonger et al. 2011) was applied using a regular, horizontally layered grid of 100 m × 100 m × 5 m sized cells (Fig. 12). The model incorporates the upper dolomitic part of the Middle Muschelkalk (Diemel formation, mmD), the Upper Muschelkalk (Trochitenkalk formation, moTK, including the Haßmersheimer layers, moH, and the Meissner formation, moM), the Lower Keuper (Erfurt formation, kuE), the gypsum layer in the Middle Keuper (Grabfeld formation, kmGr), further layers of the Middle and Upper Keuper (undifferentiated, km2345o), the Jurassic limestone (J), and the Quaternary layers. The vertical grid resolution of 5 m is mainly determined by the thickness of the moH and is just sufficient to represent this unit in the model as a continuous layer. Due to the inclination of layers the model domain spans vertically a total number of 132 layers (from −30 m to 630 m asl), each with a different number of active model cells. This regular gridding was chosen because a representation of the stratigraphy with layers of variable thickness showed to be infeasible due to the large number of faults with partly large offsets.

Fig. 12

Geometry of the 3-D finite difference groundwater flow model (active model cells) with regular discretization and assigned parameter zones. View from southwest, local coordinates in meters

The basic model setup including the zonation, i.e. the assignment of code numbers to all model cells according to stratigraphic layers is done automatically as part of the functionalities of Leapfrog. In this step, the mo and kmGr layers were further differentiated in terms of their anticipated karstification and weathering or leaching status, respectively. The zonation within the Upper Muschelkalk (mo) distinguishes two zones: one zone where karstification of mo is likely, and one zone where, due to an overburden (kuE and kmGr), no relevant karstification of mo is expected (see Fig. S2 of the ESM for an illustration). The Grabfeld formation (kmGr) was differentiated into two zones: one zone where gypsum is expected to have been leached out in the past, and which is assigned to the upper 10–15 m of kmGr, and one zone where most likely no leaching has taken place (Fig. S3 of the ESM). In total, this hydrostratigraphic zonation yielded 13 different zones (Fig. 3).

The model was calibrated at steady state in the sense of an approximate model calibration. The goal was to achieve a model that reasonably reproduces the hydraulic head measurements at the monitoring wells in the Muschelkalk aquifer (see Fig. 11) and the estimated water balance of the Ammer catchment including parts of its neighboring catchments. As data and properties of the model boundaries were largely fixed (see above), the calibration was done with respect to the spatial distribution of hydraulic conductivity only. The calibration process consisted of two steps that are iteratively run through: (1) estimation of the hydraulic conductivity values of each zone that has been differentiated in advance (as discussed previously), (2) incorporation of high-conductivity cells to mimic discrete karst conduits to ‘correct’ the model parameterization resulting from step 1. This second step has shown to be necessary to properly model the water balance in the domain, particularly the discharge of the Bronnbach spring. With a pure zonal calibration, no reasonable result could be achieved. Obviously, there is no clear evidence or knowledge about the location of such conduit-like karst features but information about karstification (Fig. 5) and about flow direction from tracer tests (see Figs. 6 and 9) gave a good orientation for reasonable guesses. Particularly, tracer tests suggest that karst conduits are more evolved along antithetic fractures (joints and faults), whose strikes are roughly parallel to the dipping strata. Furthermore, the calibration was supported by particle tracking simulations using MODPATH (Pollock 1994) to visualize catchment areas in the model domain (not shown). This way, both high-K features and zonal K values were estimated. Figure S2 of the ESM shows the discrete karst features that have been delineated in the calibration process. Estimated zonal K values (Table 1) show appropriate values that are all within the reported ranges (see the preceding). The comparison of modelled and measured hydraulic head values shows a reasonable agreement (Fig. 13) given that no further spatial differentiation of K within the hydrostratigraphic zones has been introduced during the calibration. Local-scale variations, as are reflected by the three monitoring well clusters at the quarries (see Fig. 13: Sulz1 to Sulz7, Has1 to Has6, and Moetz1 to Moetz5), cannot be properly resolved by the model as expected, but regionally the model seems to appropriately describe the hydraulic situation in the Muschelkalk aquifer. The largest deviation (~17 m) exists at Hailfingen (Hailf1) where, so far, measurements could be made available only for a half-year period in 2019. The calculated mean (406.5 m asl) may not represent average conditions (note that Villinger (1982) mentioned a value of “390–396 m asl”). At Wurmlingen (Wurm), the deviation is most likely related to the definition of the model boundary in this area. As there is no clear information available about the flow conditions in the Muschelkalk aquifer below and along the Neckar River, a no-flow boundary was set in the Muschelkalk (see the preceding)—assuming that flow is parallel to the Neckar, and groundwater may discharge only across the eastern boundary or to the Neckar after passing through the Erfurt and Grabfeld formations. A likely explanation for the deviation at Breitenholz (Brhz) is the fact that the monitoring well is very close to the groundwater supply well, a situation that cannot be properly described in the coarsely gridded model.

Fig. 13

Results of approximate model calibration

Table 1 Hydraulic conductivity values estimated in approximate model calibration

The modeled groundwater flow rates support quite well the presumed hydrological conditions. The mismatch of topographic and groundwater divides (as depicted in Fig. 6) and associated gains and losses are quantified through inter-catchment flow rates (Table 2). For example, about 80% of the groundwater that is recharged in area with ID 13, an area that belongs to the topographic catchment of the Ammer, discharges to the Neckar catchment. In contrast, the Ammer catchment gains from neighboring topographic catchments of the Nagold (about 13% of the budget of the area with ID 15) and the Würm River (about 93% of the budget of the area with ID 17). Note that inter-catchment flow predominantly occurs in those areas that have been delineated based on the tracer tests (IDs 13, 15, 17, and 35), which confirms the great significance of tracer tests for the understanding and conceptualization of karstic areas also in this study.

Table 2 Water balance of catchment areas: inter-catchment, input and output flow rates (m³/day). Negative numbers indicate losses/outputs and positive number indicate gains/inputs if the table is read column-wise. x indicates no data. Refer to Fig. 6 for the location and extent of the areas

The modeled discharge rate of the Bronnbach spring is 41,580 m³/day (~480 L/s), which agrees well with the mean of the rates measured within the period 1984–2019 (530 L/s). A similar good agreement was achieved for the Ammer spring and the Herrenberg spring (modeled discharge = 11 L/s, measured ~15 L/s) as well as for the Reusten spring (7 L/s vs. 10 L/s). Furthermore, the steady-state simulation indicates that the tributaries of the Ammer, Würm, and Neckar receive groundwater only locally (Fig. S4 of the ESM).

Based on the modelled steady-state hydraulic head distribution and the derived groundwater-flow field, MODPATH was used to assess the directions of groundwater runoff from the three limestone quarries in the study area in order to localize receptors that might be at risk in case of accidental spills. The results show that the situation is not clear for any of the quarries (Fig. 14). The runoff of the two quarries near Sulz am Eck and Herrenberg-Haslach may reach the Ammer spring, the Herrenberg spring or the Ammer main stem. For the quarry in Mötzingen, the modelled pathlines identify the Bronnbach spring in the Neckar Valley as the receptor. The model reproduces the findings of the tracer test conducted in this area and confirms the previously estimated path of tracer propagation (Fig. 12). However, the model results also indicate that a runoff towards the Reusten spring and the Ammer main stem cannot be ruled out.

Fig. 14

Modelled hydraulic head distribution in Muschelkalk aquifer and predicted groundwater runoff areas from limestone quarries

Conclusions and outlook

In this study, 3-D implicit geological modeling was used to develop a consistent conceptual and geometric model of a karst aquifer system affected by tectonic features. During the model development, using implicit-modeling techniques allowed a fast, automated and straightforward generation of bounding surfaces of geological blocks when additional or modified data in terms of contact points or fault lines shall be taken into account. This is a prerequisite for exploring competing hypotheses in the iterative development and verification of the model, which may be supported through targeted field investigations. Model results must always be critically evaluated, in order to ensure geological realism.

The 3-D geological model allowed the recognition of unmapped subsurface features and provided a better understanding of the aquifer geometry. The acquired understanding of the subsurface structure and flow regime has been essential for the quantitative modeling of groundwater flow and particle-tracking-based risk assessment. Further applications of the model will support individual projects in the collaborative research center (CRC) CAMPOS (University of Tübingen 2020).

Unavoidably, there is insufficient knowledge about the inherently complex organization of hydraulically relevant features in the Middle Triassic karst aquifer. Hence, the presented groundwater model is obviously an approximation, but, thanks to the very good data bases (including an accurate and highly resolved 3-D geological model, the results of a comparatively large number of tracer tests and a reasonably large number of hydraulic head measurements), it was possible to show that the groundwater model represents the measured water regime quite well. A single-continuum equivalent porous medium model (EPM) was employed, which is the simplest approach to incorporate conduits and the one that can be most readily applied in case studies (e.g., Green et al. 2006). While this study and others (see, e.g., Worthington 2009) could show that this approach is useful to simulate karst systems in a steady-state mode to address problems involving annual or monthly average hydrologic conditions, EPM models likely cannot mimic the dynamic behavior of the karst system in simulations of, e.g., storm events, when non-laminar, i.e. turbulent flow, in conduits affects discharge (Kuniansky et al. 2016). These issues of transient simulations will need to be examined in future studies including comparisons of the EPM model with double-porosity models (i.e. dual-continuum equivalent models) and coupled continuum pipe-flow models (see Saller et al. (2013) for an example). Both model types are implemented in MODFLOW-CFP (Shoemaker et al. 2008). Note that both the double-porosity and the coupled continuum pipe-flow models require additional parameters to characterize the flow in the conduits, which will complicate model calibration. These parameters include, depending on the model, critical Reynolds numbers, mean void diameter, equivalent porosity or roughness, diameter, and tortuosity, and an exchange coefficient for the interaction between the conduits and matrix.