# Karst Water Resources in a Changing World: Review of Solute Transport Modeling Approaches

K. Ö. Çallı, G. Chiogna, D. Bittner, Vianney Sivelle, D. Labat, B. Richieri, S.

Çallı, A. Hartmann

► To cite this version:

K. Ö. Çallı, G. Chiogna, D. Bittner, Vianney Sivelle, D. Labat, et al., Karst Water Resources in a Changing World: Review of Solute Transport Modeling Approaches. Reviews of Geophysics, 2025, 63 (1), 10.1029/2023RG000811. [https://doi.org/10.1029/2023RG000811](https://doi.org/10.1029/2023RG000811)

HAL Id: hal-04948984

https://hal.science/hal-04948984v1

Submitted on 31 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Reviews of Geophysics®

## REVIEW ARTICLE

10.1029/2023RG000811

**Key Points:**

* Modeling karst transport dynamics requires a better understanding of multiscale karstic heterogeneity

* Fundamental challenges afront the successful applications of karst transport modeling are provided

* Ensuring model reliability in the real-world application of karst transport models is essential

**Correspondence to:**

A. Hartmann,
andreas.hartmann@tu-dresden.de

**Citation:**

Çallı, K. Ö., Chiogna, G., Bittner, D., Sivelle, V., Labat, D., Richieri, B., et al. (2025). Karst water resources in a changing world: Review of solute transport modeling approaches. *Reviews of Geophysics*, 63, e2023RG000811. https://doi.org/10.1029/2023RG000811

Received 25 MAR 2024

Accepted 8 DEC 2024

**Author Contributions:**

**Conceptualization:** K. Ö. Çallı, G. Chiogna, D. Bittner, V. Sivelle, D. Labat, S. S. Çallı, A. Hartmann

**Funding acquisition:** G. Chiogna, A. Hartmann

**Methodology:** K. Ö. Çallı, G. Chiogna, D. Bittner, A. Hartmann

**Project administration:** A. Hartmann

**Supervision:** G. Chiogna, A. Hartmann

**Visualization:** K. Ö. Çallı, S. S. Çallı, A. Hartmann

**Writing – original draft:** K. Ö. Çallı, G. Chiogna, V. Sivelle, D. Labat, B. Richieri, A. Hartmann

**Writing – review & editing:** K. Ö. Çallı, G. Chiogna, D. Bittner, V. Sivelle, D. Labat, B. Richieri, S. S. Çallı, A. Hartmann

<figure><img src="image_2.png" /></figure>

## Karst Water Resources in a Changing World: Review of Solute Transport Modeling Approaches

K. Ö. Çallı¹ ID, G. Chiogna²,³ ID, D. Bittner¹,⁴, V. Sivelle⁵ ID, D. Labat⁶, B. Richieri²,³, S. S. Çallı⁷ ID, and A. Hartmann¹ ID

¹Institute of Groundwater Management, TUD Dresden University of Technology, Dresden, Germany, ²Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), GeoZentrum Nordbayern, Erlangen, Germany, ³School of Engineering and Design, Technical University of Munich, Munich, Germany, ⁴Erftverband, Department of River Basin Management, Bergheim, Germany, ⁵National Research Council, Research Institute for Geo-Hydrological Protection, Perugia, Italy, ⁶Géosciences Environment Toulouse (GET), Université Toulouse 3-Géosciences Environment Toulouse-CNRS-UPS-IRD, Toulouse, France, ⁷Faculty of Engineering, Geological Engineering Department, Ankara University, Ankara, Turkey

**Abstract** Karst water resources are valuable freshwater sources for around 10% of the world's population. Nonetheless, anthropogenic impacts and global changes have seriously deteriorated karst water quality and dependent ecosystems. Multiscale karstic heterogeneity—referring to the spatial variations of the karst aquifer's physical and chemical characteristics at varying scales—is the main challenge in describing karst flow and contaminant transport dynamics. Solute transport models are powerful tools to represent and predict the spatiotemporal behaviors of contaminant migration in karst water resources. By enhancing our understanding of the transport processes, the solute transport models enable us to explore contamination risks and potential outcomes of the contamination-related issues in karst systems. Because of that, they are often used for monitoring, controlling, and managing karst water quality and dependent ecosystem functioning. This paper reviews the current state of knowledge on the modeling of karst transport processes with a focus on single-phase solute transport. By unveiling the fundamental challenges underlying a successful real-world application of karst transport models, we discuss to what extent and how we can handle these challenges. By further deriving the key challenges afront the successful modeling applications in karst systems, we, therefore, provide directions to ensure the reliable modeling of karst transport dynamics in the present context of global changes.

**Plain Language Summary** Karst groundwaters are valuable freshwater sources. However, human interventions (i.e., dam construction, groundwater abstractions) and global changes (i.e., climate change, land-use changes) have seriously deteriorated the karst water quality and dependent ecosystem services (i.e., rivers, lakes, caves). For these reasons, assessing human impacts and global changes on the karst groundwater quality is essential. Contaminant transport models are the numerical tools often used for predicting the spatiotemporal distributions of contaminants in karst groundwaters, as well as for revealing the potential outcomes of deteriorated water quality. Since the outcomes of the contaminant transport models can be used to introduce efficient adaptation strategies and regulatory decisions to solve contamination problems in karst systems, they are often utilized as valuable decision-support systems to effectively manage and sustainably protect the karst water quality and dependent ecosystems. In this review, we provide an overview of how contaminants are transported in karst groundwater, and in what ways we can resemble the movement of contamination in karst systems. By discussing the main difficulties of the successful model applications in karst groundwaters, we provide emerging directions for better predicting the contamination-related issues in the present context of global changes.

## 1. Introduction

A safe and clean groundwater supply is necessary for our society and water-dependent ecosystems (Hartmann et al., 2021; Kalhor et al., 2019). Karst groundwater is a valuable freshwater source that provides potable water for nearly 10% of the world's population (Stevanović, 2019). However, contamination issues have seriously threatened karst water quality (Henry & Suk, 2018; Reberski et al., 2022; Savio et al., 2018) and dependent ecosystems (Goldscheider, 2019; Siegel et al., 2023). For this reason, a reliable risk assessment for the contamination-related problems in karst water resources is essential to effectively manage and sustainably protect the karst water quality and dependent ecosystem functioning.

© 2025. The Author(s). This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Deterioration of karst water quality and ecosystem functioning can be linked to two main factors including anthropogenic impacts and global changes. While the anthropogenic impacts represent the direct human intervention that modifies the physical, chemical, and biogeochemical characteristics of karst systems, such as via water management alterations (i.e., groundwater abstractions, reservoir, and dam buildings), the global changes indirectly influence the karst water quality. Urbanization, climate change, soil degradation, land use, and land cover changes are such factors that predominantly impact the karst water quality and ecosystem services (i.e., rivers, springs, caves, groundwater-fed wetlands) (Gutiérrez et al., 2014; Hartmann et al., 2014; Katz, 2019).

Karst aquifers are particularly vulnerable to contamination risks due to the substantial heterogeneity manifested by the preferential flow paths (i.e., fractures, conduits, cave streams) at differing spatial scales (e.g., Bakalowicz, 2005; Ford & Williams, 2007; Goldscheider & Drew, 2014; Padilla & Vesper, 2018; Stevanović, 2018; Quinn et al., 2006, among many others). Besides these (multiple) discrete flow paths that cause the rapid spreading of the contaminants throughout the karst systems (Field, 2018; Hartmann et al., 2014, 2021; White, 2018), the close interaction with the surface waters (i.e., rivers, lakes), thin soil thickness, and rapid water infiltration to the karst groundwater substantially limit the attenuation of the contaminants released into the karst water resources (Goldscheider, 2005; Sinreich, 2014; White, 2019).

The sources of contamination in karst aquifers can vary greatly, but the origin is either geogenic or anthropogenic. Although the most common contamination source is the anthropogenic origin, geogenic contaminants (e.g., arsenic, fluoride, uranium) resulting from the dissolution of geological media (e.g., gypsum, salt) are also naturally available in karst aquifers. As a result, a wide range of contaminants (e.g., organic, inorganic, microbiological, and radioactive contaminants) can affect the karst water quality. But among them, petroleum hydrocarbons, chlorinated solvents, in particular nonaqueous phase liquid (e.g., DNAPL), toxic metals, pesticides, biological, and microbial contaminants (e.g., bacteria, viruses, protozoa) (Auckenthaler et al., 2002; Dussart-Baptista et al., 2003; Vucinic et al., 2022), microplastics (Balestra et al., 2023; Hasenmueller et al., 2023; Valentić et al., 2022), trace organic contaminants (e.g., pesticides), and emerging contaminants (e.g., pharmaceuticals and personal-care products) (Campanale et al., 2022; Selak et al., 2022) are those threatening karst water quality, human health, and dependent ecosystem functioning.

Solute transport models with differing complexity have been used to solve complex environmental issues by various scientific communities (e.g., fractured rock, geothermal energy, mining, and soil sciences). However, as the multiscale karstic heterogeneity—referring to the spatial variations of the karst aquifer's physical and chemical characteristics at varying scales—shapes the karst system's physical, chemical, and biogeochemical characteristics, the successful application of the solute transport models in karst systems is a challenging task, particularly due to the (significant) lack of knowledge on the system's physical and (biogeo)chemical characteristics that govern the karst transport dynamics. Furthermore, since the karst complexity impedes our proper understanding of the contamination-related issues in the karst systems, finding a tradeoff between the proper representation of the multiscale karstic heterogeneity in karst models and the reliability of the model predictions remains a difficult task. But beyond that, predicting the impacts of global changes on the karst water quality is still a scientific and operational challenge.

This paper reviews the current state of knowledge on the modeling of solute transport processes in karst water resources by synthesizing the existing body of knowledge in the literature. For that, we first introduce the principle of the solute transport processes in karst aquifers, and then present the state-of-the-art for karst transport modeling approaches by providing the governing transport equations that describe the process of interest with differing complexity. We further discuss the methodological concerns about the (real-world) applications of solute transport modeling in karst water resources, and hence report the fundamental challenges afront the successful modeling practice. To move forward, we give a broader perspective on the karst transport modeling, thereby providing emerging directions with potential drawbacks for the karst community.

The focus of our review is single-phase solute transport processes in karst aquifers with a particular interest in the aqueous-phase liquid contaminants. But two other groups of contaminants relevant for the transport problems in karst systems, such as petroleum products (e.g., dense non-aqueous phase liquid DNAPL (i.e., chlorinated hydrocarbons) and light non-aqueous phase liquid, LNAPL (i.e., gasoline, diesel fuel), and particle/sediment-facilitated transport (i.e., pathogens, bacteria) are not much concern of this review, because the transport mechanisms for these two groups of contaminants are rather complicated than those of dissolved solute transport as requiring additional information about the karst system's transport characteristics. For instance, predicting the
migration of DNAPL plume involves the knowledge of multiphase transport (simultaneous water, vapor migration, and solid). Therefore, the readers are referred to find more information about the contamination problem for non-aqueous phase liquids in karst water resources in Field (2018), Loop and White (2001), and Padilla and Vesper (2018).

## 2. Principles of Solute Transport Processes in Karst Water Resources

Heterogeneity in karst water resources manifests itself at multiple spatial scales (i.e., the regional scale, aquifer scale, local scale (such as borehole), and single fracture scale (e.g., Jourde & Wang, 2023). This multiscale karstic heterogeneity ultimately governs the karst flow and solute transport processes (i.e., advection, diffusion, mechanical and turbulent dispersions, and chemical reactions) (Figure 1).

The nature of heterogeneous karst aquifer mainly consists of (a) granular rock matrix, (b) small aperture fractures/ fissures, and (c) large, cavernous conduits (Worthington, 1999) (Figure 1). As a result, the karst drainage network is characterized by a complex mixed flow characteristic in which the different sets of flow equations are required for solving different flow regimes along these multiple flow paths for example, fractures, conduit, and matrix (e.g., Kresic & Panday, 2018). In general, however, the low porous matrix (also with small fractures) is considered the region where the slow flow (i.e., laminar flow regime that obeys Darcy's Law) is particularly dominant. By comparison, the larger fractures and conduits are the primary (fast) flow regions in which a transition flow varying between the laminar (Darcian) flow and turbulent (non-Darcian) flow regimes can occur (White, 2007; Worthington, 2009). Furthermore, the (hydrological) connectivity of the karst drainage network—refers to the efficiency of water transmission along different flow paths—plays a crucial role in shaping karst transport dynamics (i.e., anomalous transport or non-Fickian transport) as it causes high variability in the magnitude and direction of the karst flow.

As compared to porous or fractured groundwater systems, the most peculiar characteristic of karst aquifers is, therefore, the complex interaction between the matrix (with small fractures) and karst conduit, which often leads to the duality in the flow and transport processes (Joodi et al., 2010; Kavousi et al., 2020; Massei et al., 2003). Since this duality is often attributed to the marked difference in the flow velocities in different flow paths—typically described by the laminar flow in the matrix (with small fractures) and turbulent flow in the karst conduits (Hauns et al., 2001; Kiraly, 1998), the karst flow is expressed by the different characteristic spatial and temporal scales in describing transport processes in matrix and conduits.

The flow velocity field in karst systems is often extremely variable, that is, in the order of meters per day in the matrix (hereinafter referred to as a porous medium or fine fracture system), and of meters per second in the karst conduits (hereinafter referred to as the channels with apertures larger than a few centimeters). As a consequence, advection is the most influenced process by the flow variations in the karst systems, while the diffusive processes resulting from the random motion of the solute molecules are not significantly affected by the heterogeneity in the karst system's hydraulic properties. In comparison, the (biogeo)chemical reactions (i.e., sorption, chemical transformation) are mostly affected by the heterogeneity in the hydrogeochemical conditions within the matrix and conduits.

In the following, we describe the main transport processes occurring separately in the karst matrix and conduit, as well as the karst-specific transport dynamics based on the system's duality. We then briefly define the concepts of spreading, dilution, and mixing because there is no a general agreement on the use of this terminology in the karst literature. By defining sorption, retardation, and chemical reactions, we overview the reactive transport in karst aquifers with a focus on the carbonate dissolution/precipitation processes. Furthermore, we provide information on tracers as diagnostic tools to describe and quantify the karst transport processes, thereby only touching upon groundwater age concept and its primary role in identifying solute mixing and transport characteristics in karst systems.

### 2.1. Transport in the Matrix

There is a vast literature describing transport processes in the karst matrix. Recently, several review papers covered the topic by emphasizing the different characteristic aspects. Dentz et al. (2023) review the concepts and approaches for the quantification of passive, non-reactive solute mixing in steady uniform porous media flows across scales. Rolle and Le Borgne (2019) focus both on the description of reactive mixing, on displacement
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Conceptual representations of solute transport processes in karst water resources at differing spatial scales. (a) 3D block diagram of a karst aquifer scale, (b) aquifer scale, (c) borehole scale, (d) single-fracture scale, (e) pore-scale level (described at the Representative Elementary Volume, REV, known also as Darcy scale). Here, solute attenuation is described by the chemical transformation influenced by physical, chemical, and (biogeo)chemical processes. The figure only describes anthropogenic contamination by indicating diffuse (areal) and point sources as both are relevant to the contamination problems in karst systems.</figcaption></figure>

scenarios and continuously released plumes approaching the transport problem from the pore scale to the aquifer scale. Valocchi et al. (2018) describe mixing-limited reactions and upscaling processes from the pore scale to the continuum scale. Berkowitz et al. (2016) emphasize the complex interplay between reactive mixing and spreading, as well as the differences among the models to describe it.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Description of different spatial scales relevant for karst flow and transport mechanisms. Here, the pore-scale level is schematized at Representative Elementary Volume, REV where the individual pore geometry does not need to be solved. In the figure, the relative effect of flow and transport dynamics are prescribed, although a wide range of flow regimes (i.e., laminar, turbulent, and transition flows), transport, and reactive processes are present at all scales. Here, non-Darcian flow refers to the turbulent flow conditions where the Darcy's Law is not valid (i.e., spring discharge is not linearly proportional to the karst groundwater head drawdown). Similarly, non-Fickian transport denotes the anomalous transport (e.g., heavy tailing in solute breakthrough curve) where the advective-dispersive equation is not necessarily sufficient to solve the spatiotemporal distributions of the contamination plume (see further discussion in Section 3.2).</figcaption></figure>

In our review, despite the undoubtful relevance of the transport processes at pore scale (wherein transport pro-
cesses and chemical reactions are fundamentally governed, and thus verifiable results can be obtained from a
theoretical point of view), we will describe transport mechanisms starting from the continuum scale—generally
known as the Darcy-scale—(wherein flow, transport processes, and reactions are averaged/upscaled over a
certain Representative Elementary Volume, REV, and thus more applicable for practical applications larger than
centimeters) (Figure 2). Therefore, we refer the readers to the aforementioned review papers for the upscaling of
transport processes from the pore scale to the continuum scale.

2.1.1. Darcy Scale: Pore Diffusion and Mechanical Dispersion

Aqueous diffusion results from the random motion of the molecules of a solute in water and controls the occurrence of complete mixing at the pore scale (Rolle et al., 2012). Aqueous diffusion coefficients are compound-specific and range from 10⁻⁹ to 10⁻¹⁰ m²/s. When aqueous diffusion occurs in a porous material, it is often called pore diffusion at the Darcy scale (Bear, 1972). Since the solute is not free to move in all directions due to the presence of the solid material, pore diffusion coefficients are smaller than aqueous diffusion coefficients (Ghanbarian et al., 2013). When the (grain) Pèclet number—defined as a ratio between the time required for the particles in the fluid to move along a characteristic length by advection and the time needed for the particles's travel along the same length by diffusion (Table 1)—is smaller than 1, pore diffusion is the dominant transport process in the karst matrix.

Mechanical dispersion caused by the velocity variations at different flow pathways also contributes to the displacement of the solute of interest. Solute particles experience different velocities for different flow paths, and mechanical dispersion reflects the effect of (random) velocity variations at the pore scale during advective transport (Bear, 1972). Generally, the advective velocity at the pore scale is often laminar, and the (grain) Reynolds number (the ratio between inertial and viscous forces) is small, typically below one. If aqueous diffusion eliminates concentration gradients occurring in pore channels, mechanical dispersion at the Darcy scale
<table><caption>**Table 1**
The Formulation of the Descriptive (Dimensionless) Numbers to Represent Flow Regimes and Relative Effect of Transport Processes at Different Characteristic Length Scales (i.e., Pore-Scale Level at REV and Aquifer Scale)</caption><thead><tr><th>Scale</th><th>Parameter [dimensionless]</th><th>Mathematical definition</th><th>Equation</th></tr></thead><tbody><tr><td rowspan="2">Pore scale at REV</td><td>(Grain) Reynolds number</td><td>characteristic flow velocity [LT<sup>-1</sup>] × grain size diameter [L]<br>kinematic viscosity of water [L<sup>2</sup>T<sup>-1</sup>]</td><td>v×d/v<sub>f</sub></td></tr><tr><td>(Grain) Peclet number</td><td>groundwater flow velocity [LT<sup>-1</sup>] × grain size diameter [L]</td><td>v×d/D<sub>aq</sub></td></tr><tr><td rowspan="4">Aquifer scale</td><td>Damköhler number</td><td>aqueous diffusion coefficient [L<sup>2</sup>T<sup>-1</sup>]</td><td>k<sub>aq</sub></td></tr><tr><td>Reynolds number</td><td>local reaction rate<br>transport rate</td><td>v<sub>aq</sub>×H/v<sub>f</sub></td></tr><tr><td>Peclet number</td><td>(mean) flow velocity [LT<sup>-1</sup>] × conduit diameter [L]<br>kinematic viscosity water [L<sup>2</sup>T<sup>-1</sup>]</td><td>v<sub>h</sub>×L/D<sub>L</sub></td></tr><tr><td>Damköhler number</td><td>(mean) flow velocity [LT<sup>-1</sup>] × correlation length [L]<br>(longitudinal) hydrodynamic dispersion [L<sup>2</sup>T<sup>-1</sup>]</td><td>r<sub>c</sub>/τ</td></tr></tbody></table>

Note. The Damköhler number can be defined either based on the advective transport (i.e., reaction rate to advective transport rate) or based on the diffusive transport (i.e., reaction rate to diffusive transport rate).

is parameterized as the product of the velocity and a constant property of the porous material defined as dispersivity (Scheidegger, 1961). While dispersivity values differ in the longitudinal and transverse directions, their values are proportional to the grain diameter (Scheidegger, 1961; Xu & Eckstein, 1997). If concentration gradients at the pore scale still exist, then mechanical dispersion is compound-specific, and depends non-linearly on the grain Péclet number (Bear, 1972; Chiogna et al., 2010).

At Darcy scale, the sum of pore diffusion and mechanical dispersion is called hydrodynamic dispersion.

### 2.1.2. Aquifer Scale: Effective and Macro Dispersion

Correct characterization of the pore diffusion and mechanical dispersion is fundamental to reproducing the transport processes at the Darcy scale. Since the heterogeneity in the matrix hydraulic conductivity is hardly measurable, we are not able to assign spatially distributed values for the pore diffusion and mechanical dispersion to the karst systems. For this reason, to move from the Darcy scale to the aquifer scale (Figure 2), the concepts of effective dispersion and macro dispersion are introduced as an additional upscaling mechanism. Here, the macro dispersion is a concept for explaining “the scale effect observed in experimental data of dispersion coefficients measured in the heterogeneous porous media at different spatial scales” by Dentz et al. (2023). However, we can distinguish between the definitions of ensemble dispersion coefficients and effective dispersion coefficients. While the ensemble dispersion defines the solute distribution due to the advective spreading, effective dispersion characterizes the width of the solute plume originating from a point-like injection (Dentz et al., 2023). As a typical characteristic of both parameters, they are all scale-dependent, which means that the larger the area of interest, the larger the coefficients (e.g., Gelhar et al., 1992). Furthermore, their values are often orders of magnitude larger than the hydrodynamic dispersion coefficients.

### 2.2. Transport in the Conduits

Transport in the karst conduits shows similarities with the transport in rivers and pipes (or drains). Depending on the saturation of the karst system, the conduits can behave either as a free surface flow (under non-pressured conditions) or as a (saturated) full-pipe flow (under pressurized conditions), or as a transition flow between a free surface and full-pipe flow (under partially saturated conditions) (Liedl et al., 2003; Loop & White, 2001; Raeisi et al., 2007; Reimann et al., 2014). These variations in the conduit's saturation largely shape the solute transport mechanism. In general, solute transport in karst conduits is mainly controlled by advection and aqueous diffusion. The typical flow velocities encountered in the karst conduits are orders of magnitude larger than the ones presented in the matrix, while the turbulent flow conditions in the conduits make aqueous diffusion a dominant process along the conduit, although it is generally neglected (Hauns et al., 2001). At the aquifer scale, by comparison, the advective flow field is often not known due to the multiscale karstic heterogeneity, as well as due to the lack of knowledge about the shape of conduits and conduit network geometry.
Here, it is important to acknowledge that in preferential flow-dominated natural systems, like karst aquifers, a small hydraulic perturbation in the flow domain (i.e., storm event) can result in different flow characteristics (e.g., McDonnell & Beven, 2014), which also largely affects the karst transport processes. Celerity is such a characteristic that is different from the flow velocity (i.e., advective movement of water itself) as it mainly shapes the hydraulic response of the system, rather than defining the system's transport mechanism (McDonnell & Beven, 2014; Worthington, 2019). Celerity refers to the speed of a pressure pulse (e.g., flood wave) which is caused by a change in water level (Covington et al., 2009; Worthington & Foley, 2021). Therefore, the evolution of celerity enables a better understanding of the karst system's internal hydraulic structure, whereas the advective velocity often frames the solute transport processes (Koit et al., 2022; Worthington, 2019). Nonetheless, the measurement of the celerity of karst conduits is rarely considered in field-scale experiments, except the study by Worthington, 2019 demonstrated that celerity is often much greater than the flow velocity.

### 2.2.1. Aqueous Diffusion

Aqueous diffusion in the conduits is the same process as in the matrix. However, we can highlight three main differences here. First, in karst conduits, we do not have to reduce the aqueous diffusion coefficient by the presence of the porous medium and its tortuosity. Second, as we often move from laminar to turbulent flow conditions, the magnitude of the aqueous diffusion coefficients is much smaller than the turbulent dispersion coefficient. Third, incomplete mixing is generally not controlled by aqueous diffusion along the conduit system, but it is mainly controlled by the conduit network topology and the turbulence occurring in the conduits.

### 2.2.2. Turbulent Dispersion

The distinguishing characteristics of turbulent flow are the rapid changes in flow velocity over time and space, and the randomness of the velocity and pressure fields (Fischer et al., 1979). Therefore, the rapid irregular and chaotic fluctuations in flow velocity and pressure lead to the formation of swirling vortices, eddies (random complex patterns), and turbulent structures that interact with each other. Because of these random motions, turbulent flows occur at several characteristic scales. Reynolds number ($R_e$) is such a characteristic that characterizes the critical velocity at which a laminar flow can transform into a turbulent flow regime (Table 1). For instance, at large Reynolds number (e.g., $R_e > 500$ [-] in karst conduits), the strong nonlinearity of the equation of fluid motion results in the spreading of the kinetic energy over a range of eddy sizes, and thus causes the interaction between the large and small scales of motion (Fischer et al., 1979). These varying flow velocity fields result in the spreading and dispersion, and each of which is quantified by the turbulent dispersion coefficients. Therefore, being process scale-dependent, the definition of the dispersion coefficient needs to be done considering the specific spatial scale (Kinzelbach, 1986).

The turbulent dispersion coefficient is commonly assumed to be proportional to the flow velocity, where the proportionality constant is called dispersivity. Taylor (1953, 1954) developed a more process-based approach by correlating dispersion to the differences in flow velocity within the flow cross-section and the consequent mass transfer processes. Under turbulent flow conditions, however, the dispersion coefficient for a pipe is proportional to the pipe radius and friction velocity (Taylor, 1954). In the aquifer scale where the karst conduits exist, the dispersivity is often required to adjust during karst transport modeling by model parametrization (see Sections 3 and 4).

## 2.3. Sorption, Retardation, and Chemical Reactions in Karst Systems

Solutes that are neither affected by sorption nor undergo chemical reactions are defined as conservative (Figure 3). In our review, we overview the single-phase flow and transport processes, which are much relevant to assessing natural attenuation processes (Einsiedl et al., 2009; Flynn & Sinreich, 2010; Geyer et al., 2007; Schiperski et al., 2016; Valiente et al., 2020), although the natural attenuation is often very limited in karst aquifers (Kaçaroğlu, 1999). Therefore, we only touch upon the transport of reactive solutes with a focus on carbonate dissolution and precipitation kinetics (see Section 2.5).

Sorption refers to the process that transfers the solute concentration from the dissolved aqueous phase into an immobile phase attached to the solid geological matrix. This process is particularly relevant for solute transport in the karst matrix and at the interface between the matrix and conduit because it requires direct contact between the solute and the porous material. Petroleum hydrocarbons (i.e., DNAPL, LNAPL) and radionuclides (i.e., uranium)
<figure><img src="image_5.png" /><figcaption>b) Spatiotemporal distribution of a reactive solute</figcaption></figure>

<figure><img src="image_6.png" /></figure>

**Figure 3.** (a) Spatiotemporal distributions of a solute plume for conservative transport. (b) Spatiotemporal distributions of a solute plume for reactive transport. Here, the acronyms *Adv*, *Dis*, and *Diff* refer to advection (or advective flow), dispersion, and diffusion processes, respectively. *Sp* (sorption) and *Rc* (chemical reaction) indicate the impact of retardation and reactive processes on the movement of solute plume. In the figure, *C*₀ and *C* refer to the initial solute concentration and the concentration of solute at a given time, respectively. Here, *t*ᵢ indicates the first detection time of solute of interest (e.g., at the observation well) and *t*obs refers to the observed concentrations at the time of interest. In the figure, the two-way red arrows indicate the solute/mass exchange between the conduit and the matrix.

are such contaminants often influenced by the sorption processes in karst groundwater. Sorption can also occur on suspended particles/sediments, such as polycyclic aromatic hydrocarbons, thus leading to the so-called particle-facilitated transport (Lan et al., 2018; Massei et al., 2002; Schwarz et al., 2011). In general, sorption has two main effects on the solute concentration. First, the solute concentration in the water is reduced due to the mass transfer toward the immobile sorbed phase. Second, the advective velocity of the solute is reduced, and hence the mass exchange between the dissolved and the sorbed phase leads to an effective solute retardation (Figure 3).

Retardation is described as the delay of the solute transport relative to the groundwater flow. In karst aquifers, retardation is often assumed to be negligible due to the rapid transport in the karst conduits (Kresic & Stevanovic, 2009). By comparison, sorption reactions (adsorption or absorption) and ion-exchanges (combined sorption-desorption processes) can occur at the conduit walls/surfaces, as well as within the surrounding carbonate matrix (Field & Pinsky, 2000; Geyer et al., 2007; Morales et al., 2010). However, both processes are generally more relevant to local-scale transport problems, such as in the vicinity of an injection well (Hauns et al., 2001; Jeannin & Maréchal, 1997).

Chemical reactions can lead to the production and/or degradation (e.g., degradation of hydrocarbons such as crude oil and gasoline). A prerequisite for a chemical reaction is that all reaction partners are present in the same location at the same time while the diffusive and dispersive processes are mainly responsible for this. Nonetheless, since the diffusive and dispersive processes are slow in the karst matrix, the chemical reactions are often limited by transport mechanisms. On the contrary, reactive transport in the karst conduit is often limited by the reaction rate constant of the chemical reaction due to the high turbulent dispersion coefficient. For this reason, to reflect the relationship between the characteristic times of the chemical reaction and transport rate, the Damköhler number (*D*<sub>a</sub>) is generally used (see also Table 1). *D*<sub>a</sub> is a dimensionless (descriptive) number that compares the roles of chemical reactions and transport mechanisms at the scale of interest. *D*<sub>a</sub> is particularly useful in defining the hydro(geo)chemical processes which are functions of chemical reactions and advective-dispersive transport. For instance, *D*<sub>a</sub> is in the range of 10<sup>-3</sup>–1 for limestone fractures (Szymczak & Ladd, 2011; Wang et al., 2017), which defines the ratio of the characteristic time of chemical reaction to that of transport rate.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Schematic representations of the movement of solute plume with respect to the mass exchange rate between the (saturated) karst conduit and surrounding matrix (with small fractures). Here, the variability of mass transfer rate is controlled by the hydraulic properties of the matrix (i.e., matrix porosity) and the karst conduit (i.e., conduit diameter). For instance, due to the greater hydraulic conductivity of the karst conduits, advection and dispersion processes dominate the physical transport along the conduit whereas the matrix acts here as a sink (not source) considering the concentration gradient at the domain of interest.</figcaption></figure>

The relative effect of advection, diffusion, and reactions in karst systems at the scale of interest is often evaluated by comparing the corresponding timescales (based either on Peclet number or Damköhler numbers, or both) (see Table 1). For instance, a wormhole (an initial stage for developing a preferential flow path in karst systems due to the unstable carbonate dissolution) occurs when $D_a > 0.01$ and $10 < P_e < 1,000$ (Dreybrodt & Gabrovšek, 2018). Nonetheless, the timescales of the (biogeo)chemical reactions such as degradation, nitrification, dissolution, mineralization, and weathering are relatively slow in the matrix (i.e., seasonal, decadal, or multidecadal timescales) as compared to the advective transport in karst conduit (i.e., days, months). For that reason, the (biogeo)chemical reactions are generally overlooked in karst systems. But when considered, a (biogeo)chemical interaction term is often coupled with the advective-dispersive transport processes (e.g., Dvory, Kuznetsov, et al., 2018; Dvory, Livshitz, et al., 2018; Pastore et al., 2021) (see Section 3.2).

Overall, due to the multiscale variabilities in (hydro)geochemical conditions modified both by the karst matrix and by the composition of mobile water, the proper description of the complex reactions is almost impossible in karst systems. For this reason, sorption and (biogeo)chemical processes are crudely represented in the karst literature, while there is a general lack of knowledge of the hydrogeochemical compositions of karst systems.

## 2.4. Dynamics of Solute Transport Processes in Karst Systems

The spatiotemporal dynamics of solute transport in karst aquifers are controlled by the relative contributions of the rapid (i.e., conduits) and slow (i.e., matrix or storage zones due to the vortices and eddies along the conduits) flow components (Ewers et al., 2012; Field, 2002, 2018; Geyer et al., 2007; Hauns et al., 2001) (Figure 4). As a result, solute transport in karst aquifers is often characterized by anomalous transport which is controlled both by
the preferential flow paths (mostly by the conduits and conduit networks) and by the rate-limited solute/mass transfer between the slow and fast flow components. Overall, the key transport mechanisms in karst systems are defined as (a) physical transport along the main flow direction controlled both by advection and dispersion processes, (b) matrix diffusion and mechanical dispersion along the conduit surfaces in and out of the carbonate matrix, and (c) mass exchange between mobile water in karst conduit and matrix components, such as via sorption, ion exchange, radioactive decay, and (biogeo)chemical reactions (see also Figure 4) (Field & Leij, 2014; Field & Pinsky, 2000; Sinreich, 2011; Sinreich et al., 2007).

The nature of the transport along the conduit is strongly dependent on the flow regimes ranging from laminar to turbulent (Field, 2002; Green et al., 2006; White, 2002). The advective transport in the karst conduit are mainly influenced by the conduit's structural properties (e.g., conduit diameter, aperture, sinuosity) (Aliouache et al., 2019; Field & Pinsky, 2000; Ronayne, 2013), as well as by the local and/or regional recharge characteristics that mainly adjust the flow magnitude along the conduit and conduit networks (Doumar et al., 2018; Göppert & Goldscheider, 2008; Ji et al., 2022; Luo et al., 2024). Similarly, the rate of mass transfer at the conduit-matrix interface is shaped by the flow velocity (Bajracharya & Barry, 1997; Lehmann et al., 2022). But, the magnitude and direction of the mass exchange widely vary depending on the saturation of the karst systems, its hydraulic functioning (i.e., hydrological connectivity of the karst drainage network), and antecedent aquifer conditions (i.e., soil and/or epikarst wetness). For instance, the high velocity might reduce the contact time of the contaminant with the matrix, thus limiting the dilution of mass concentration with the lower rate of mass exchange between the conduit and matrix (Figure 4). In comparison, the mass exchange is generally more favorable under laminar flow conditions (Faulkner et al., 2009; Martin & Dean, 1999; Martin et al., 2012).

Matrix is characterized as the region where solutes can be stored for a long period, and then slowly released into the conduit and conduit network (Berglund et al., 2020; Ji et al., 2022; Nativ et al., 1999; White & White, 2005; Worthington, 2011). Since the mass transfer from the conduit to the matrix is often controlled by the (matrix) pore diffusion and mechanical dispersion, these processes can be important during the mass/solute storage in the matrix, especially for the low permeable karst media, such as chalk (Polak et al., 2002; Worthington, 2015; Zuber et al., 2011). Therefore, the rate of the mass transfer in and out of the matrix is regulated by several other factors including matrix permeability, conduit flow velocity, irregularities of conduit surface, solute retardation factor, as well as the (solute) matrix pore diffusion and mechanical dispersion coefficients (Field, 2002; Joodi et al., 2010; Kresic & Stevanovic, 2009; Luo et al., 2024; Mosthaf et al., 2018) (Figure 4).

## 2.5. Spreading, Dilution, and Reactive Mixing

A large amount of literature deals with the proper distinction among three different concepts to quantify transport processes: spreading, dilution, and mixing (e.g., Dentz et al., 2023; Kitanidis, 1994).

Spreading refers to a change of the solute plume shape due to the spatially varying velocity, thus leading to an increase in the average distance of the particles. This process is quantified by computing the second spatial moments of the solute present in the system. If we neglect aqueous diffusion, mechanical and turbulent dispersion, then spreading, if only caused by advective transport, is a reversible process and does not lead to a change in the concentration of the solute, but only to a different spatial distribution of the solute concentration.

Dilution, on the contrary, quantifies the increase in volume of the fluid occupied by the solute, which is typically used for conservative solutes. Due to dilution, the solute plume tends to occupy a larger volume leading to a decrease in the solute concentration (Kitanidis, 1994). This process is irreversible and is controlled both by diffusive and dispersive processes. For instance, a solute in karst conduit can be diluted either by variable dispersion along the conduit (passive dilution) or the (clean) water seepage from the surrounding matrix (active dilution) (Li & Loper, 2011). For that, we can often find the term mixing to refer to dilution in the literature, leading to mechanical dispersion and diffusion. Therefore, the term “mixing” is widely established for the description of conservative transport, because mixing is mainly caused by the velocity difference along the different flow paths. But “mixing” can also be used for identifying the contributions of the different water and/or solute sources in karst systems, thus often relating the source contribution with groundwater age dynamics (see Section 3.2.2).

Reactive mixing is the process that brings two different chemical species together as the result of molecular diffusion, mechanical dispersion, and turbulent dispersion (Neupauer et al., 2014), which allows a chemical reaction to occur. Here, it is important to acknowledge that dilution (mixing) is not necessarily affected by
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Dissolution and precipitation processes on the carbonate mineral surfaces for CaCO₃. (a) Equilibrium conditions at which dissolution and precipitation processes are in equilibrium, (b) Carbonate dissolution along the fractures (i.e., fracture aperture less than 1 mm), (c) (Secondary) calcite precipitation along the fractures on the CaCO₃ surface. Here, although the dissolution and precipitation are described separately, the precipitation of carbonate minerals occurs either simultaneously by the carbonate dissolution or in response to the dissolution. In the figure, *K₀* refers to the initial permeability of the carbonate matrix at the scale of interest whereas *K* is the resultant permeability after the dissolution and precipitation processes. In the formulas, H₂CO₃⁻, HCO₃⁻, H₂O, and CO₂ refer to carbonic acid, bicarbonate, water, and carbon dioxide, respectively.</figcaption></figure>

incomplete mixing at the pore scale, whereas reactive mixing is significantly controlled by incomplete mixing (Berkowitz et al., 2016), particularly in the karst matrix. Therefore, the presence of incomplete mixing very often leads to the need for case study-specific reaction rate constants to match the experimental observations.

Reactive transport couples flow, transport (i.e., advection, diffusion, dispersion), and the evolution of chemical reactions, and/or a multitude of biogeochemical reactions (Seigneur et al., 2019) (Figure 3). In karst systems, reactive transport involves a wide range of chemical reactions (e.g., acid-base and redox reactions, aqueous complexation, gas-aqueous phase exchange, and microbial reactions) and hydrogeochemical processes (e.g., dissolution of carbonate matrix, dedolomitization (calcitization of dolomite), ion exchanges). However, a prime example of reactive transport in karst systems can be the evolution of the carbonate matrix (i.e., via the dissolution process that is often induced by the water enriched with carbon dioxide) (Kaufmann et al., 2016) (see Figure 5). The dissolution and precipitation processes predominantly occur in the carbonate minerals, such as calcite (CaCO₃), dolomite (CaMg(CO₃)₂), and gypsum (CaSO₄·2H₂O), and they lead to the evolution of karst hydrogeology (i.e., development of conduit and/or conduit networks, cave genesis) and geomorphology (Kaufmann, 2009; Raines & Dewers, 1997; White, 2020)—commonly termed karstification. The dissolution/precipitation occurs along the carbonate mineral surfaces at the pore-scale level, such as at CaCO₃ surface as described in Figure 5. However, for karst genesis, both processes are considered as an initiation from a narrow single fracture (i.e., the aperture is typically less than 1 mm) or from a wormhole formation to the dissolution-enlarged karst conduits for rather larger scales ranging from a couple of meters to a few kilometers (Dreybrodt & Gabrovšek, 2018).

Given the fact that carbonate dissolution kinetics and resultant features (i.e., secondary calcite, dedolomitization) can modify hydraulic and chemical properties of karst systems at the scale of interest (Figure 5) (Birk et al., 2002, 2005; Guo et al., 2015), they can lead to the substantial changes in the matrix porosity/permeability, thereby largely shaping karst flow and transport dynamics at differing spatial scales. While the carbonate dissolution rate is associated with several factors including flow rate (e.g., laminar or turbulent), carbonate surface area, temperature, solute saturation, and the lithology where the dissolution occurs (e.g., dolomite, limestone), positive feedback between the dissolution reactions and karst flow is often available. This means that a high flow leads to higher dissolution rate along the flow paths (i.e., fracture, conduit, matrix), which result in higher porosity/permeability. In the same vein, further enlargement in the apertures of the flow paths leads to even greater flow rates (Worthington, 2009; Xiao & Jones, 2007). This (positive) feedback is, therefore, responsible for the preferential enlargements of the flow paths at larger scales (Jiang et al., 2022; Kaufmann et al., 2016). Apart from that,
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Sketch for idealized groundwater age distributions in two different karst systems. (a) Conduit-dominated karst system resembling advection-dominated transport dynamics—characterized by relatively younger water age (i.e., shorter travel time), (b) Matrix-dominated transport dynamics, which reflects diffusive-dispersive dominant transport characteristics—characterized by older groundwater age (i.e., longer residence time).</figcaption></figure>

the carbonate dissolution/precipitation control the chemical composition of the karst groundwater such that dominant water type (also referred to as hydrochemical facies) is the calcium bicarbonate (Ca-HCO₃) and calcium-magnesium bicarbonate (Ca-Mg-HCO₃). Besides, the ratio of Calcium (Ca²⁺) and Magnesium (Mg²⁺) can be used as a monitoring tracer to assess the (spatial) hydrochemical evolution of the karst water quality. For instance, the evolution of high Mg²⁺/Ca²⁺ to low Mg²⁺/Ca²⁺ can account for the dedolomitization process (e.g., Wu et al., 2009).

## 2.6. Tracers as Diagnostic Tool in Describing Transport Processes in Karst Systems

Describing flow and transport processes in karst systems is difficult mainly due to the unknown configuration and geometry of conduit network (Borghi et al., 2012; Fandel et al., 2022). Tracers are frequently used to define the karst system's flow and transport dynamics (Geyer et al., 2007; Goldscheider et al., 2008; Martín-Rodríguez et al., 2023; Mudarra et al., 2014). Various substances (i.e., solutes, particles, solids, and gasses) and physical quantities (i.e., temperature or heat) can be used for quantifying flow and transport parameters (i.e., porosity, dispersivity, flow velocity) at the spatiotemporal scales of interest (Benischke, 2021; Field, 2020; Massei et al., 2006). Additionally, tracers are frequently utilized as additional information sources to confirm or invalidate karst simulation models (Çallı, Bittner, Liu, et al., 2023; Hartmann et al., 2014, 2017; Mudarra et al., 2019; Rusjan et al., 2019; Zhang et al., 2021).

Environmental tracers and artificial tracers are two main groups commonly used to describe the karst transport processes. While the environmental tracers such as stable isotopes, radionuclides, dissolved gasses (i.e., chlorofluorocarbons and Sulfur hexafluoride), hydro-chemicals, and sediments already exist in the karst groundwater, artificial tracers (i.e., fluorescent dyes, salts, and bacteria) are introduced into the karst systems under pre-defined experimental conditions where the substance is supposed to be nonexistent under natural conditions (Benischke, 2021; Goldscheider et al., 2008).

Among environmental tracers, water isotopes such as Oxygen-18 and Deuterium are frequently used for defining the hydrological system's water age and water age distributions (i.e., travel time distributions (TTDs) and residence time distributions (RTDs) (Jasechko et al., 2017; Musgrove et al., 2023). As a terminology, water “age” refers to the time elapsed since the meteoric water entered the system (Bethke & Johnson, 2008; Long & Putnam, 2009; Suckow, 2014). In groundwater studies, water age is often denoted as residence time which indicates the time between the water entrance and leave (Aquilina et al., 2003; Maloszewski & Zuber, 1993; Maloszewski et al., 1992; Plummer, Busenberg, Drenkard, et al., 1998; Plummer, Busenberg, McConnell, et al., 1998). Whereas travel time is defined as the time required for a volume of water moving across the specified flow paths. The knowledge of water age and/or water age distributions is useful for revealing the vulnerability of karst aquifers to contamination risks (Çallı, Bittner, Liu, et al., 2023; Malík et al., 2016; Musgrove et al., 2023) as it provides valuable insight into the karst system's underlying mixing (here refers to the contributions of different sources of solute of interest), storage, and transport characteristics (i.e., piston flow, complete mixing) (Figure 6). For instance, the karst aquifers characterized by (relatively) younger water age (i.e., shorter travel time) are often at risk of contamination problems due to the system's rapid response. By comparison, the karst systems represented with (relatively) older groundwater ages (i.e., longer residence time) are often an indication of the sufficient time for the contaminant attenuation processes (Figure 6).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** A conceptual description of the (typical) anomalous behaviors of the (saturated) karst systems. (a) Solute plume migration along the different karst compartments under changing times. (b) Solute Breakthrough Curve (BTC) characterized by the dominant transport characteristics at the domain of interest within the selected time-periods given in Panel-a.</figcaption></figure>

The artificial tracer test is a powerful tool to estimate the karst system's hydraulic properties (i.e., porosity, hydraulic conductivity, conduit volumes) and transport parameters (i.e., dispersion, mass exchange rate, sorption) (Aquilanti et al., 2013; Lauber et al., 2014; Mudarra et al., 2019; Peely et al., 2021; Perrin & Luetscher, 2008). In these tests, the investigation scale at which the tracer experiment is applied is often a key factor in characterizing and quantifying the karst system's flow and transport dynamics. For this reason, the observations/data obtained by these tests are rather site-specific, and hence the estimated parameters are frequently dependent on the initial and boundary conditions of the experimental works (i.e., climatic, hydrological conditions) (Duran et al., 2016; Stevanović & Stevanović, 2021) (see Section 4.1).

The solute breakthrough curve (BTC) is the primary production of artificial tracer tests (Figure 7), which reflects solute transport behavior at the scales of interest. BTC is often used as a karst system's (characteristic) response obtained at the given spatial and/or temporal scales (Field, 2002; Labat & Mangin, 2015). However, due to the multiscale karstic heterogeneity, BTC is characterized by various anomalous characteristics such as power-law tailing, high skewness (positive or negative), and multi-modality (Figure 7), while each of which can resemble the karst system's flow and transport dynamics. In general, the tailing effects in BTC are associated with dispersion and (matrix) diffusion processes, as well as with the presence of conduit pools where the eddies and vortices form the (transient) storage zones that retain the solute of interest for a while (Bodin et al., 2022; Hauns et al., 2001; Li & Loper, 2011; Morales et al., 2010). In comparison, the multi-peaked BTC is characterized by a resemble of anastomoses (or channeled) of the conduits and/or conduit networks (i.e., auxiliary conduits) along which advective transport is often the dominant mechanism (Cen et al., 2021; Deleu et al., 2023; Dewaide et al., 2018; Goldscheider et al., 2008).

# 3. The Current Status of Solute Transport Modeling in Karst Water Resources

## 3.1. A Generic Framework for a Successful Transport Modeling in Karst Systems

To present the current status of solute transport modeling in karst water resources, we first define a workflow for a successful transport modeling application. A successful transport modeling practice starts with the definition of the karst transport problem for which some critical information about the properties of the source of nutrient or
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Schematic representations of knowledge transfer, starting from the karst system's conceptualization to the assessments of model results including model parameters and predictions. (a) System conceptualization, (b) Development and parametrization of the solute transport model (e.g., lumped transport model), and (c) Evaluation of the transport model results. Here, the dashed black line along the modeling stages indicates the feedback mechanisms across the different stages of the successful transport modeling practice, which ensures the (robust and reliable) model outcomes (i.e., model parameters and predictions). Here, the acronyms for E (Epikarst), M (Matrix), C (Conduit). *Q*(t) and *C*(t) denote karst spring discharge at a given time *t* and corresponding solute concentration, respectively.</figcaption></figure>

contaminant (i.e., source types such as point or areal sources, source locations) and the karst system's physical and chemical boundaries are mainly gathered. Similar to the development of a flow model (Beven, 2002, 2011; Wagener et al., 2001, 2021), a successful implementation of a transport model in karst systems is a strategy that contains three main successive stages, including (a) system conceptualization, (b) development of mathematical model and model parametrization, and (c) model evaluation (Figure 8). Here, each stage feeds one another with multiple feedback loops to retrieve more robust and reliable results for the successful application of a karst transport model. The overall success in this practice, therefore, reflects the level of accumulated knowledge and observational inputs for the proper descriptions of the complex physical, chemical, and (biogeo)chemical processes in karst systems.

A successful transport model application usually requires an a priori definition of a flow model that solves the water mass balance equations based on the changes of the variable of interest (i.e., hydraulic head, drawdown, or water pressure) within the flow domain where the dominant hydrodynamic conditions (i.e., regional and/or local flow directions and magnitudes) govern the karst transport processes, which has been extensively discussed in previous reviews (e.g., Hartmann et al., 2014).

### 3.1.1. Model Conceptualization for Karst Transport Processes

Perceptual understanding of karst transport processes is a key stage of successful transport modeling practice. A perceptual model is a qualitative description of the karst system's physical, chemical, and (biogeo)chemical transport characteristics, and each of which forms the definition of dominant physicochemical processes at the domain of interest (e.g., Hartmann, 2018; Hartmann, Kralik, et al., 2012; Hartmann, Wagener, et al., 2013). At this stage, the understanding of a dominant process is critical as it mainly conditions the concept of karst transport characteristics (e.g., advection-dominated transport, diffusion-dominated transport) during the conceptual model development. However, since the perceptual model is just a mental model amongst the multiple other possibilities, it is more likely to be updated by the accumulated scientific knowledge and collected data (Beven & Chappell, 2021; Wagener et al., 2021) before evolving into a conceptual model (Enemark et al., 2019; Gupta et al., 2012).
The development of a conceptual model for representing karst transport processes is a critical stage for the successful transport modeling practice (Figure 8). This model is structured based on generic scientific knowledge and a set of decisions as to how the karst system functions and/or it might function over the changing time and space (Chang et al., 2017, 2019; Hartmann et al., 2014). Model conceptualization helps to determine and/or select a suitable mathematical formulization (hereinafter referred to as mathematical model) for karst transport processes, thereby refining the governing equations within the given range of spatiotemporal scales. While these equations represent the dominant (physicochemical) processes, the mathematical formulation can be adjusted by which type of contaminant is of interest (either reactive transport or conservative transport, see also Figure 3). In most cases, the formulation for the conservative transport is expanded by coupling reactive and/or chemical processes (e.g., such as by the organic chemicals that are sorbed onto the carbonate matrix and/or fracture filling) where these processes can play a role in karst transport dynamics. In such a process-coupling, the system conceptualization can be extended by the development of a karst “hydrochemical conceptual model”. Here, a hydrochemical conceptual model refers to the spatial variations of different hydrochemical facies—or water types —by which similar geochemical characteristics are grouped (Güler & Thyne, 2004; Kresic & Stevanovic, 2009; Thyne et al., 2004). At this stage, the hydrochemical conceptual model can be seen as a complementary step for a successful transport modeling practice because it provides valuable information on modeling chemical and/or reactive transport processes, particularly during the representation of the (different) origins of contaminant sources, and their mixing characteristics at the domain of interest.

### 3.1.2. Model Parameterization in Karst Transport Models

Mathematical models translate the system conceptualization into a set of governing equations that are formulated based on the karst system's physical and (biogeo)chemical properties under (specified) boundary and initial conditions (Figure 8). The mathematical models solve solute mass balance to define the spatiotemporal dimension of the process of interest. In general, the structure of the mathematical model (also referred to as model structure) is formed and/or chosen by the modeler himself/herself depending on the adequate representation of the system and process of interest. But the selection of an appropriate model structure often depends on data availability and information collected. For this reason, different models with differing complexity can be successful for solving (different) karst transport problems at various spatial and/or temporal scales. To solve the solute mass balance, the mathematical models attribute specific values to the coefficients of the governing equations, while these coefficients are called model parameters (e.g., hydraulic conductivity, dispersivity, sorption constant). The model parameters are the main integrators of available data and information content, and hence they are often formulated based on the karst system's physical and chemical properties (see further discussion in Section 4.3).

Model parametrization is the stage for the (realistic) estimation of the parameter values to capture the spatial and/ or temporal dimensions of the karst transport dynamics. The parameter values can be obtained in three common ways: (a) direct measurements (field and/or lab-scale measurements), (b) model calibration (i.e., inverse solution, history matching), and (c) the combination thereof. Since the direct field measurements and/or laboratory scale experiments often resemble the spatial and temporal variabilities of the process of interest at the domain of interest, they are favorable for improved model parameterization due to the well-defined parameter values. Nonetheless, the characterization methods for the data acquisition and/or interpretation techniques in karst aquifers are rather appropriate for the continuum (Darcy) scale applications where the effective parameters can be readily used for the entire model domain (Faulkner et al., 2009; Huntoon, 1995; Kresic & Stevanovic, 2009; Reimann, Rehrl, et al., 2011; Zhang et al., 2021). For this reason, when the direct measurements of transport parameters (e.g., dispersivity, mass transfer coefficient) are not practical and/or their values might not be representative in resembling the process of interest, model parameters are frequently estimated via model calibration (Bittner et al., 2020; Cinkus et al., 2023; Çallı, Bittner, Liu, et al., 2023; Dewaide et al., 2016; Sivelle et al., 2023).

Model calibration refers to model tuning during which the parameter values are frequently adjusted to reduce the misfit between the observations and model simulations (Figure 8). But in successful model calibration, the measured variables (i.e., solute concentrations for the transport models, spring discharge or groundwater heads for the flow models) are required to reasonably meet the simulated counterparts (Dietrich & Uniyal, 2020; Neven & Sorab, 2021; Oehlmann et al., 2015; Sullivan et al., 2019). Such a parameterization requires a calibration metric (termed also as the objective function, performance criteria, or likelihood function) to indicate how well and to what degree the transport model represents the system's reality (Çallı, Bittner, Liu, et al., 2023; Gharari
et al., 2021; Hartmann, 2018). A single variable (i.e., spring discharge, groundwater level) is often used to quantify the misfit between the simulated and observed variables for the calibration of karst simulation models (see Figure 8c). But when it is possible, introducing multiple sources of information (i.e., isotope concentrations, water age information, gravimetric information) and/or soft data (i.e., knowledge and experience on the modeled area) into the model calibration phase is generally preferable, particularly for the improvement of model robustness and parameter reliability (Banusch et al., 2022; Çallı Bittner, & Hartmann, 2023; Çallı et al., 2022; Hartmann, Barberá, et al., 2013; Hartmann, Wagener, et al., 2013; Mazzilli et al., 2013). To a lesser extent, the information content of hydro(geo)chemical and water quality data (e.g., electrical conductivity, $NO_3^-$ and $SO_4^{2-}$) can provide distinctive chemical signatures about the karst system's flow and transport processes (e.g., Hartmann, Barberá, et al., 2013). For that reason, they are often used to reduce model uncertainty by distinguishing different states at different model compartments, for example, epikarst, conduit, and matrix (e.g., Hartmann & Baker, 2017; Hartmann et al., 2017), thus revealing the compartment-specific flow and transport characteristics. For instance, introducing $NO_3^-$ information into the model calibration marks the unsaturated zone flow and transport dynamics because it mainly resembles the (water) infiltration processes throughout the soil and epikarst (i.e., nitrogen cycling). Similarly, $SO_4^{2-}$ can constrain the saturated zone flow and transport dynamics (i.e., dissolution of evaporites) (e.g., Hartmann et al., 2017).

The number of parameters required for modeling of transport processes is typically greater than the flow models (Konikow, 1996, 2010; Reilly, 1987). As a result, multiple parameters (i.e., degrees of freedom or parameter space) are often necessary to represent the system's physical and chemical behaviors in karst transport models. The number of parameters is further increased when the multispecies and/or multiphase (e.g., air and water phases) transport is considered in the modeling framework (Dewaide et al., 2016; Faulkner et al., 2009; Loop & White, 2001; Zhang et al., 2017). For that, the model calibration for solute transport models is often framed by how the transport model is coupled with the flow model. Here, this coupling can be done at two different levels, either fully coupled (flow and transport embedded in each other) or partly/loosely coupled. Often depending on the model coupling, two approaches are frequently applied for the parametrization of karst transport models: (a) sequential calibration and (b) simultaneous calibration (Schilling et al., 2019). In sequential calibration, the (calibrated) parameters of the karst flow model are served for the parametrization of the transport model (Çallı, Bittner, Liu, et al., 2023; Husic et al., 2021; Sullivan et al., 2019), whereas the simultaneous calibration is performed for a fully coupled model, and hence whole parameter set is considered against the available observations (e.g., both karst spring discharge and $\delta^{18}O$ composition of spring discharge) (Hartmann, Kralik, et al., 2012; Hartmannet al., 2017; Kavousi et al., 2023; Zhang et al., 2020a, 2020b).

### 3.1.3. Evaluation of Karst Transport Model Results

The successful application of a karst transport model requires a robust model evaluation to ensure the reliability of the model predictions (Çallı, Bittner, Liu, et al., 2023; Hartmann, Wagener, et al., 2013) (Figure 8). In general, the tools used for the model evaluation are (a) model performance assessment (similarity of the observed variable), (b) model plausibility (alignment with the conceptual model), and (c) uncertainty estimation (Dietrich & Unyal, 2020; Hartmann, 2018; Jeannin et al., 2021; Jukić & Denić-Jukić, 2009). The model predictions might be greatly varied during the modeling practice depending on the modeler's target, system conceptualization, model selection, and applied calibration scheme (i.e., selection of objective function, calibration-validation periods, multi-objective, and/or multivariable calibration) (Bittner et al., 2020; Çallı, Bittner, Liu, et al., 2023; Mai, 2023; Mazzilli et al., 2012; Zhang et al., 2020b). However, as the karst transport model often relies on the accuracy of a karst flow model, the transport model's predictions are generally less reliable as compared to the flow models (Konikow, 1996, 2011; Scanlon et al., 2003). Besides, the incorporation of solute concentration data in the flow models often means that new model parameters are required. Since this increases the model complexity, it likely results in less identifiable parameters, as well as less confidence in the model results. Therefore, the predictions of transport models are always uncertain to some degree, thus requiring diagnosing and quantifying sources of uncertainty in the model outcomes.

Uncertainty propagates and amplifies throughout the modeling practice (Figure 8). For that reason, the exploration of the principal sources of uncertainty is strictly limited (Beven, 2005; Hartmann, 2018). Similarly, the translation of each uncertainty source into a quantitative form—refers to the attribution of a numerical value for each uncertainty source—might change from one modeler to another and is mostly dependent on the progressive decisions and choices over the modeling practice, from the system conceptualization to the model assessment as
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Generic classification of karst simulation models based on the model parametrization considering process complexity and data requirement. (a) Conceptualization of the karst aquifer physical boundaries depicted by the gray-shaded area with a blue-indicated karst conduit and conduit network (the blank circles also describe the swallets/sinkholes along the conduit network), (b) Spatially lumped karst simulation models depicted based on the solute concentration distribution over different karst compartments including epikarst, matrix, and matrix, (c) Spatially distributed karst simulation models described considering the spatial distribution of the solute concentration. The classification is adapted from Hartmann et al., 2014 (see further details in Hartmann et al., 2014). Herein, an artificial tracer test is described only to demonstrate the spatial distributions of contamination plume across two main karst simulation approaches.</figcaption></figure>

described by Figure 8. Therefore, different sources of uncertainty in model results, particularly data uncertainty, model structure, and parameter uncertainties could be reduced to a certain degree, only if a systematic model assessment of the model results is carried out. The model assessment can identify what we need for the reduction of model uncertainty, and how much we shall put effort to reduce the different sources of uncertainty. Although this is not an easy task, some techniques such as Generalized Likelihood Uncertainty Estimation (GLUE) (Beven & Binley, 1992), Bayesian model averaging (BMA) (Hoeting et al., 1999; Neuman, 2003), and multi-model ensemble method (MME) can be useful for the systematic and robust uncertainty assessment.

## 3.2. Description and Classification of Solute Transport Models in Karst Water Resources

### 3.2.1. Generic Classification of Karst Simulation Models

To deal with a wide variety of flow and transport problems in karst water resources, two fundamental approaches have been used for representing the physical heterogeneities of karst systems: (a) spatially lumped karst simulation models and (b) spatially distributed karst simulation models (Hartmann et al., 2014) (Figure 9), while various forms of models in between (called semi-distributed or hybrid models) are available for the simulations of the karst systems (Gill et al., 2013, 2021; Hartmann, Barberá, et al., 2013; Kovács & Sauter, 2014). Since the distributed models represent the physical processes more explicitly based on the physical laws by discretizing the karst systems into several grid cells (Figure 9c), the data requirement is often too high to adequately represent the karst flow and transport processes at varying spatial scales. By comparison, the spatially lumped karst simulation models do not explicitly consider the karst system's spatial variabilities (i.e., geometric boundaries, flow pathways), thereby implicitly representing the physical processes across the domain of interest (for further discussions, see Hartmann et al., 2014).
For karst systems, the spatially distributed karst simulation models are classified as (a) Equivalent Porous Medium Approach (EPM), (b) Double Continuum Approach (DC), and (c) Combined Discrete-Continuum approach (CDC) (Ghasemizadeh et al., 2012; Hartmann et al., 2014). EPM is a continuum-based modeling approach that considers the Representative Elementary Volume (REV) and groundwater flow equations (Hartmann et al., 2014; Scanlon et al., 2003) (Figure 9c). REV is defined as the smallest volume over which the averaged flow and transport parameters can be represented by the spatial scale where we can still apply the Darcy law based on the laminar flow conditions. REV represents a threshold volume in which the hydrological properties are assumed to be invariant after then (Medici & West, 2021). Nonetheless, due to the presence of the karst conduit and/or conduit networks, the flow in karst conduits is often characterized by turbulent flow characteristics, and hence representing the whole karst system with a single REV is not necessarily possible (Geyer et al., 2013; Kiraly, 2002). For that reason, EPM is often (more) suitable for simulating larger-scale model applications where the karst flow and transport processes can be described by a single REV (Anderson et al., 2015; Scanlon et al., 2003). In DC approach, the karst system consists of (at least) two overlapping continua (i.e., karst conduit and matrix) where the physical and/or chemical interaction between these two continua is addressed by an exchange rate (i.e., flow exchange and/or mass exchange) (Cherubini, 2008; Kiraly, 1998; Maréchal et al., 2008; Robineau et al., 2018; Teutsch & Sauter, 1998). Here, the rate of exchange is proportional to the differences in hydraulic head and/or concentration gradient between the karst conduit and matrix (Kavousi et al., 2020; Oehlmann et al., 2015; Reimann, Geyer, et al., 2011; Reimann, Rehrl, et al., 2011; Reimann et al., 2014).

By comparison, the CDC treats the karst conduits as the main flow and transport regions while an additional parameter set (e.g., mass exchange parameter) is defined for the matrix component. According to the CDC approach, karst conduits can be either treated as high permeable regions or characterized as drains or pipes (Figure 9c). This approach considers the discrete conduit and/or conduit networks, and it is often more reasonable choice for the well-defined sites and/or small-scale karst modeling applications, due to extensive data need and high computational effort for representing karst conduits and their topology (e.g., conduit networks, fracture intersection) (Kovács, 2003; Ravbar et al., 2011; Xanke et al., 2016). Here, the turbulent flow characteristics in karst conduits and/or discrete conduits are often simulated by Conduit Flow Process (CFP) (developed by Shoemaker et al., 2008 and later extended by Reimann et al., 2013, 2018, 2023). CFP has been successfully applied for well-characterized karst systems to simulate conduit flow and transport characteristics (Chang et al., 2019; Shirafkan et al., 2023; Spiessl et al., 2007; Sullivan et al., 2019; Xu et al., 2015). Therefore, CFP can consider both laminar and turbulent flow regimes by coupling solute and/or heat transport processes (Kavousi et al., 2020; Reimann et al., 2018).

Spatially lumped karst simulation models (i.e., bucket-type models, cascading reservoirs) conceptualize the entire karst systems based on the karst system's compartments (i.e., matrix, conduit, epikarst) in which average physical properties of the karst system are represented (Charlier et al., 2012; Geyer et al., 2008; Jukić & Denić-Jukić, 2009; Mazzilli et al., 2019; Orban et al., 2010; Tritz et al., 2011) (Figure 9b). By quantifying the vertical and/or horizontal water exchange between these compartments, the lumped transport models are often coupled with the “mixing” concept that distinguishes the different contributions of solute concentration from the different compartments (i.e., epikarst, matrix, conduit).

### 3.2.2. Classification of Solute Transport Models in Karst Systems

As there is not a generic classification for the solute transport models in karst water resources, we suggest the following classification for modeling approaches of solute transport processes in karst systems. This proposed classification is categorized based on the complexity of karst system conceptualization, (model) data requirement, and (model) numerical formulation. Therefore, the solute transport models are divided into two main classes: Lumped transport models (LTMs) and Advanced transport models (ATMs) (Figure 10).

Overall, LTMs consider the karst flow field as a single non-distributed black box by assuming the karst flow is constant at the domain of interest (Bayari, 2002; Ozyurt & Bayari, 2005, 2008). Since LTMs do not require specific knowledge about the various transport processes (e.g., advection, dispersion, mass exchange) within the flow domain, they are often preferable for karst modeling applications, in particular at the aquifer and regional scales (Maloszewski & Zuber, 1996; Toth & Katz, 2006). ATMs, by comparison, explicitly consider the impact of the karst heterogeneities both on the flow velocity field and on the anomalous transport (e.g., heavy tailing, and early arrivals in BTC). Therefore, the underlying assumptions of ATMs are physically better than LTMs, whereas
<figure><img src="image_13.png" /><figcaption>**Figure 10.** The conceptual descriptions of the solute transport models often used in karst systems. (a) Lumped transport models (LTMs) including complete instantaneous mixing model (CIMM) and lumped parameter model (LPM), (b) Advanced transport models (ATM) such as Multi-Rate Mass-Transfer Model (MRMT) that resemble the solute transport along one mobile region (i.e., conduit, cave streams) interacting with the multiple immobile regions (i.e., matrix, eddies, vortices). Here, *C<sub>out</sub>* refers to solute concentration at time *t* at the outlet of the system. The two-way arrow in MRMT indicates the mass exchange between mobile and immobile zones whereas the single-way arrows demonstrate the diffusion between the immobile zones. Here *c* and *m* refer to conduit and matrix, respectively.</figcaption></figure>

the mathematical solutions are more complex and rather prone to numerical errors, in particular the numerical dispersion caused by truncation errors, and computational instabilities.

### 3.2.2.1. Lumped Transport Models (LTMs)

Lumped transport models (LTMs) comprise the compartmental mixing models—often referred to as complete and instantaneous mixing models (CIMMs)—and lumped parameter models (LPMs) (Figure 10a).

CIMM is a compartmental model approach in which a solute of interest is mixed in a the model compartment (usually the bucket of a lumped model) based on an assumption of instantaneous and complete mixing. The instantaneous and complete mixing can be applied to each of model compartments (i.e., epikarst, matrix, and conduits as depicted in Figure 9b) or to the entire simulated dynamic water storage. Besides the complete-mixing assumption, partial mixing within the same storage has been also applied in some studies (e.g., Fenicia et al., 2010; McMillan et al., 2012). In CIMMs, the structure of the lumped karst flow model specifies the flow paths between the compartments, whereas the solute mass balance between those compartments is solved by the solute mixing ratios. Although simple, this approach brings the advantage of requiring no (complete mixing, Hartmann, Lange, et al., 2012; Hartmann, Wagener, et al., 2013; Hartmann et al., 2016) or few (partial mixing, Zhang et al., 2020a) parameters for the transport models. However, CIMMs ignore the advection and time-variant mixing and transport characteristics of the karst systems. Therefore, their application for karst transport modeling is very often limited.

LPMs were a common modeling approach in karst aquifers mainly due to their simplicity (significant simplifications with fewer parameters) and suitability for the data-rare sites (e.g., Bayari, 2002; Małoszewski et al., 1992, 2002; Zuber, 1986; Zuber & Małoszewski, 2001, among many others). LPM was first introduced by Małoszewski and Zuber (1982) to solve the transport problems in karst systems.
<table><caption>**Table 2**
Different Approaches Depending on the Mixing Assumption for LPMs</caption><thead><tr><th>Approach</th><th>Mixing assumption</th><th>Response function, g(τ)</th><th>Parameters</th></tr></thead><tbody><tr><td>Piston Flow Model (PFM)</td><td>No mixing</td><td>g(τ) = δ(t - τ)</td><td>τ</td></tr><tr><td>Exponential Model (EM)</td><td>Complete mixing (known also perfect mixing or well-mixed)</td><td>g(τ) = τ<sup>-1</sup> exp<sup>-t/τ</sup></td><td>τ</td></tr><tr><td>The combination of PFM and EM:</td><td></td><td></td><td></td></tr><tr><td>(1) Exponential-Piston Flow Models (EPM):</td><td>η = 1 → perfect mixing<br>η → ∞ → piston flow model</td><td>g(τ) = {<br>τ/η<sup>-1</sup> exp(-ηt/τ + η - 1), t ≥ τ(1 - η<sup>-1</sup>)<br>0, t < τ(1 - η<sup>-1</sup>)</td><td>τ and η</td></tr><tr><td>- Perfect-Partial-Piston flow model</td><td>η = V<sub>τ</sub>/V<sub>m</sub> = (V<sub>m</sub> + V<sub>p</sub>)/V<sub>m</sub></td><td></td><td></td></tr><tr><td>- Partial-Piston Flow model</td><td></td><td></td><td></td></tr><tr><td>(2) Partial-Mixing model</td><td>V<sub>m</sub> : perfect mixing reservoir volume</td><td></td><td></td></tr><tr><td>(3) Partial Exponential Model (PEM): a derivative of the exponential model (EM)</td><td>V<sub>p</sub> : piston flow reservoir volume</td><td></td><td></td></tr><tr><td>Dispersion Model (DM)</td><td>Partial mixing based on the Peclet number (P<sub>e</sub>) by defining relative magnitude of advection and diffusion</td><td>g(τ) = <sup>1</sup>⁄<sub>√(4πPe<sup>-1</sup>(t/τ))</sub> exp[-(1 - (t/τ)<sup>2</sup>)/(4Pe<sup>-1</sup>(t/τ))]</td><td>τ and 1/P<sub>e</sub></td></tr><tr><td>Gamma Distribution (GM) (Pearson Type 3)</td><td>No assumption, depends on the shape parameters (α and β)</td><td>g(τ) = τ<sup>α-1</sup>/β<sup>α-1</sup>(α) e<sup>-α/β</sup> = <sup>α</sup>⁄<sub>(τ/α)<sup>α</sup>Γ(α)</sub> e<sup>-α/β</sup><br>or g(τ) = τ<sup>β-1</sup>/β<sup>β-1</sup>(α) e<sup>-β/τ</sup><br>Γ(α) : gamma function</td><td>α, β,<br>and<br>τ = α, β</td></tr><tr><td></td><td></td><td>α → 1 : refers to the exponential distribution</td><td></td></tr></tbody></table>

Note. Here, τ [T] denotes MTT, η [-] is percentage of the contribution of the different flow regions (e.g., conduit, matrix), α [-] and β [-] are the shape and scale parameters, respectively.

LPM is a groundwater age-based modeling approach that solves age distributions (e.g., TTDs, RTDs) of solute of interest at the karst system outlet (i.e., karst spring, observation well). In doing so, LPM converts the input solute concentrations (i.e., δ¹⁸O concentrations of atmospheric precipitation) to the output solute concentrations (i.e., δ¹⁸O concentration of karst spring) using a convolution integral in Equation 1.

$$C_{out}(t) = \int_{0}^{\infty} C_{in}(t - \tau) g(\tau) \exp(-\eta\tau) d\tau \quad (1)$$

where $C_{out}$ [ML⁻³] and $C_{in}$ [ML⁻³] are the input and output solute/tracer concentrations, respectively. τ [T] is the travel time (or residence time is often used in groundwater studies). g(τ) is the impulse response function (known also as transfer function) that is characterized by the system's TTDs. Here, t refers to the calendar (or current) time and η is a decay constant accounting for the non-conservative transport mechanisms. In Equation 1, when a conservative tracer is considered for the model application, the decay term is then set to zero (Małoszewski & Zuber, 1982; Zuber & Małoszewski, 2001). But the critical parameter for karst transport modeling in Equation 1 is the mean transit time (MTT) of the solute of interest, which is obtained for the time-invariant hydrological conditions as followed by,

$$MTT = \int_{0}^{\infty} tg(\tau) d\tau \quad (2)$$

In LPMs, g(τ) plays an important role in quantifying the physical spreading of solutes of interest at the transport domain. Since it describes how different (groundwater) ages in the convolution integral contribute to water age distributions at the karst system's outlet (i.e., springs, observation well), a wide range of g(τ) can be defined (see also Table 2). For that, a proper definition or selection of a suitable response function for the LPMs is a critical
<figure><img src="image_14.png" /><figcaption>b) Multi-dispersion Model (MDM)</figcaption></figure>

<figure><img src="image_15.png" /></figure>

**Figure 11.** The conceptual descriptions of (a) the single-fissure dispersion model (SFDM), and (b) the multi-dispersion model (MDM). Here, (c) refers to a single conduit while $c_1$ and $c_2$ are denoted as the solute concentration in Conduit-1 and Conduit-2. $C_{out}$ refers to the solute concentration at time $t$ at the outlet of the system. Here, the single-way arrows demonstrate the (matrix) diffusion to the conduit walls and/or karst matrix.

step in representing the underlying flow and transport processes (Małoszewski & Zuber, 1993; Zuber, 1986). This functional form can either be defined by a priori understanding or by comparing the simulations of different response functions to the historical records of solute concentrations at a particular sampling location (e.g., springs, observation wells). For instance, the simplistic LPMs resemble the individual flow and transport paths with no mixing in karst systems (i.e., for the conduit-dominated karst aquifers) or well-mixed karst systems (i.e., for the matrix-dominated karst aquifers) (Table 2). Here, while the former would be represented by the Piston Flow Model (PFM), the latter would be represented by the Exponential Model (EM) (Małoszewski & Zuber, 1982, 1993; Zuber, 1986). The different combinations of the response functions can also be used to represent the karst transport dynamics, such as the combined exponential piston models (EPM) and dispersion models (DM) which resemble the solute mixing that varies between perfect mixing and no mixing (Table 2) (Małoszewski & Zuber, 1982, 2002; Maloszewski et al., 2002; Zuber & Maloszewski, 2001).

Furthermore, as $g(\tau)$ in Equation 1 defines groundwater age distributions in the karst systems (e.g., TTDs or RTDs), it also enables us to link between the LPMs and groundwater age models—known also as Direct Age Models (DAMs)—(Cornaton & Perrochet, 2006a, 2006b; Ginn, 1999; Ginn et al., 2009; Goode, 1996; Sanford, 2011; Suckow, 2014). These models predict the water age distributions by solving the Advection-Dispersion Equation in which the “solute mass” concentration is replaced by a groundwater “age term” (see more details in Ginn, 1999, Ginn et al., 2009; Goode, 1996). Nonetheless, solving the age (master) equation is often challenging as it requires specific numerical schemes to calculate the variations in the water age within the reservoir's storage, thereby involving the consideration of the effect of aging in the storage. However, this is not entirely similar to solving solute concentration because the water-age information is a derived-information content from the concentration of interest, such as from the $\delta^{18}$O concentration of spring discharge (e.g., Çallı et al., 2023a; Rusjan et al., 2019; Zhang et al., 2019).

Single-fissure dispersion model (SFDM) and Multi-dispersion model (MDM) (known also as the multi-flow dispersion model) are two common types used for estimating karst flow and transport parameters (Figure 11) by including more process understanding compared to simple LPMs. These models differ from LTMs in that additional parameters (i.e., dispersion parameter) are required for solving dispersive-diffusive processes in the karst matrix. In general, as SDFM and MDM analytically solve the karst transport problems, they are more suitable for the well-defined aquifer conditions where the initial and boundary conditions are known (Frank et al., 2021; Göppert & Goldscheider, 2008; Kübeck et al., 2013; Vincenzi et al., 2009). These approaches are derived from the earliest analytical solutions of one-dimensional ADE by Lenda and Zuber (1970) and Kreft and Zuber (1978), and the general form of the advection-dispersion model (ADM) is as follows by;
$$C_f(x,t) = \frac{M}{Q \cdot t_0 \sqrt{4\pi \cdot P \cdot e^{-1} \cdot \left(\frac{t}{t_0}\right)^3}} \exp\left[\frac{\left(1 - \frac{t}{t_0}\right)^2}{4 \cdot P \cdot e^{-1} \cdot \frac{t}{t_0}}\right] \quad (3a)$$

$$P_e = \frac{v,x}{D_L}; \quad D_L = \alpha_L \cdot x \quad (3b)$$

where M [M] is the solute mass, $C_f(x,t)$ [ML⁻³] is the concentration of solute of interest at the longitudinal transport distance, x [L] when the time is t [T]. v [MT⁻¹] refers to the mean flow velocity. Q [M³T⁻¹] is discharge, and t₀ [T] is the time of mean flow (or refers also to the mean residence time) [T]. D_L [L²T⁻¹] is the longitudinal dispersion, and α_L [L] is the longitudinal dispersivity. 1/P_e is called as a dispersion parameter for describing the (matrix) diffusion effect on the solute transport behavior (Małoszewski & Zuber, 1982, 1992; Zuber, 1986).

In Equation 3a, P_e refers to the (dimensionless) Péclet number defined at the aquifer-scale, which scales the rate of advection to diffusion, thus formulated by the characteristic flow length multiplied by the mean flow velocity (θ) and divided by the diffusivity (Labat & Mangin, 2015) (see also Table 1). In Equation 3a, P_e represents the ratio of the advective transport to the diffusive transport. Along the characteristic length, P_e is assumed to be constant and homogeneous, although, this rarely holds the reality in karst aquifers. For a aquifer scale, often under laminar flow conditions, when P_e is below 0.4 [-], diffusion controls the solute transport. By comparison, P_e > 6.0 [-] indicates the advection-dominated transport characteristics, while 0.4 < P_e < 6.0 demonstrates the (transition) transport characteristics varying between the diffusion-dominated and advection-dominated transport (Field & Nash, 1997). But, when P_e exceeds 1,000–2,000 [-], the flow in the conduit becomes turbulent (Kaufmann et al., 2014).

SFDM or Parallel fissure dispersion model (PFDM) by Małoszewski and Zuber (1985) is based on the theory of the parallel-plate modeling approach introduced by Sudicky and Frind (1982). This approach considers the solute transport along the parallel fissure plates and into the matrix coupled by the matrix diffusion (Małoszewski & Zuber, 1985). For karst aquifers, SFDM is first applied to the fissured chalk aquifers by Małoszewski and Zuber (1990) and modified later by Witthüser et al. (2003) to simulate transversal dispersion in chalk aquifers. To solve the karst transport problem, SFDM requires three parameters including the dispersion parameter, diffusion parameter, and mean transit time (Małoszewski & Zuber, 1985; Małoszewski et al., 1994; Witthüser et al., 2003).

MDM is an extension of ADM, and it considers the multiple discrete flow and transport pathways. While MDM solves ADM in each flow path based on the specified parameters (i.e., flow rate, transit time, and dispersivity) (Małoszewski et al., 1992, 1994), it is often used for (a large-scale) multi-peaked BTC analysis in karst systems (Goldscheider et al., 2008; Richter et al., 2022).

### 3.2.2.2. Advanced Transport Models (ATMs)

The anomalous transport is ubiquitous in karst systems and frequently characterized by the rapid spread of the font of the contaminant plume due to the existence of preferential paths, as well as by the slow diffusion caused by the mass exchange between the karst conduit and surrounding carbonate matrix (Berkowitz et al., 2006; Ewers et al., 2012; Goeppert et al., 2020; Wang et al., 2022). Diffusion and retardation processes often result in the late-time heavy tailing in BTC, and therefore refer to time-nonlocality in karst transport behavior. By comparison, highly skewed and/or multi-peaked BTC characteristics are defined by the space-nonlocality in solute transport behavior (see also Figure 7) (Haggerty & Gorelick, 1995; Haggerty et al., 2000; Zhang et al., 2009, 2020). In this context, ATMs intend to capture the spatiotemporal dynamics of the anomalous transport in karst systems.

In this review, ATMs are divided into three main modeling concepts: (a) Continuous-time Random Walk (CTRW) that discretizes the solute concentration into the particles to describe the spatiotemporal evolution of each particle, (b) Fractional Advection-Dispersion Equation (FADE) wherein the solute movement obeys the non-Gaussian distribution, and (c) Multi-rate mass transfer (MRMT) models that divide flow and transport domain into (at least) matrix and conduits. Broadly speaking, MRMT and FADE represent the Eulerian-type (Continuum) modeling approaches, while CTRW is the Lagrangian-based transport modeling approach. Here, the Eulerian-type models consider the volume-based flow and transport by discretizing the problem domain into
the different volumes, whereas the Lagrangian-based approaches, such as Particle Tracking (PT), Random Walk Particle Tracking (RWPT) discretize the solute concentration into the particles, thereby considering the solute particles as a discrete phase without an explicit need for a flow field definition. CTRW and FADEs are robust transport modeling approaches, especially for capturing the time-nonlocalities in solute transport behaviors in heterogeneous groundwater systems (Benson et al., 2000; Meerschaert et al., 1999; Metzler & Klafter, 2000; Schumer et al., 2009) although the applications of both modeling approaches in solving karst transport problems are rather new. In the following, we will provide fundamental knowledge about these advanced karst transport models and their primer applications in modeling of karst solute transport processes.

The fundamental governing equation for one-dimensional contaminant transport with a source term in karst groundwater is traditionally modeled by the advection-dispersion equation (ADE). ADE assumes that dispersive transport is Fickian, meaning that dispersion is linearly dependent on the groundwater flow velocity. Therefore, ADE (with a reactive source term) is formulated as follows by:

$$R \frac{\partial C}{\partial t} = -\nu \frac{\partial C}{\partial x} + D \frac{\partial^2 C}{\partial x^2} + rC \quad (4)$$

in which $C$ [ML⁻³] is the solute concentration in time (t) [T]. $ν$ [LT⁻¹] is the flow velocity and $x$ [L] is the (longitudinal) transport distance in the main flow direction. $D$ [L²T⁻¹] is the (longitudinal) hydrodynamic dispersion coefficient along the main flow direction. $rC$ accounts for the gain or loss contributions of the chemical reactions defined by a reaction rate (r). For conservative single dissolved solute transport, $rC$ is set to 0. In Equation 4, $R$ [-] is denoted as a retardation factor and estimated by the following formula,

$$R = 1 + \rho_s K_d \quad (5)$$

where $K_d$ [L³M⁻¹] is denoted as the distribution (or partition) coefficient that is described as a ratio of the sorbed concentration into (geological) material and the dissolved concentration in groundwater at (thermodynamically) equilibrium conditions. In karst aquifers, $K_d$ is often used to describe the sorption degree of any species in karst groundwater based on the distribution of solute between solid matrix and groundwater (e.g., Pastore et al., 2021). Therefore, when $K_d$ is equivalent to 0, it equals 1, meaning that $R$ [-] is negligible during the transport at the scale of interest. To account for sorption and retardation processes in karst systems, the sorption isotherm models (i.e., linear, Langmuir, Freundlich isotherms) that mainly ignore the complex sorption kinetics are often preferred (e.g., Frank et al., 2021; Hillebrand et al., 2015; Luhmann et al., 2015; Morales et al., 2010; Schiperski et al., 2016; Tinet et al., 2019, among others).

Among the Advance Transport Models (ATMs), FADE uses an additional “fractional” term in the classical ADE given in Equation 4. The “fractional” term in FADE refers to the fractional derivatives of dispersion that are incorporated by time and/or space for solving time and/or space non-localities in BTC (see for latest discussions on FADEs by Sun et al. (2020)). For one-dimensional conservative solute transport form of FADE is followed by Meerschaert et al. (1999) and Kelly and Meerschaert (2019);

$$\frac{\partial C}{\partial t} = -\nu \frac{\partial C}{\partial x} + D \frac{\partial^\alpha C}{\partial x^\alpha} + D(1 - \rho) \frac{\partial^\alpha C}{\partial (-x)^\alpha} \quad (6)$$

where $C$ [ML⁻³] is the solute concentration in time (t) [T] and at the transport distance (x) [L] in main flow direction. $ν$ [LT⁻¹] is the mean flow velocity, and $D$ [L²T⁻¹] is a (longitudinal) hydrodynamic dispersion coefficient in the main flow direction. Here, $0 ≤ ρ ≤ 1$ describes the skewness of the transport processes while $α$ [-] defines the order of fractional derivatives.

Depending on which nonlocality problem (time or space nonlocalities) is being solved, FADE is either called the time-in-fractional advection-dispersion equation or the space-in-fractional advection-dispersion equation. Here, the time-nonlocality is typically related with the diffusion and retardation processes resulting in late-time heavy tailing in breakthrough curve (BTC), while the space-nonlocality leads to highly skewed and/or multi-modality on the BTC (e.g., Feehley et al., 2000; Haggerty & Gorelick, 1995; Haggerty et al., 2000; Zhang et al., 2020, among others). For solving karst transport problems, Yin et al. (2022) proposed the efficient fractional-in-time
transient storage model for simulating the multi-peaked BTC in karst aquifers which considers the karst solute transport in two domains by solving the tailing effects in BTC due to the mass exchange between these regions. Besides, Zhang et al. (2024) introduced an improved version of the FADE to solve the late-time behaviors of BTC characterized by the multiple-conduit flow within an intermediate-scale karst aquifer, named as the impulsive tempered FADE approach.

CTRW considers the solute (particle) transport as a random walk in spatial and temporal scales (Berkowitz et al., 2006; Dentz & Berkowitz, 2003). It is a generalization of the Random walk methods in that the spatial displacement and temporal resolution (or time-step) are counted as random variables (Dentz et al., 2004). The formalism of CTRW originates from Einstein's Brownian motion approach which considers the random motions of the solute particles correlated with time or space. While CTRW has been extensively applied to the fractured rocks and aquifers (Berkowitz et al., 2001, 2006; Hansen & Berkowitz, 2020), the first application of CTRW in karst systems was done by Goeppert et al. (2020) for the simulation of the flow and transport processes in an alpine karst aquifer. CTRW is particularly useful for simulating the late-time and heavy-tailed characteristics of the solute transport behavior because it allows dispersivity to change in time or in space (Berkowitz et al., 2006; Berkowitz & Scher, 1995; Dentz & Berkowitz, 2003). For this reason, CTRW better resembles the dispersion and diffusion processes as compared to the particle tracking techniques (Tsang & Tsang, 2001). The applications of the Particle-tracking methods in karst aquifers are generally preferable for the delineation of capture zone characteristics in the well-head protection studies (e.g., Klaas et al., 2017; Qiao et al., 2015; Rayne et al., 2001).

The CTRW framework is defined in the Laplace domain wherein many complex integrals and differential operations to represent the system's non-linearity are converted into the linear operations that are easy to solve, as followed by (based on the formulation by Zhang et al., 2020):

$$uC(x, u) - C_0(x) = -M(u) \left[ v \frac{\partial C(x, u)}{\partial x} - D \frac{\partial^2 C(x, u)}{\partial x^2} \right] \quad (7a)$$

where the time, $t$ [T] → $u$ [T] in the Laplace domain. Here, $v$ [LT⁻¹] and $D$ [L²T⁻¹] refer to the effective velocity and effective dispersion coefficient, respectively. $M(u)$ denotes the memory function expressed by Equation 7b:

$$M(u) = t_1 u \frac{\psi(u)}{1 - \psi(u)} \quad (7b)$$

in which $t_1$ is the (mean) transient time (or mean waiting time). $\psi(u)$ is the probability density function (PDF) of the solute particles' waiting time in the Laplace domain (for further details see Berkowitz et al., 2006; Goeppert et al., 2020; Zhang et al., 2020).

MRMT models are the generalization of the mobile-immobile model (MIM) concept—known also as the classical dual-domain approach—(Carrera et al., 1998; Coats & Smith, 1964; Gerke & Van Genuchten, 1993; Haggerty & Gorelick, 1995). MIM simulates the flow and solute transport processes based on two distinct flow and transport regions denoted as mobile zone (e.g., conduit, fracture, fissure) and immobile zone (e.g., matrix, voids, vugs, dead-end passages along conduits and/or cave streams) (see Figure 12). Here, MRMT differs from the MIM approach in that it considers a single mobile domain interacting with multiple immobile regions with differing physical and chemical characteristics (see also Figure 10b). For this reason, solving MRMT requires multi-rate mass transfer exchange terms, each of which describes its interaction zone with the mobile region such as karst conduit.

In karst aquifers, adaptations of the MIM approach have been extensively used to simulate karst transport dynamics. In the MIM approach, the dominant transport is described by the advective-dispersive processes along the mobile zone while the mass transfer between the mobile and immobile regions is only controlled by the (matrix) diffusion without any fluid exchange (Figure 12) (Becker & Bellin, 2013; Maloszewski & Zuber, 1990; Toride et al., 1995). Therefore, solute transport along the karst conduit is set by the classical ADE with an additional term that describes the mass transfer rate between the conduit and matrix. For one-dimensional solute transport with the MIM approach is then written as follows:
<figure><img src="image_16.png" /><figcaption>**Figure 12.** The schematic description of the Mobile-Immobile (MIM) Models. Here, $C_{out}$ refers to solute concentration at time $t$ at the system's outlet (e.g., karst spring, observation well). In the figure, the two-way arrows demonstrate the (matrix) diffusion to the conduit walls and/or to the karst matrix.</figcaption></figure>

$$ \beta \frac{\partial C_m}{\partial t} + (1 - \beta) \frac{\partial C_{im}}{\partial t} = D_m \frac{\partial^2 C_m}{\partial x^2} - v_m \frac{\partial C_m}{\partial x} \quad (8a) $$

$$ \frac{\partial C_{im}}{\partial t} = \omega(C_m - C_{im}) \quad (8b) $$

where $C_m$ [ML$^{-3}$] and $C_{im}$ [ML$^{-3}$] refer to the solute concentrations in mobile and immobile regions, respectively. $D_m$ [L$^2$T$^{-1}$] refers to the (longitudinal) hydrodynamic dispersion in the mobile region while $v_m$ [LT$^{-1}$] is the mean flow velocity in the mobile region. Here, $\beta$ [-] is the partition coefficient that describes the proportion of water in the mobile region. Higher values of $\beta$ indicate that the substantial amount of water influences the solute transport ($0 \le \beta \le 1$). When $\beta$ is equal to 0, however, the MIM model turns out to be the classical ADE (see Equation 4). $\omega$ [T$^{-1}$] is the first-order mass transfer coefficient that characterizes the mass exchange rate between the mobile and immobile regions. This coefficient is often assumed to be linearly proportional to the concentration differences between these regions. To a lesser extent, increasing values for $\omega$ [T$^{-1}$] implies that the intensive (mass) transfer takes place between the mobile and immobile regions. In Equation 8b, when the multiple-rate mass transfer exchanges is considered, however, the term $\frac{\partial C_{im}}{\partial t}$ is defined based on the memory function that characterizes the (particle) travel time distributions over the immobile regions.

To simulate karst transport, the two-region non-equilibrium (2RNE) approach is frequently used for simulating transport dynamics in karst aquifers (e.g., Assunção et al., 2023; Barberá et al., 2018; Birk et al., 2005; Ender et al., 2018; Geyer et al., 2007; Göppert & Goldscheider, 2008; Goeppert et al., 2020). Here, 2RNE only considers the flow along the karst conduit in which flow velocity is assumed to be constant whilst no water exchanges are considered in between the matrix and conduit components (Field & Li, 2011; Li & Liu, 2014). For the simulation of reactive processes in karst aquifers, 2RNE was extended by the inclusion of the retardation factor by Geyer et al. (2007). Furthermore, similar approaches such as the physical non-equilibrium transport (PNE), the chemical non-equilibrium transport (CNE), and the physical and chemical non-equilibrium transport (PCNE) models by Field and Leij (2014) are also available. For a more detailed representation of karst transport dynamics, the dual-domain based models (frequently the dual-porosity models over the dual-permeability models) have also been extended by a triple porosity concept where the conduit flow is formulated based on the turbulent flow properties (Faulkner et al., 2009; Reimann, Rehrl, et al., 2011; Reimann et al., 2014).
<figure><img src="image_17.png" /><figcaption>**Figure 13.** (a) A perceptual schematic of the dual transport domain along the (reach) of the karst conduit and an overlapping matrix. (b) The transformation of the conceptual model described in Panel-a into the mathematical concept of the Transient Storage model (TSM).</figcaption></figure>

Other adaptations of the MIM approach for karst transport comprise the Dual-advection dispersion (DADE) model and the Transient storage model (TSM). Developed by Field and Leij (2012), DADE considers the solute transport characteristics between physically interacting two conduits (named for *i* and *j* in Equation 9) whilst the mass exchange is approximated by a first-order mass exchange between the main and auxiliary conduits (Bodin, 2020; Dewaide et al., 2016, 2018; Zhao et al., 2017, 2021) formulated as followed by;

$$ \theta_i \frac{\partial C_i}{\partial x} = \theta_i D_i \frac{\partial C_i}{\partial x} - \theta_i \nu_i \frac{\partial C_i}{\partial x} + \omega(C_j - C_i) \quad (9) $$

where $\theta_i$ [$L^3L^{-3}$] is the volumetric water content in *Conduit-i*, while $C_i$ and $C_j$ are the solute concentrations in *Conduit-i* and *Conduit-j*, respectively. $D_i$ [$L^2T^{-1}$] is the (longitudinal) hydrodynamic dispersion coefficient while $\nu_i$ [$LT^{-1}$] is the mean flow velocity in the *Conduit-i*. $\omega$ [$T^{-1}$] is an effective first-order rate coefficient for the mass exchange between *Conduit-i* and *Conduit-j*. Here, $\theta_i = \frac{A_i}{A}$ and $A_i = A - A_j$ while $A_i$ represents the water-filled area for the *Conduit-i* (see more explanations in Field and Leij (2012)).

TSM is built upon the dual-domain model concept in which the preferential flow paths are considered the main flow and transport regions (Figure 13). Originally developed for solving transport problems in surface streams, TSM was adapted to simulate karst transport problems (e.g., Hensley & Cohen, 2012; Morales et al., 2010). TSM was then reformulated by Dewaide et al. (2016) for the karst aquifers in a way that karst conduits are counted as the main flow and transport regions, whereas the (transient) storage zones are either represented by the pools in which eddies, vortices, and recirculation zones frequently occur or resembled by the cave sediments and carbonate matrix with small fractures (see Figure 13) (e.g., Deleu et al., 2023; Dewaide et al., 2018; Tinet et al., 2019; Zhao et al., 2017, 2019, 2021). TSM is formulated as followed,

$$ \frac{\partial C}{\partial t} = -\frac{Q}{A} \frac{\partial C}{\partial x} + \frac{1}{A} \frac{\partial}{\partial x} \left( A D_L \frac{\partial C}{\partial x} \right) + \alpha(C_s - C) \quad (10a) $$

$$ \frac{\partial C_s}{\partial x} = \alpha \frac{A}{A_s} (C - C_s) \quad (10b) $$

where $Q$ [$L^3T^{-1}$] is the volumetric flow rate and $D_L$ [$L^2T^{-1}$] is the (longitudinal) hydrodynamic dispersion coefficient. $C$ and $C_s$ [$ML^{-3}$] represent the solute concentrations in the main flow channel (i.e., karst conduit, cave stream) and in the storage zone (i.e., matrix), respectively. While $A$ [$L^2$] and $A_s$ [$L^2$] are the cross-sectional areas of the main flow channel and storage zones, the term $\frac{Q}{A}$ [$LT^{-1}$] refers to the mean flow velocity along the main flow channel. $\alpha$ [$T^{-1}$] is the first-order mass exchange term, termed also as storage zone exchange coefficient in TSM (see Figure 13).
Overall, the MIM-based modeling approaches have strong limitations in representing the physical transport re-
gion in karst systems, particularly by representing the eddy mixing and inertial flow (Anwar & Sukop, 2009). For
that reason, there are newly proposed advanced transport models to solve the anomalous transport such as late-
time, early arrivals, and/or multi-peaked BTC in karst systems while some of them were originally developed for
the heterogeneous porous media and fractured aquifers. The fundamental concept of these advanced modeling
approaches is (again) the consideration of two-region transport. However, they further require more sophisticated
parameters in describing the mass exchange between the karst conduit and matrix. For instance, the dual het-
erogeneous domain model (DHDM) developed by Yin et al. (2022) solves the large-scale anomalous transport in
karst systems that are characterized by multiple peaks and strong tailing in BTC (Ji et al., 2022; Yin et al., 2022).
DHDM considers the combination of MRMT and DADE (Leij et al., 2012) based on the memory function and
fractional derivatives (see further details in Yin et al., 2022; Ji et al., 2022). That way, two main parameter sets in
DHDM describe the individual flow velocity (advection) and (hydrodynamic) dispersion for each mobile region.
Similarly, the multi-region advection-dispersion (MRAD) and Dual Region Advection Dispersion (DRAD)
models proposed by Majdalani et al. (2018) consider the karst flow and transport process for (at least) two mobile
regions in parallel, and the mass exchange resulted from the concentration gradient (Luo et al., 2024; Majdalani
et al., 2018).

4. Lessons Learned: Challenges in the Successful Transport Modeling Practice

Considerable efforts have been devoted to representing transport processes in karst flow models over the last
couple of decades (see Section 3). However, great challenges remain to successfully model solute transport
dynamics in karst systems and to provide robust transport predictions. The challenges spread over all parts of the
karst transport modeling practice (Section 3.1, see Figure 8) and will be elaborated in the following.

4.1. System Conceptualization Considering Data Collection and System Understanding

System conceptualization for modeling karst transport processes requires a proper understanding of the system's
physical, chemical, and bio(geo)chemical characteristics. This requires extensive data sets (Huntoon, 1995;
Jeannin et al., 2013; Parise et al., 2015; Perrin et al., 2003a). However, given that the availability of such massive
data is very limited, a proper conceptualization of the spatiotemporal variabilities of karst transport processes (i.e.,
advection, dispersion, sorption, and/or chemical transformations) is limited, too (e.g., Covington et al., 2011;
Field & Pinsky, 2000; Hartmann et al., 2014; Quinlan et al., 1996).

Knowledge accumulation by the compiled data feeds the perceptual understanding of the karst system's char-
acteristics, and thus allows the modelers to map the system's conceptual definition out. For this reason, the
required level of details for the karst transport problem often shapes our data collection strategy, and hence
designs the monitoring networks (e.g., Hunkeler & Mudry, 2014; Nativ et al., 1999; Panno et al., 2019; Steva-
nović & Stevanović, 2021; Vucinic et al., 2022, among others). In practice, however, there are certain difficulties
during the karst system conceptualization including the (a) right definition of spatial scale for the transport
problem, (b) right definition of the temporal scale of karst transport dynamics with their temporal variabilities, (c)
realistic definitions of parameter values, and (d) process dependency to the complex karst flow, and (e) issue of
parameter upscaling for the reactive transport. Each of these will be pointed out in the following.

The right definition of the scale of the karst transport problem is vital for the system conceptualization. However,
due to the karstic heterogeneity at multiple scales, the investigated/modelled domain might substantially differ
from the actual transport domain in which physical and (biogeo)chemical processes predominantly occur. This is
often the case in karst systems in that both physical and chemical parameters controlling flow and transport
processes rely on the scale being investigated (Göppert & Goldscheider, 2008; Medici & West, 2021). The spatial
extent of the transport problem is generally defined by the field experts and/or modelers themselves. As an
implicit rule, the scale of interest should not be less or greater than the scale of method (i.e., artificial tracer tests,
pumping tests) applied for the characterization of that particular process (Bakalowicz, 2005; Geyer, 2008;
Ghasemizadeh et al., 2016; Maréchal et al., 2008). Nonetheless, the methods for data acquisition and/or inter-
pretation techniques are often (more) appropriate for the continuum scale applications where the effective pa-
rameters can often be used for the entire domain (Faulkner et al., 2009; Huntoon, 1995; Reimann, Geyer,
et al., 2011; Zhang et al., 2021). For this reason, more sophisticated equipment and research techniques might be
required for the descriptions of the karst transport dynamics. For instance, the characterization of nitrate or
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Schematic representations of the model conceptualization for solute transport dynamics in karst system considering the spatial extent and temporal scale of the karst transport. The schematic also demonstrates what affects the temporal variability of the transport processes under different hydrological conditions, for example, dry and wet conditions. (a) Contaminant plume migration under the high-water flow conditions, (b) Contaminant plume migration under the low-water flow conditions. Here, $t_p$ depicts the time to the peak solute concentrations with respect to the hydrological variability and sampling frequency. $t_{p1}$ (high resolution) and $t_{p2}$ (low resolution) describe the time of the peak concentrations depending on sampling frequency.</figcaption></figure>

(dissolved) organic carbon migrations in karst systems is rather challenging due to the (contaminant) input variability (e.g., Husic et al., 2019; Opsahl et al., 2017; Orban et al., 2010; Yue et al., 2019, 2020) and less-known source history (e.g., Liu et al., 2021; Lu et al., 2023; Ren et al., 2022; Sullivan et al., 2019). Therefore, it might also be difficult to successfully capture the spatiotemporal dynamics of the contaminant migration depending on the contaminant of interest with the regular measurement devices and/or monitoring systems (Finsterle et al., 2012; Zhu et al., 2007).

The right identification of the transport timescale with its temporal resolution is crucial for the proper representation of the karst transport dynamics (e.g., Bailly-Comte et al., 2011; Covington et al., 2012; Zhang et al., 2019). However, given that karst conduit and matrix might have order-of-magnitude differences in their hydraulic and physicochemical properties (i.e., porosity, hydraulic conductivity, reactivity constants), a wide range of timescales can be identified for the karst transport (e.g., Robinson & Hasenmueller, 2017; Sivelle et al., 2019). But in the context of experimental characterization, the (karst) transport timescale is related to the process timescale along the karst conduits (i.e., hours, days), for which shorter-timescale observations with higher sampling resolution (i.e., in minutes, hourly, daily) are usually enough for capturing the karst transport dynamics (Figure 14). Without a long-term monitoring and (further) data collection strategy, however, the chemical and (biogeo)chemical processes occurring at (rather) larger timescales (i.e., annual, seasonal, decadal) might be underestimated and/or overlooked in such a data collection (e.g., Mathias et al., 2005; Tran et al., 2020; Wynn et al., 2013; Zhang et al., 2020b). Therefore, proper understanding of the (actual) mass transfer between the karst conduit and matrix is strictly limited (Brettmann et al., 1993; Cherubini et al., 2013; Zhao et al., 2017).

The temporal variability of the data is another complement to the successful conceptualization because the karst transport dynamics are often controlled by the hydrological conditions being investigated (Figure 14). For
instance, BTC resembles one of the solute transport behaviors under certain flow conditions at certain time periods (Figure 14). However, actual transport can be broadly varied by the different timescales ranging from an event to seasonal and/or inter-annual timescales (e.g., Doummar et al., 2018; Luo et al., 2024; Remondi et al., 2018, 2019; Sivelle & Labat, 2019; Zhang et al., 2020a). As a matter of fact, the same shape for BTC might not actually be reproduced again if those particular conditions (both hydrological and physicochemical) at that particular timescale are not met (e.g., Aydin et al., 2014; Chu et al., 2021; Kogovsek & Petric, 2014; Morales et al., 2007, 2018; Ravbar et al., 2012, among others). Therefore, the system conceptualization is restricted to certain timescales and might not properly represent the overall karst transport dynamics especially when the hydrological or (biogeo)chemical conditions change or vary temporarily.

Another limitation in the system conceptualization is the identification of the realistic values for the system's parameters since both hydraulic and transport parameters can vary by orders of magnitude at differing spatial scales—from a few centimeters to hundreds of meters or more (see also Figure 2)—(Berkowitz, 2002; Giese et al., 2018; Medici et al., 2019; Medici & West, 2021). Hydraulic conductivity and dispersivity are such parameters very sensitive to scale variations, thus causing substantial variations in the spreading of the contamination plume in the karst aquifer. Likewise, chemical parameters (i.e., sorption parameters) obtained under the laboratory conditions might largely differ from those measured at the field scale (Faulkner et al., 2009; Witthüser et al., 2003). For instance, the partition coefficient ($K_d$) and mass transfer coefficient (i.e., first-order mass transfer coefficient) could not be generalized (Stockmann et al., 2017; Swami et al., 2018) because their values are extremely varied (often by order of magnitudes) at the field-scale measurements depending on the matrix hydraulic properties (e.g., surface area, matrix porosity, permeability) and the fracture/conduit characteristics (Howroyd & Novakowski, 2021; Lehmann et al., 2022; Maloszewski & Zuber, 1993; Zhou et al., 2007). Similarly, the values for reactive transport parameters (e.g., dissolution reaction rate) could also be misleading particularly when they are scaled up or extrapolated, because of the oversimplification of the multifaceted processes that are varied both by chemical heterogeneities in the matrix and local/global groundwater chemistry (Howroyd & Novakowski, 2021; Jiang et al., 2023; Katz et al., 1998; Kaufmann et al., 2014; Mohammadi et al., 2021; Maqueda et al., 2023). Therefore, a realistic definition of the reactive transport parameters is strictly limited by the well-defined karstic sites or site-specific applications (i.e., single-well push-pull test) (Hillebrand et al., 2012a, 2012b, 2015; Priebe et al., 2022; Schiperski et al., 2016; Tran et al., 2020, 2021).

Field-scale experiments, particularly artificial tracer tests, are frequently performed under highly vulnerable aquifer conditions in which sinkholes and/or swallets are frequently utilized as tracer injection points (Jones, 2019; Li et al., 2008; Morales et al., 2007; Sinreich, 2011). In such experimental designs, as the preferential flow paths predominantly shape the spatiotemporal behavior of the contamination plume, karst transport dynamics are often characterized by the advection-dominated transport between the injection and sampling points (i.e., between sinkhole and observation well, Figure 14). As a result, the role of matrix diffusion and mass exchange processes between the karst conduit and the surrounding matrix might be underestimated. Furthermore, due to the randomly varying flow velocities in karst conduits, physical parameters such as hydraulic conductivity, velocity, dispersivity, solute travel time (Kavousi et al., 2020; Rehrl & Birk, 2010; Zhao et al., 2017) and chemical parameters including sorption and mass exchange parameters are rarely understood (Barberá et al., 2018; Field & Pinsky, 2000; Haggerty et al., 2004). Consequently, the observed data might not provide a complete picture of the karst system's (overall) hydrochemical functioning and transport characteristics due to the limited transport distance and (certain) flow conditions under specific experimental conditions.

Geochemical reactions (i.e., complex physicochemical alterations) and dissolution/precipitation processes occur at the mineral surface at the pore-scale level (Seignur et al., 2019) (see also Figure 5). As a result, the definition of reactive transport processes for karst aquifers (from a couple of meters to a few kilometers) requires upscaled parameters (i.e., permeability, effective diffusion coefficient, reaction rate) where the interferences of karst flow, transport, and chemical parameters (i.e., evolving reactive surface area, reactivity rate, tortuosity, porosity) can be developed and assumed to be invariant. For instance, incipient karst generation is often coupled by flow, reactive transport of Calcium ($Ca^+$), and the evolution of fracture and/or conduit (Aliouache & Jourde, 2024; Dreybrodt, 1996; Jiang et al., 2022). Nevertheless, introducing chemical reactions that are fundamentally formulated at the pore scale complicates the upscaling of effective parameters because the reaction length scale is often smaller than advective or dispersive length scales (Ladd & Szymczak, 2021; Yu et al., 2022). Given the complexity of the karst aquifers, a wide range of chemical reactions and/or perturbations can occur from microscale (i.e., pore scale) to macroscale (i.e., aquifer scale). But the karst flow models coupled with the carbonate dissolution kinetics
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Model selection based on the dominant processes across the different karstic environments (Here, the karst systems are classified based on the study by Goldscheider et al. (2020)). (a) Generic perceptual descriptions of the broad range of karstic environments under differing climatological, hydrological, and hydrogeochemical conditions, together with the (possible) karst transport dynamics, (b) The selected (conceptual) model structures with respect to the right choice of the dominant process associated with corresponding environment. Here, only lumped simulation models with their (possible) structures are depicted to emphasize that a wide range of the (parsimonious) model structure is possible depending on the definition of the dominant processes and the problem of interest. Here, the acronyms indicate the different components of karst systems such as *E* (Epikarst), *C* (Conduit), and *M* (Matrix), as well as the model state variables including *P* (Precipitation or Rainfall), *ET* (Evaporation/Evapotranspiration), *S*<sub>b</sub> (Sublimation), and *SWE* (Snow Water Equivalent).</figcaption></figure>

remain a challenge (e.g., Guo et al., 2021). This, therefore, may not allow for a proper definition of the effective model parameters for reactive transport at the karst aquifer scale.

All these difficulties explain why a robust conceptualization for the karst transport dynamics is inherently limited, and why there are only a few real-world applications of the karst transport models that solve the classical ADE (see Section 3.2). In this matter, the limitations in system conceptualization for karst transport point out the need for in-depth discussions on how to better align the investigated (transport) domain and the model domain in karst systems, as well as for finding better ways to collect physical and (biogeo)chemical data across the scales.

## 4.2. Selection and Parameterization of the Karst Transport Model

### 4.2.1. Model Selection Based on the (Right) Choice of Dominant Process

Model selection with a detailed understanding of its assumptions and limitations is one of the key steps for successful transport modeling in karst systems. However, such a choice from a model plethora is not an easy task (Horton et al., 2022; Masciopinto et al., 2021) as each model has its own features (i.e., numerical formulations, parameters) that may be appropriate for a specific application and/or specified conditions (Figure 15). Certainly, the modeler's perceptual understanding of the karst flow and transport processes can also be judicious in selecting a candidate model structure (e.g., Cousquer & Jourde, 2022; Hartmann, Wagener, et al., 2013; Hartmann et al., 2014). Of the numerous criteria, however, the objective of the study, data availability, and model structural adequacy generally define our most suitable choice in simulating karst transport dynamics (see Section 3).
In general, an (accurate) representation of the flow and transport processes in the model structure can be seen as a crux of the robust model selection (Figure 15). But, as we have often a poor understanding about the complicated physical, chemical, and (biogeo)chemical processes, special care is required for the representation of the karst transport dynamics in the (selected) model structure with adequate numerical formulations and parameters. From a modeling point of view, however, the representations of multiple processes at differing spatiotemporal scales under a broad range of karstic environments—in which different transport problems might be encountered—with similar model structures is neither possible nor enough to realistically simulate the karst transport dynamics (e.g., Chang et al., 2017, 2021; Cen et al., 2021; Hublart et al., 2015; Makropoulos et al., 2008; Massmann et al., 2005, among others). For this reason, the overall success in the robust model selection can be attributed to the “right” choice of the key process that predominantly controls the karst system’s hydrological and/or physicochemical behaviors by the available data (Hartmann, 2018; Hartmann, Barberá, et al., 2013).

In karst systems, since the (accurate) process description is restricted by data availability, the principle of parsimony—the minimum level of the model complexity—(Schoups et al., 2008; Wagener et al., 2001, 2003) has been recommended to ensure the adequate model complexity in terms of process representation with the data we have (Hartmann, 2018; Hartmann et al., 2014; Kavousi et al., 2023). In general, based on the data/knowledge on the spatial variabilities of the dominant processes, the model structure can be implemented in a lumped (parsimonious), distributed (more process-based), and semi-distributed way. Amongst them, the spatially lumped karst simulation models (i.e., CIMMs, LTMs) are often preferable to simulate the karst transport dynamics in data-scarce karst regions (Barrett, 1996; Barrett & Charbeneau, 1997; Hartmann et al., 2014, 2016; Scanlon et al., 2003), while spatially distributed process-based models (i.e., ATMs) require a large type of observations to resemble varying physical and chemical processes with a high dimensional of parameter sets (Assari & Mohammadi, 2017; Chang et al., 2019; Kavousi et al., 2020; Xu et al., 2015). For instance, the mass transfer between the karst conduit and matrix needs to be considered with additional parameter sets that resemble the mass exchange in such models that is, 2RNE, TSM, DHDM (Assunção et al., 2023; Ender et al., 2018; Schiperski et al., 2022; Yin et al., 2022) (see Subsection 3.2.2). As a consequence, the successful applications of the advanced transport models are strictly limited by the well-defined sites, small-scale field, and/or lab-scale investigations where the rather site-specific parameters (i.e., mass transfer rate, sorption constants) are obtained (Lehmann et al., 2022; Schiperski et al., 2022; Shu et al., 2020; Wang et al., 2022). Because of this limitation, the semi-distributed karst simulations models, like Varkarst (Hartmann, Barberá, et al., 2013) and TCD-Dublin (Gill et al., 2013), are frequently served as the optimal model structures to ensure the balance between the process representation and data availability (Borghi et al., 2016; Hartmann et al., 2014). The complexity of such models is often enough to represent karst transport dynamics as compared to the parsimonious lumped models, as well as better choice with the lower/adapted data need than that of distributed counterparts (Hartmann, 2016; Teixeira Parente et al., 2019; Zoghbi & Basha, 2020). But despite the semi-distributed, physically oriented flow path simulation of these models, a big gap remains in their presentation of transport dynamics which has been still relying on simple CIMM approaches.

Overall, finding the most adequate transport model structure and degree of complexity is a particular challenge for karst transport modeling. But the optimal structure should be as simple as possible and as complex enough to better represent the karst system's underlying key physical and chemical phenomena in accordance with the available data. In principle, the complexity of the selected model structure can only be increased when the data amount (i.e., multiple observations) and its quality ensure the proper representations of the multiple processes in the model structure (Hartmann, 2018; Hartmann, Barberá, et al., 2013; Hartmann et al., 2014; Milanović & Vasić, 2015), thereby increasing the model realism in terms of process representation. Nonetheless, seeking a trade-off between the (robust) representation of karst transport processes and their proper transcriptions in the model structures based on the principle of parsimony has not yet been fully explored, yet. Therefore, which type of model we could use, and how complex it would be are those questions still open to discussion.

### 4.2.2. Model Parametrization: Parameter Non-Uniqueness and Transferability

A successful transport model parametrization must ensure robust parameter estimation with reliable model predictions (Chang et al., 2017; Çallı, Bittner, Liu, et al., 2023; Hartmann et al., 2017; Sivelle et al., 2022). Yet, a common tenet in karst modeling studies is to obtain physically reasonable parameters in comparison with the actual site characteristics. For that, a successful model parametrization is required. However, due to the limited system knowledge with the insufficient field estimates, a robust parameterization in karst transport models is a
a) An ensemble of behavioural
model predictions

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

Equifinal solutions

Uncertainty bound

b) Equifinality concept

<figure><img src="image_22.png" /></figure>

v₁

κ₁

d₁

<figure><img src="image_23.png" /></figure>

0.5v₁

κ₁

2d₁

<figure><img src="image_24.png" /></figure>

0.5v₁

2κ₁

d₁

**Figure 16.** Description of the model calibration results for the karst transport model based on the concept of model equifinality, both in the model predictions and in model parameters. (a) Calibration of a karst transport model ended up with the equifinal simulations within a certain range of uncertainty (depicted as an ensemble of the behavioral simulations with the shaded area), (b) The (hypothetical) description of the equifinality concept in the karst transport model parameters resulting from the equifinal solutions as depicted in Panel (a). Here, multiple feasible parameter-sets provide equally-well hydrological responses (i.e., spring discharge, *Q*) and physicochemical (solute concentration, *C*<sub>out</sub>) with respect to the (selected) calibration metric. In the equifinal model results, finding the “best” or “most acceptable” solution is not necessarily possible, and hence the results might be misleading. In the figure, the acronyms are *v*₁ (flow velocity), *K*₁ (hydraulic conductivity), and *d*₁ (cross-sectional area).

challenging task, particularly in two ways. First, how to better translate the karst transport dynamics into the model structure with sufficient parameters is a serious matter because additional parameters are involved in the karst transport model (i.e., dispersivity, mass transfer coefficient) compared to the flow models. Second, since most parameter values cannot be assessed from direct measurements, inverse modeling (i.e., model calibration) is necessary to describe the key physical and (biogeo)chemical processes (Figure 16). In general, the former task could be handled by sticking to the principle of parsimony where the minimum number of parameters meet with observational data sets (Hartmann, 2018; Kavousi et al., 2020, 2023), whilst the latter is often rather difficult because it will possibly go along with overparameterization. This thus limits the physical interpretation of the parameter values, as well as reduces the model reliability with a certain degree of uncertainty (Figure 16) (Çallı, Bittner, Liu, et al., 2023; Hartmann, 2016; Teixeira Parente et al., 2019; Zhang et al., 2020b; Zoghbi & Basha, 2020).

The direct field measurements and/or lab-scale experiments are favorable for the improved parameterization of the karst transport models in that the spatial and temporal variabilities of the transport processes could be defined by the well-defined parameter values. In principle, since the (more) physically based models (i.e., ATMs) consider the parameter heterogeneities along the domain of interest, they are often built upon the experimental data that describes the physical/chemical significance of the transport processes with the effective values (Dentz et al., 2011; Li et al., 2008; Orban et al., 2010; Pérotin et al., 2021; Zhang et al., 2021). Nevertheless, it might not be still not possible to define the specified values for some of the parameters, especially for the larger scale transport modeling (e.g., aquifer scale or regional scale) such as conduit/fracture geometry (i.e., orientation,
spacing, aperture) and conduit/fracture networks (i.e., distributions, intersections, tortuosity factor). Consequently, such parameters are frequently unknown, and thus required to be calibrated (e.g., Barberá et al., 2018; Dreybrodt, 2012; Faulkner et al., 2009; Florea & Wicks, 2001; Geyer et al., 2013; Robineau et al., 2018, among others). For instance, to parametrize the Transient Storage Model (TSM) where transport distance is divided into several smaller reaches needs several effective parameter sets to represent the transport processes along the karst conduit. Here, the length of each reach can be defined by the (extensive) field works, whereas dispersion and cross-sectional areas for the main flow channel and transient storage are still obtained by the model calibration (Cholet et al., 2017; Deleu et al., 2023; Dewaide et al., 2016; Tinet et al., 2019; Zhao et al., 2021). Likewise, the mass transfer coefficient is often considered as a calibration parameter, especially for the dual-domain transport models (i.e., MIM, 2RNE, TSM). Considering all, model calibration is an inevitable part of successful model applications, both for the lumped and advanced karst transport models.

As an inverse solution of the karst system's transport dynamics, a successful calibration seeks to infer a (unique) parameter set that resembles the process of interest at the spatiotemporal domain of interest. Yet, as there is no unique solution for an inverse problem (e.g., Hartmann, 2018; Sivelle et al., 2020; Valota et al., 2002), it is often difficult to link the karst transport parameters with the actual transport dynamics. For that reason, karst transport models are more prone to suffering from the issue of non-uniqueness as compared to the flow models, particularly in describing and (physically) interpreting transport parameters (i.e., flow velocity, dispersivity) (e.g., Adinehvand et al., 2017; Chang et al., 2017; Field & Li, 2011; Johnston et al., 2009; Massmann et al., 2005, among others). This often results in overparameterization and equifinality as exemplified in Figure 16. While the overparameterization is often a serious problem for the advanced/distributed models where the multiple parameters are represented in the model structure despite the less information available (Fischer et al., 2024; Kavousi et al., 2020, 2023), it might be rather less of a problem for the lumped transport models (i.e., CIMMs, LPMs) where the direct observations are not possible due to the aggregation of the karst flow and transport processes along the domain of interest (Hartmann et al., 2014; Sivelle et al., 2022).

Model overparameterization or equifinality manifests itself in larger uncertainties in the model predictions and parameters (Figure 16). For this reason, finding the “best” or “most acceptable” solution is not necessarily possible. Therefore, model predictions (or equifinal solutions) can be strongly misleading (Beven, 2006; Husic et al., 2020; Wang et al., 2020; Zhang et al., 2019, 2020b). In theory, parameter identifiability tests can quantify the equifinality and help to reduce it with additional information or model structure simplifications. However, these improvements may be limited because the physical and (biogeo)chemical parameters that control the migration of the contamination plume are highly sensitive to the karst system's physical and (biogeo)chemical heterogeneities (Bonanno et al., 2022; Jia et al., 2018; Zhang et al., 2020b; Zoghbi & Basha, 2020). For instance, the interpretation of mass exchange rate, spatially varying sorption properties, and reactive kinetic surface areas for the dual-domain transport models (e.g., 2RNE, MIM, TSM) are often problematic and less reliable for predicting the karst transport dynamics due to the lack of exhaustive information. Therefore, the non-identifiability problem in such transport parameters might be inherent.

A (successful) application of the karst transport model is often unique in solving (particular) flow and transport conditions in the domain of interest. For this reason, the transferability of the (model-specific) calibrated parameters for the (rather) larger and/or smaller transport domains, for various karstic environments, and for different models under different concepts is strictly limited, too. As a result, using the literature values for some of the transport parameters might not be entirely possible simply because karst transport parameters and their values are mostly unique in place and time.

Overall, it is particularly important to emphasize that the substantial limitations behind the successful model parametrization of karst transport models reside especially in the governing solutions (i.e., ADE) that are particularly developed for the porous aquifer-specific transport problems (Bear et al., 2012; Faulkner et al., 2009; Field & Pinsky, 2000; Zoghbi & Basha, 2020). For that reason, the targeted solutions cannot be easily transferred for the highly heterogeneous karst aquifers in which the additional mathematical complexities (e.g., the inclusion of mass transfer coefficients, sorption parameters in classical ADE) are not easily incorporated and over-parameterization is likely to arise.
### 4.3. Uncertainty of the Model Simulations

Characterization and reduction of the uncertainty is one of the priorities for the evaluation of transport model applications as it ensures the robustness and reliability of the model predictions (Chang et al., 2021; Çallı, Bittner, Liu, et al., 2023; Hartmann, 2018; Hartmann, Wagener, et al., 2013). For this reason, treating uncertainty with its (potential) sources requires special care throughout the karst transport modeling practice, especially for identifying the ways for the model improvements.

Taxonomically, the larger source of uncertainty in the model predictions is attributed to the epistemistic uncertainty due to the lack of system knowledge and/or ignorance of some key processes in the formal model structure, rather than the stochastic uncertainty arising from the natural variability of the hydro(geo)logical systems (Beven, 2023; Beven et al., 2011; Reichert et al., 2021). This uncertainty type is particularly relevant for the karst flow and transport models because they are mostly built upon insufficient knowledge and/or site-specific information, and thus often shaped by the expert elicitations on the complex physical and chemical properties of the karst systems (e.g., Hartmann et al., 2014; Padilla & Vesper, 2018; Zhang et al., 2017). In general, uncertainty in karst transport modeling is raised from the different stages of the karst transport modeling practice (see also Figure 8) while the main sources of uncertainty in the model predictions can be attributed to the following sources: (a) data uncertainty, (b) model structural uncertainty (both conceptual and numerical), (c) parameter uncertainty, and (d) uncertainty associated with the unknown future behavior of the karst system.

Despite clear evidence that the uncertainties in groups (a)-(c) are substantial (see previous subsections), only a few karst transport studies with LTMs have quantified the uncertainty of their transport simulations (Çallı, Bittner, Liu, et al., 2023; Chang et al., 2021; Mudarra et al., 2019). Amongst them, Chang et al. (2021) were the only ones trying to reveal structural uncertainty by comparing different model structures. For ATMs, different modeling approaches (i.e., ADE, MIM, PCNE) were also tested for the simulations of the shape characteristics of BTCs (e.g., multi-modality and heavy tailing) obtained from different karst conduits by Field and Leij (2014). This model comparison research provided a better understanding of the suitable model structure and parameterization of the karst transport process (Massmann et al., 2005; Zhao et al., 2019). Therefore, the disagreements in the model predictions with differing complexity for the same applications might provide a complete understanding about the suitable model structure and improved parameterization (Assunção et al., 2023; Dong et al., 2023; Massmann et al., 2005; Zhao et al., 2019), as well as indicate the scale dependency of the model parameters and process interactions (Goeppert et al., 2020; Kaufmann, 2003; Liu et al., 2021).

Concerning the uncertainty associated with the unknown future behavior of the karst system, the realistic representations of the karst system's current and historical transport conditions may not guarantee that future conditions (e.g., unknown future state variables) will successfully be met by the chosen model structure, as well as by the obtained (behavioral) parameter set (e.g., Beven & Chappell, 2021; Oreskes et al., 1994; Wagener, 2022). Regarding this, one of the critical tasks for the karst modelers is to assess the model results by seeking (more) reasonable answers for to what extent the model simulates the process behavior under (globally) changing hydrological and physicochemical conditions. Improving our physical understanding of karst systems and implementing it into new and better model structures is an ongoing challenge for the successful modeling practice.

### 5. Ways Ahead of Modeling Karst Transport Dynamics and Future Perspectives

Studies on modeling of karst flow and transport processes have considerably increased over the last decades. Several factors can play a role in such a progress including the variety of collected data (i.e., multiple observations), its accuracy, and increments in computational power and capacity (both accuracy and efficiency). But the driving force is the need for solving transport problems in karst water resources across the world. Nonetheless, as a community, we are still at the stage of model improvements for the successful real-world applications of the karst transport models by confronting the main challenges (i.e., system conceptualization, model parametrization, uncertainty quantification) to ensure model robustness and reliability (see Section 4). Along with these (long-lasting) challenges, however, new challenges and ways are ahead of us, especially considering the deterioration of karst water quality in a changing world (e.g., climate change impacts on karst water quality and ecosystem services). For that, in this review, we highlight four key directions that could be addressed for better representing and predicting the solute transport behavior in karst systems, as well as for exploring karst water quality risks and potential outcomes.
## 5.1. Collect and Provide More Data for Improved Karst Transport Understanding

To reasonably predict karst transport dynamics, the solute transport models often require spatially and temporally continuous and consistent data sets. However, the collection of chemical and (biogeo)chemical data (i.e., nitrate, arsenic, biological and microbial contaminants, microplastics) is still sparse, and mostly inconsistent with the temporal coverage of the physical data (i.e., spring discharge, groundwater head, water pressure). As a result, our proper understanding of the karst system-specific contaminant transport characteristics, as well as the knowledge of reactive transport processes and (biogeo)chemical transformations are substantially limited. Nonetheless, taking long-term precautions to protect karst water quality from anthropogenic impacts and global changes (i.e., climate, land cover, and land use changes, soil degradation) often rely on a better understanding of chemical and biogeochemical processes that influence the long-term behavior of the solute movement. For this reason, a particular attention is required for an improved understanding of the (multiple) physicochemical interactions between the preferential flow paths (e.g., conduits, fractures) and overlapping matrix, thereby developing more integrated karst modeling approaches that simultaneously consider flow, transport, and (bio)geochemical reactions in karst systems.

Karst scientists and water authorities are collecting more data than ever before. However, the collected data are rarely published and shared with the community, especially those used for flow and transport modeling studies (i.e., karst spring $δ^{18}$O isotope ratios). Some working groups under specified goals, for instance, SISAL (Speleothem Isotope Synthesis and Analysis) provide a readily available global database for the speleothem isotopes collected from individual speleothems and grouped them by the cave locations across the world (Atsawawaranunt et al., 2018; Comas-Bru et al., 2020; Kaushal et al., 2023). There is, therefore, a need for a common database in which we can gather the original data on karst transport dynamics from different experimental and modeling works (i.e., tracer tests, hydrochemical studies, hydraulic tests, cave isotopic and hydrochemical experiments) similar to the few and valuable attempts to collect and provide karst data by Jourde et al. (2018), Petrič et al. (2020), Olarinoye et al. (2020).

## 5.2. Fill the Gap Between Oversimplified LTM and Data-Demanding ATMs

The critical discussion on how to find a tradeoff between the model complexity and the choice of the right process controlling the system's physicochemical behavior (with a reasonable number of parameters) with respect to the data availability is a requirement for the completion of the successful transport modeling practice in karst systems. Due to suitability in data-rare sites, mathematically simpler LTM have outpaced their sophisticated ATM fellows that are represented by more variables with multiple parameters. But although the LTM are often served as valuable heuristic tools for modeling karst transport processes, they are still the most simplistic approach in representing rather complicated karst transport dynamics. By comparison, due to the (substantial) data need, ATMs are often inapplicable for resembling the karst transport dynamics, especially for the larger scale field applications. Therefore, we see a particular need for representing karst transport processes with more process-oriented and well-constrained models that can be used for alternative modeling approaches to better resemble the spatiotemporal dynamics of contaminant migration in karst systems.

In this direction, this review reveals that there is an emerging modeling concept that is coming outside of the karst. StorAge Selection (SAS) Function Approach (hereinafter referred to as the SAS model) can be an alternative modeling approach to simulate the karst transport dynamics (Figure 17). The SAS theory was built by Botter et al. (2011), Harman (2015), and Van Der Velde et al. (2012) for the catchment-scale flow and transport modeling applications, while the successful applications were also extended to simulate the karst transport processes by Zhang et al. (2021) and Çallı, Bittner, Liu, et al. (2023). The SAS model defines the relationship between the distribution of water ages in the (different) reservoir storage (i.e., soil, groundwater) and the ages exiting as the system's outflux, either discharge or evapotranspiration, or both (Figure 17). Compared to LTM, the SAS model considers the time-variant characteristics of TTDs by revealing the time-dependent mixing, storage, and discharge characteristics (e.g., Benettin & Bertuzzo, 2018; Botter, 2012; Botter et al., 2011; Harman, 2015; Rinaldo et al., 2015; Schwemmle & Weiler, 2024, among others). Therefore, the SAS model is structured by rather fewer parameters as compared to ATMs, and thus allows better identification of the model parameters while reducing overparameterization and equifinality issues.
<figure><img src="image_25.png" /><figcaption>**Figure 17.** Schematic representations of (a) StorAge Selection functions, (b) Precipitation and discharge age distributions in the SAS model conceptualization, (c) Selected functional form for the description of water age distributions in the SAS model (here, younger water age preference is described, often known as the power-law distribution), (d) Hypothetical demonstration of karst groundwater age distributions at a given time, *t*<sub>c</sub>, which only resembles a possible age distribution (e.g., TTDs) among multiple other possibilities because the SAS model only derives the system's age structure at which we can measure the solute concentration (for our case description only at karst spring's outlet). Here, the figures (a), (b), and (c) were adapted for the karst systems based on the SAS model descriptions in Benettin et al. (2022), Harman (2015), and Harman (2019). In the depicted figures, *S*<sub>T</sub> refers to the age-ranked storage that describes the amount of water in storage whose age is younger than (or equal to) a certain age *T* at the current time *t*<sub>c</sub>. *ω*(*S*<sub>T</sub>, *t*) is the probability density function of the travel time expressing each water age (*T*) removed from the storage (for further details see Benettin and Bertuzzo (2018) and Benettin et al. (2022)).</figcaption></figure>
### 5.3. Develop Adequate Approaches for Reactive and Sediment Transport in Karst Systems

Karst transport model is frequently applied for the predictions of solute transport dynamics in karst systems. However, some agents such as vapor, heat, and (suspended) sediments (i.e., particles and colloids) can also play a role in defining solute transport behaviors. Suspended sediments are such agents that are often induced by the turbulent flow along the karst conduits and cave streams (Jukić et al., 2022; Mueller et al., 2023; Richter et al., 2022). In this respect, understanding the (suspended) sediment transport is particularly important for karst transport processes in two ways. First, since (suspended) sediments can adsorb a wide range of contaminant types (i.e., organic substance, hydrocarbons, heavy, and trace metals, (micro)biological contaminants), they often act as storage and/or transformation regions for the contaminants (Herman et al., 2012; Mahler et al., 1999, 2000). Second, they can attenuate the contaminant of interest as a zone of remobilization (Bettel et al., 2022; Nerantzaki et al., 2015). Furthermore, the origin of the sediments provides evidence of transport conditions during the contaminant events taking place (Bosch & White, 2018; Çelik et al., 2022). Besides, the microbial transport in karst aquifers is frequently approximated by the suspended sediment transport (e.g., Goldscheider et al., 2010; Göppert & Goldscheider, 2008). However, the transportation and transformation of the contaminants with the suspended materials (e.g., particle-facilitated transport) are not directly measurable and observable. For these reasons, their characterizations require high-resolution monitoring and suitable technologies. Therefore, future research is recommended to address contaminant fate and transport processes due to the (suspended) sediment and/or particles at high temporal resolution.

Modeling of reactive transport processes in karst aquifers is demanding mainly due to the multiple controlling factors that result in complex geochemical and/or biogeochemical reactions (Xiong et al., 2021; Xu et al., 2018). Such complexities are further magnified by which multiple reaction constituents and/or multi-phases (i.e., DNAPL, LNAPL) are involved in the karst flow and transport at the domain of interest. As a result, modeling of reactive transport in karst aquifers is still very limited. Despite the valuable attempts at the aquifer scale application of reactive transport models (e.g., Priebe et al., 2022), the karst applications are frequently for the specific modeling purpose (e.g., either site-specific or system-specific, or both) and they generally neglect certain (biogeo)chemical processes to simplify the complex chemical processes (e.g., Tran et al., 2020, 2021). Nonetheless, understanding reactive transport is critical as they often shape the (long-term) spatiotemporal behavior of the contaminant plume. Herein, we believe that reactive transport models including biogeochemical transport models might provoke and/or necessitate cross-disciplinary research to realistically model reactive processes in karst aquifers. This could, therefore, motivate the karst researchers and modelers to explore the boundaries and threshold(s) of the successful modeling practice based on the karst system's complexity and process representation in the formal model structure. Even though this attempt is accompanied by various scientific and operational difficulties to be handled (e.g., upscaling issues both for parameters and processes).

### 5.4. Estimate the Impact of Climate and Land Use Changes on Karst Water Quality

Karst water resources are highly regulated, stressed, and managed by the human dimension, especially by the operational water management (i.e., dam and reservoir capacity) and increased water consumption (i.e., pumping rate, water needs). For this reason, greater uncertainty is available for the karst system's (natural) hydrological and hydrochemical functioning. In this matter, the impacts of human-induced factors on karst flow and transport dynamics (i.e., pumping rate, sectoral water needs) need to be considered. But up until now, the human dimension on karst transport dynamics is overlooked and/or underestimated in the real-world applications of the karst transport models while the focus is often more on naturalized aquifer conditions. Therefore, coupling the knowledge of human-induced factors with the (naturalized) karst system dynamics at differing scales under different scenarios (i.e., the impact of groundwater abstraction rates on karst water quality) is one of the key research areas to provide more realistic predictions, as well as to support decision-makers for sustainable water management in karst systems.

Contaminant transport problems in karst water resources are frequently denoted as (saturated) groundwater-zone specific issues, and the role of epikarst in karst transport dynamics is significantly simplified or completely ignored from the (formal) model structures. Nonetheless, epikarst is a critical zone whereby (contaminant) mixing, storage, and transport characteristics are largely controlled (e.g., Dal Soglio et al., 2020; Flynn & Sinreich, 2010; Li et al., 2017; Perrin et al., 2003b; Zhang et al., 2017, among others). Furthermore, the epikarst forms the reactive and biogeochemical processes (i.e., sorption, microbial activities, and oxidation), and significantly
modifies the (initial) composition of the contaminants released from the ground surface before reaching the saturated groundwater table (Druhan et al., 2021; Dubois et al., 2014; Oster et al., 2021; Zeqiraj, 2022), thereby predominantly controlling the karst system's (biogeo)chemical functioning (Anna-Neva et al., 2021; Laincz & Hays, 2014; Panno et al., 2001). Considering all these, a better understanding of epikarst transport dynamics is a fundamental need to protect and manage karst water quality and dependent ecosystem functioning. Nevertheless, due to the unbearable complexity and substantial heterogeneity of epikarst, scientific attempts to develop epikarst-dedicated transport models are quite rare.

Deterioration of karst water quality is not only a local scale issue in which we can only apply a site-specific solute transport model to predict the spatiotemporal behavior of contaminant migration, but it is also a global phenomenon that we shall further consider global changes and their potential impacts on the karst water quality and dependent ecosystems. However, to the best of our knowledge, the effects of the global changes on the karst water quality in the context of the transport modeling are not adequately addressed at local/global scales. For that reason, the impact of local/global changes on the karst water quality brings us a fundamental responsibility to make robust predictions in a changing world, despite the substantial uncertainty. Therefore, future works call attention to what extent (and how) karst transport dynamics can (or could be) linked to the altered boundary conditions due to human interventions and global changes, in particular, due to climate change, land-use and/or land-type changes, soil degradation, water management alterations, and groundwater depletion.

## 6. Conclusions

The threat of contamination in karst water resources has come to light over the last couple of decades. Nonetheless, the multiscale physical and (biogeo)chemical heterogeneities still pose an enormous challenge for properly characterizing, understanding, and modeling flow and transport dynamics in karst water resources. In this manner, we still have substantial limitations in providing more reasonable and reliable transport model predictions, thus handling the modeling issues such as parameter scaling, model non-uniqueness, and uncertainty. But beyond that, both anthropogenic impacts and global changes (i.e., climate change, land use and land cover changes, soil degradation, urbanization, agricultural practices, and water management alterations) have been seriously deteriorating the quality of karst water resources and dependent ecosystem functioning. Since the karst system's functioning both hydrodynamically and biogeochemically is also controlled by such global changes in a way that they can shape the karst flow and transport characteristics, understanding karst transport dynamics is critical to successfully predict the impacts of such factors on the karst water quality and dependent ecosystems. With this review, we hope to provide the necessary overview of available methods, challenges, and directions to support and improve the present and future management of water quality in different karst regions across the globe.

## Data Availability Statement

No original/new data were collected and/or developed for this review. The paper is a synthesis of the published works. All the figures were created in CorelDRAW Essentials 2021.

## Acknowledgments

Kübra Özdemir Çallı, Daniel Bittner, and Andreas Hartmann (TUD) and Gabriel Chiogna and Beatrice Richieri (TUM) are supported by the German Research Foundation (DFG), project titled "Robust Conceptualisation of KArst Transport (ROCKAT)," Grants HA 8113/6-1 and CH 981/6-1, respectively. Süleyman Selim Çallı was supported by the Scientific and Technological Research Council of Türkiye (TÜBİTAK, Grant 1059B142000592). Open Access funding enabled and organized by Projekt DEAL.

## References

Adinehvand, R., Raeisi, E., & Hartmann, A. (2017). A step-wise semi-distributed simulation approach to characterize a karst aquifer and to support dam construction in a data-scarce environment. *Journal of Hydrology*, 554, 470-481. https://doi.org/10.1016/j.jhydrol.2017.08.056

Aliouache, M., & Jourde, H. (2024). Influence of structural properties and connectivity of initial fracture network on incipient karst genesis. *Journal of Hydrology*, 640, 131684. https://doi.org/10.1016/j.jhydrol.2024.131684

Aliouache, M., Wang, X., Jourde, H., Huang, Z., & Yao, J. (2019). Incipient karst formation in carbonate rocks: Influence of fracture network topology. *Journal of Hydrology*, 575, 824-837. https://doi.org/10.1016/j.jhydrol.2019.05.082

Anderson, M. P., Woessner, W. W., & Hunt, R. J. (2015). *Applied groundwater modeling: Simulation of flow and advective transport*. Academic press.

Anna-Neva, V., Lehmann, M. F., Hermann, R., Grathwohl, P., Nia, B., Kappler, A., et al. (2021). Fate of nitrate during groundwater recharge in a fractured karst aquifer in Southwest Germany. *Hydrogeology Journal*, 29(3), 1153-1171. https://doi.org/10.1007/s10040-021-02314-2

Anwar, S., & Sukop, M. C. (2009). Lattice Boltzmann models for flow and transport in saturated karst. *Groundwater*, 47(3), 401-413. https://doi.org/10.1111/j.1745-6584.2008.00514.x

Aquilanti, L., Clementi, F., Landolfo, S., Nanni, T., Palpacelli, S., & Tazioli, A. (2013). A DNA tracer used in column tests for hydrogeology applications. *Environmental Earth Sciences*, 70(7), 3143-3154. https://doi.org/10.1007/s12665-013-2379-y

Aquilina, L., Ladouche, B., Dörfliger, N., & Bakalowicz, M. (2003). Deep water circulation, residence time, and chemistry in a karst complex. *Groundwater*, 41(6), 790-805. https://doi.org/10.1111/j.1745-6584.2003.tb02420.x

Assari, A., & Mohammadi, Z. (2017). Assessing flow paths in a karst aquifer based on multiple dye tracing tests using stochastic simulation and the MODFLOW-CFP code. *Hydrogeology Journal*, 25(6), 1679-1702. https://doi.org/10.1007/s10040-017-1595-z
Assunção, P., Galvão, P., Lucon, T., Doi, B., Fleming, P. M., Marques, T., & Costa, F. (2023). Hydrodynamic and hydrodispersive behavior of a highly karstified neoproterozoic hydrosystem indicated by tracer tests and modeling approach. *Journal of Hydrology*, **619**, 129300. https://doi.org/10.1016/j.jhydrol.2023.129300

Atsawawanunt, K., Comas-Bru, L., Amirnezhad Mozhdehi, S., Deininger, M., Harrison, S. P., Baker, A., et al. (2018). The SISAL database: A global resource to document oxygen and carbon isotope records from speleothems. *Earth System Science Data*, **10**(3), 1687–1713. https://doi.org/10.5194/essd-10-1687-2018

Auckenthaler, A., Raso, G., & Huggenberger, P. (2002). Particle transport in a karst aquifer: Natural and artificial tracer experiments with bacteria, bacteriophages and microspheres. *Water Science and Technology*, **46**(3), 131–138. https://doi.org/10.2166/wst.2002.0072

Aydin, H., Ekmekçi, M., & Soylu, M. E. (2014). Effects of sinuosity factor on hydrodynamic parameters estimation in karst systems: A dye tracer experiment from the Beyyayla sinkhole (Eskişehir, Turkey). *Environmental Earth Sciences*, **71**(9), 3921–3933. https://doi.org/10.1007/s12665-013-2777-1

Bailly-Comte, V., Martin, J. B., & Screaton, E. J. (2011). Time variant cross correlation to assess residence time of water and implication for hydraulics of a sink-rise karst system. *Water Resources Research*, **47**(5). https://doi.org/10.1029/2010wr009613

Bajracharya, K., & Barry, D. A. (1997). Nonequilibrium solute transport parameters and their physical significance: Numerical and experimental results. *Journal of Contaminant Hydrology*, **24**(3–4), 185–204. https://doi.org/10.1016/s0169-7722(96)00017-4

Bakalowicz, M. (2005). Karst groundwater: A challenge for new resources. *Hydrogeology Journal*, **13**(1), 148–160. https://doi.org/10.1007/s10040-004-0402-9

Balestra, V., Vigna, B., De Costanzo, S., & Bellopede, R. (2023). Preliminary investigations of microplastic pollution in karst systems, from surface watercourses to cave waters. *Journal of Contaminant Hydrology*, **252**, 104117. https://doi.org/10.1016/j.jconhyd.2022.104117

Banusch, S., Somogyvári, M., Sauter, M., Renard, P., & Engelhardt, I. (2022). Stochastic modeling approach to identify uncertainties of karst conduit networks in carbonate aquifers. *Water Resources Research*, **58**(8), e2021WR031710. https://doi.org/10.1029/2021wr031710

Barberá, J. A., Mudarra, M., Andreo, B., & De la Torre, B. (2018). Regional-scale analysis of karst underground flow deduced from tracing experiments: Examples from carbonate aquifers in Malaga province, southern Spain. *Hydrogeology Journal*, **26**(1), 23–40. https://doi.org/10.1007/s10040-017-1638-5

Barrett, M. E. (1996). *A parsimonious model for simulation of flow and transport in a karst aquifer*. The University of Texas at Austin.

Barrett, M. E., & Charbeneau, R. J. (1997). A parsimonious model for simulating flow in a karst aquifer. *Journal of Hydrology*, **196**(1–4), 47–65. https://doi.org/10.1016/s0022-1694(96)03339-2

Bayari, S. (2002). TRACER: An EXCEL workbook to calculate mean residence time in groundwater by use of tracers CFC-11, CFC-12 and tritium. *Computers & Geosciences*, **28**(5), 621–630. https://doi.org/10.1016/s0098-3004(01)00094-2

Bear, J. (1972). Dynamics of fluids in porous media. *Courier Corporation*.

Bear, J., Tsang, C. F., & De Marsily, G. (2012). *Flow and contaminant transport in fractured rock*. Academic Press.

Becker, M., & Bellin, A. (2013). A reservoir model of tracer transport for karstic flow systems. *Hydrogeology Journal*, **21**(5), 1011–1019. https://doi.org/10.1007/s10040-013-0991-2

Benettin, P., & Bertuzzo, E. (2018). tran-SAS v1.0: A numerical model to compute catchment-scale hydrologic transport using StorAge Selection functions. *Geoscientific Model Development*, **11**(4), 1627–1639. https://doi.org/10.5194/gmd-11-1627-2018

Benettin, P., Rodriguez, N. B., Sprenger, M., Kim, M., Klaus, J., Harman, C. J., et al. (2022). Transit time estimation in catchments: Recent developments and future directions. *Water Resources Research*, **58**(11), e2022WR033096. https://doi.org/10.1029/2022wr033096

Benischke, R. (2021). Advances in the methodology and application of tracing in karst aquifers. *Hydrogeology Journal*, **29**(1), 67–88. https://doi.org/10.1007/s10040-020-02278-9

Benson, D. A., Wheatcraft, S. W., & Meerschaert, M. M. (2000). Application of a fractional advection-dispersion equation. *Water Resources Research*, **36**(6), 1403–1412. https://doi.org/10.1029/2000wr900031

Berglund, J. L., Toran, L., & Herman, E. K. (2020). Can karst conduit models be calibrated? A dual approach using dye tracing and temperature. *Groundwater*, **58**(6), 924–937. https://doi.org/10.1111/gwat.12988

Berkowitz, B. (2002). Characterizing flow and transport in fractured geological media: A review. *Advances in Water Resources*, **25**(8–12), 861–884. https://doi.org/10.1016/s0309-1708(02)00042-8

Berkowitz, B., Cortis, A., Dentz, M., & Scher, H. (2006). Modeling non-Fickian transport in geological formations as a continuous time random walk. *Reviews of Geophysics*, **44**(2). https://doi.org/10.1029/2005rg000178

Berkowitz, B., Dror, I., Hansen, S. K., & Scher, H. (2016). Measurements and models of reactive transport in geological media. *Reviews of Geophysics*, **54**(4), 930–986. https://doi.org/10.1002/2016rg000524

Berkowitz, B., Kosakowski, G., Margolin, G., & Scher, H. (2001). Application of continuous time random walk theory to tracer test measurements in fractured and heterogeneous porous media. *Groundwater*, **39**(4), 593–604. https://doi.org/10.1111/j.1745-6584.2001.tb02347.x

Berkowitz, B., & Scher, H. (1995). On characterization of anomalous dispersion in porous and fractured media. *Water Resources Research*, **31**(6), 1461–1466. https://doi.org/10.1029/95wr00483

Bethke, C. M., & Johnson, T. M. (2008). Groundwater age and groundwater age dating. *Annual Review of Earth and Planetary Sciences*, **36**(1), 121–152. https://doi.org/10.1146/annurev.earth.36.031207.124210

Bettel, L., Fox, J., Husic, A., Zhu, J., Al Aamery, N., Mahoney, T., & Gold-McCoy, A. (2022). Sediment transport investigation in a karst aquifer hypothesizes controls on internal versus external sediment origin and saturation impact on hysteresis. *Journal of Hydrology*, **613**, 128391. https://doi.org/10.1016/j.jhydrol.2022.128391

Beven, K. (2002). Towards a coherent philosophy for modelling the environment. *Proceedings of the royal society of London. Series A: Mathematical, Physical and Engineering Sciences*, **458**(2026), 2465–2484. https://doi.org/10.1098/rspa.2002.0986

Beven, K. (2005). On the concept of model structural error. *Water Science and Technology*, **52**(6), 167–175. https://doi.org/10.2166/wst.2005.0165

Beven, K. (2006). A manifesto for the equifinality thesis. *Journal of Hydrology*, **320**(1–2), 18–36. https://doi.org/10.1016/j.jhydrol.2005.07.007

Beven, K. (2023). Benchmarking hydrological models for an uncertain future. *Hydrological Processes*, **37**(5), e14882. https://doi.org/10.1002/hyp.14882

Beven, K., & Binley, A. (1992). The future of distributed models: Model calibration and uncertainty prediction. *Hydrological Processes*, **6**(3), 279–298. https://doi.org/10.1002/hyp.3360060305

Beven, K., Smith, P. J., & Wood, A. (2011). On the colour and spin of epistemic error (and what we might do about it). *Hydrology and Earth System Sciences*, **15**(10), 3123–3133. https://doi.org/10.5194/hess-15-3123-2011

Beven, K. J. (2011). *Rainfall-runoff modelling: The primer*. John Wiley & Sons.

Beven, K. J., & Chappell, N. A. (2021). Perceptual perplexity and parameter parsimony. *Wiley Interdisciplinary Reviews: Water*, **8**(4), e1530. https://doi.org/10.1002/wat2.1530
Birk, S., Geyer, T., Liedl, R., & Sauter, M. (2005). Process-based interpretation of tracer tests in carbonate aquifers. *Groundwater*, *43*(3), 381–388.
https://doi.org/10.1111/j.1745-6584.2005.0033.x

Birk, S., Liedl, R., & Sauter, M. (2002). Characterisation of gypsum karst aquifers by heat and solute transport simulations. *Acta Geologica Polonica*, *52*(1), 23–29.

Bittner, D., Parente, M. T., Mattis, S., Wohlmuth, B., & Chiogna, G. (2020). Identifying relevant hydrological and catchment properties in active subspaces: An inference study of a lumped karst aquifer model. *Advances in Water Resources*, *135*, 103472. https://doi.org/10.1016/j.advwatres.2019.103472

Bodin, J. (2020). MFIT 1.0: Multi-Flow Inversion of Tracer breakthrough curves in fractured and karst aquifers. *Geoscientific Model Development*, *13*(6), 2905–2924. https://doi.org/10.5194/gmd-13-2905-2020

Bodin, J., Porel, G., Nauleau, B., & Paquet, D. (2022). Delineation of discrete conduit networks in karst aquifers via combined analysis of tracer tests and geophysical data. *Hydrology and Earth System Sciences*, *26*(6), 1713–1726. https://doi.org/10.5194/hess-26-1713-2022

Bonanno, E., Blöschl, G., & Klaus, J. (2022). Exploring tracer information in a small stream to improve parameter identifiability and enhance the process interpretation in transient storage models. *Hydrology and Earth System Sciences*, *26*(23), 6003–6028. https://doi.org/10.5194/hess-26-6003-2022

Borghi, A., Renard, P., & Cornaton, F. (2016). Can one identify karst conduit networks geometry and properties from hydraulic and tracer test data? *Advances in Water Resources*, *90*, 99–115. https://doi.org/10.1016/j.advwatres.2016.02.009

Borghi, A., Renard, P., & Jenni, S. (2012). A pseudo-genetic stochastic model to generate karstic networks. *Journal of Hydrology*, *414*, 516–529.
https://doi.org/10.1016/j.jhydrol.2011.11.032

Bosch, R. F., & White, W. B. (2018). Lithofacies and transport for clastic sediments in karst conduits. In *Karst groundwater contamination and public health: Beyond case studies* (pp. 277–281). Springer International Publishing.

Botter, G. (2012). Catchment mixing processes and travel time distributions. *Water Resources Research*, *48*(5). https://doi.org/10.1029/2011wr011160

Botter, G., Bertuzzo, E., & Rinaldo, A. (2011). Catchment residence and travel time distributions: The master equation. *Geophysical Research Letters*, *38*(11). https://doi.org/10.1029/2011gl047666

Brettmann, K. L., Jensen, K. H., & Jakobsen, R. (1993). Tracer test in fractured chalk: 2. Numerical analysis. *Hydrology Research*, *24*(4), 275–296. https://doi.org/10.2166/nh.1993.0008

Çallı, K. Ö., Bittner, D., & Hartmann, A. (2023). On the choice of a performance metric for model calibration scheme using discharge-age information. In *EuroKarst 2022, Málaga: Advances in the hydrogeology of karst and carbonate reservoirs* (pp. 123–129). Springer International Publishing.

Çallı, K. Ö., Bittner, D., Liu, Y., Çallı, S. S., Melsen, L. A., Bense, V., & Hartmann, A. (2023). Revealing the positive influence of young water fractions derived from stable isotopes on the robustness of karst water resources predictions. *Journal of Hydrology*, *621*, 129549. https://doi.org/10.1016/j.jhydrol.2023.129549

Çallı, S. S., Çallı, K. Ö., Yılmaz, M. T., & Çelik, M. (2022). Contribution of the satellite-data driven snow routine to a karst hydrological model. *Journal of Hydrology*, *607*, 127511. https://doi.org/10.1016/j.jhydrol.2022.127511

Campanale, C., Losacco, D., Triozzi, M., Massarelli, C., & Uricchio, V. F. (2022). An overall perspective for the study of emerging contaminants in karst aquifers. *Resources*, *11*(11), 105. https://doi.org/10.3390/resources11110105

Carrera, J., Sánchez-Vila, X., Benet, I., Medina, A., Galarza, G., & Guimerà, J. (1998). On matrix diffusion: Formulations, solution methods and qualitative effects. *Hydrogeology Journal*, *6*(1), 178–190. https://doi.org/10.1007/s100400050143

Çelik, M., Çallı, S. S., & Karakaş, Z. S. (2022). The role of mineralogical studies in delineating the recharge area and groundwater circulation of Susuz springs, Central Taurus Belt, Turkey. *Hydrogeology Journal*, *30*(8), 2399–2415. https://doi.org/10.1007/s10040-022-02561-x

Cen, X., Xu, M., Qi, J., Zhang, Q., & Shi, H. (2021). Characterization of karst conduits by tracer tests for an artificial recharge scheme. *Hydrogeology Journal*, *29*(7), 2381–2396. https://doi.org/10.1007/s10040-021-02398-w

Chang, Y., Hartmann, A., Liu, L., Jiang, G., & Wu, J. (2021). Identifying more realistic model structures by electrical conductivity observations of the karst spring. *Water Resources Research*, *57*(4), e2020WR028587. https://doi.org/10.1029/2020wr028587

Chang, Y., Wu, J., Jiang, G., & Kang, Z. (2017). Identification of the dominant hydrological process and appropriate model structure of a karst catchment through stepwise simplification of a complex conceptual model. *Journal of Hydrology*, *548*, 75–87. https://doi.org/10.1016/j.jhydrol.2017.02.050

Chang, Y., Wu, J., Jiang, G., Liu, L., Reimann, T., & Sauter, M. (2019). Modelling spring discharge and solute transport in conduits by coupling CFPv2 to an epikarst reservoir for a karst aquifer. *Journal of Hydrology*, *569*, 587–599. https://doi.org/10.1016/j.jhydrol.2018.11.075

Charlier, J. B., Bertrand, C., & Mudry, J. (2012). Conceptual hydrogeological model of flow and transport of dissolved organic carbon in a small Jura karst system. *Journal of Hydrology*, *460*, 52–64. https://doi.org/10.1016/j.jhydrol.2012.06.043

Cherubini, C. (2008). A modeling approach for the study of contamination in a fractured aquifer. *Geotechnical & Geological Engineering*, *26*(5), 519–533. https://doi.org/10.1007/s10706-008-9186-3

Cherubini, C., Giasi, C. I., & Pastore, N. (2013). Evidence of non-Darcy flow and non-Fickian transport in fractured media at laboratory scale. *Hydrology and Earth System Sciences*, *17*(7), 2599–2611. https://doi.org/10.5194/hess-17-2599-2013

Chiogna, G., Eberhardt, C., Grathwohl, P., Cipka, O. A., & Rolle, M. (2010). Evidence of compound-dependent hydrodynamic and mechanical transverse dispersion by multitracer laboratory experiments. *Environmental Science & Technology*, *44*(2), 688–693. https://doi.org/10.1021/es023964

Cholet, C., Charlier, J. B., Moussa, R., Steinmann, M., & Denimal, S. (2017). Assessing lateral flows and solute transport during floods in a conduit-flow-dominated karst system using the inverse problem for the advection–diffusion equation. *Hydrology and Earth System Sciences*, *21*(7), 3635–3653. https://doi.org/10.5194/hess-21-3635-2017

Chu, X., Ding, H., & Zhang, X. (2021). Simulation of solute transport behaviors in saturated karst aquifer system. *Scientific Reports*, *11*(1), 1–17. https://doi.org/10.1038/s41598-021-94950-7

Cinkus, G., Wunsch, A., Mazzilli, N., Liesch, T., Chen, Z., Ravbar, N., et al. (2023). Comparison of artificial neural networks and reservoir models for simulating karst spring discharge on five test sites in the Alpine and Mediterranean regions. *Hydrology and Earth System Sciences*, *27*(10), 1961–1985. https://doi.org/10.5194/hess-27-1961-2023

Coats, K. H., & Smith, B. D. (1964). Dead-end pore volume and dispersion in porous media. *Society of Petroleum Engineers Journal*, *40*(1), 73–84. https://doi.org/10.2118/647-pa

Comas-Bru, L., Rehfeld, K., Roesch, C., Amirnezhad-Mozhdehi, S., Harrison, S. P., Atsawawaranunt, K., SISAL Working Group members, et al. (2020). SISALv2: A comprehensive speleothem isotope database with multiple age–depth models. *Earth System Science Data*, *12*(4), 2579–2606. https://doi.org/10.5194/essd-12-2579-2020
Cornaton, F., & Perrochet, P. (2006a). Groundwater age, life expectancy and transit time distributions in advective-dispersive systems: 1. Generalized reservoir theory. *Advances in Water Resources*, 29(9), 1267-1291. https://doi.org/10.1016/j.adwatres.2005.10.009

Cornaton, F., & Perrochet, P. (2006b). Groundwater age, life expectancy and transit time distributions in advective-dispersive systems: 2. Reservoir theory for sub-drainage basins. *Advances in Water Resources*, 29(9), 1292-1305. https://doi.org/10.1016/j.adwatres.2005.10.010

Cousquer, Y., & Jourde, H. (2022). Reducing uncertainty of karst aquifer modeling with complementary hydrological observations for the sustainable management of groundwater resources. *Journal of Hydrology*, 612, 128130. https://doi.org/10.1016/j.jhydrol.2022.128130

Covington, M. D., Doctor, D. H., King, J. N., & Wicks, C. M. (2011). Research in karst: A model for future directions in hydrologic science? *From the section President, 18*.

Covington, M. D., Luhmann, A. J., Wicks, C. M., & Saar, M. O. (2012). Process length scales and longitudinal damping in karst conduits. *Journal of Geophysical Research*, 117(F1). https://doi.org/10.1029/2011jf002212

Covington, M. D., Wicks, C. M., & Saar, M. O. (2009). A dimensionless number describing the effects of recharge and geometry on discharge from simple karstic aquifers. *Water Resources Research*, 45(11). https://doi.org/10.1029/2009wr008004

Dal Soglio, L., Danquigny, C., Mazzilli, N., Emblanch, C., & Massonnat, G. (2020). Modeling the matrix-conduit exchanges in both the epikarst and the transmission zone of karst systems. *Water, 12*(11), 3219. https://doi.org/10.3390/w12113219

Deleu, R., Poulain, A., Rochez, G., Soares-Frazao, S., Van Rentergem, G., De Poorter, E., & Hallet, V. (2023). Multipeaked breakthrough curves in karstic rivers: Effects of a diffuence-confluence system. *Hydrogeology Journal*, 32(2), 1-19. https://doi.org/10.1007/s10040-023-02728-0

Dentz, M., & Berkowitz, B. (2003). Transport behavior of a passive solute in continuous time random walks and multirate mass transfer. *Water Resources Research*, 39(5). https://doi.org/10.1029/2001wr001163

Dentz, M., Cortis, A., Scher, H., & Berkowitz, B. (2004). Time behavior of solute transport in heterogeneous media: Transition from anomalous to normal transport. *Advances in Water Resources*, 27(2), 155-173. https://doi.org/10.1016/j.adwatres.2003.11.002

Dentz, M., Hidalgo, J. J., & Lester, D. (2023). Mixing in porous media: Concepts and approaches across scales. *Transport in Porous Media*, 146(1-2), 5-53. https://doi.org/10.1007/s11242-022-01852-x

Dentz, M., Le Borgne, T., Englert, A., & Bijeljic, B. (2011). Mixing, spreading and reaction in heterogeneous media: A brief review. *Journal of Contaminant Hydrology*, 120, 1-17. https://doi.org/10.1016/j.jconhyd.2010.05.002

Dewaide, L., Bonniver, I., Rochez, G., & Hallet, V. (2016). Solute transport in heterogeneous karst systems: Dimensioning and estimation of the transport parameters via multi-sampling tracer-tests modelling using the OTIS (One-dimensional transport with inflow and storage) program. *Journal of Hydrology*, 534, 567-578. https://doi.org/10.1016/j.jhydrol.2016.01.049

Dewaide, L., Collon, P., Poulain, A., Rochez, G., & Hallet, V. (2018). Double-peaked breakthrough curves as a consequence of solute transport through underground lakes: A case study of the Furfooz karst system, Belgium. *Hydrogeology Journal*, 26(2), 641-650. https://doi.org/10.1007/s10040-017-1671-4

Dietrich, J., & Uniyal, B. (2020). Modeling interbasin groundwater flow in karst areas: Model development, application, and calibration strategy. *Environmental Modelling & Software*, 124, 104606. https://doi.org/10.1016/j.envsoft.2019.104606

Dong, P., Yin, M., Zhang, Y., Chen, K., Finkel, M., Grathwohl, P., & Zheng, C. (2023). A Fractional-order dual-continuum model to capture non-Fickian solute transport in a regional-scale fractured aquifer. *Journal of Contaminant Hydrology*, 258, 104231. https://doi.org/10.1016/j.jconhyd.2023.104231

Doummar, J., Margane, A., Geyer, T., & Sauter, M. (2018). Assessment of key transport parameters in a karst system under different dynamic conditions based on tracer experiments: The Jeita karst system, Lebanon. *Hydrogeology Journal*, 26(7), 2283-2295. https://doi.org/10.1007/s10040-018-1754-x

Dreybrodt, W. (1996). Principles of early development of karst conduits under natural and man-made conditions revealed by mathematical analysis of numerical models. *Water Resources Research*, 32(9), 2923-2935. https://doi.org/10.1029/96wr01332

Dreybrodt, W. (2012). *Processes in karst systems: Physics, chemistry, and geology* (Vol. 4). Springer Science & Business Media.

Dreybrodt, W., & Gabrovšek, F. (2018). Dynamics of wormhole formation in fractured karst aquifers. In *Hydrology and Earth system sciences discussion*. https://doi.org/10.5194/hess-2018-275

Druhan, J. L., Lawrence, C. R., Covey, A. K., Giannetta, M. G., & Oster, J. L. (2021). A reactive transport approach to modeling cave seepage water chemistry I: Carbon isotope transformations. *Geochimica et Cosmochimica Acta*, 311, 374-400. https://doi.org/10.1016/j.gca.2021.06.041

Dubois, C., Quinif, Y., Baele, J. M., Barriquand, L., Bini, A., Bruxelles, L., et al. (2014). The process of ghost-rock karstification and its role in the formation of cave systems. *Earth-Science Reviews*, 131, 116-148. https://doi.org/10.1016/j.earscirev.2014.01.006

Duran, L., Fournier, M., Massei, N., & Dupont, J. P. (2016). Assessing the nonlinearity of karst response function under variable boundary conditions. *Groundwater*, 54(1), 46-54. https://doi.org/10.1111/gwat.12337

Dussart-Baptista, L., Massei, N., Dupont, J. P., & Jouenne, T. (2003). Transfer of bacteria-contaminated particles in a karst aquifer: Evolution of contaminated materials from a sinkhole to a spring. *Journal of Hydrology*, 284(1-4), 285-295. https://doi.org/10.1016/j.jhydrol.2003.08.007

Dvory, N. Z., Kuznetsov, M., Livshitz, Y., Gasser, G., Pankratov, I., Lev, O., et al. (2018). Modeling sewage leakage and transport in carbonate aquifer using carbamazepine as an indicator. *Water Research*, 128, 157-170. https://doi.org/10.1016/j.watres.2017.10.044

Dvory, N. Z., Livshitz, Y., Kuznetsov, M., Adar, E., Gasser, G., Pankratov, I., et al. (2018). Caffeine vs. carbamazepine as indicators of wastewater pollution in a karst aquifer. *Hydrology and Earth System Sciences*, 22(12), 6371-6381. https://doi.org/10.5194/hess-22-6371-2018

Einsiedl, F., Maloszewski, P., & Stichler, W. (2009). Multiple isotope approach to the determination of the natural attenuation potential of a high-alpine karst system. *Journal of Hydrology*, 365(1-2), 113-121. https://doi.org/10.1016/j.jhydrol.2008.11.042

Ender, A., Goeppert, N., & Goldscheider, N. (2018). Spatial resolution of transport parameters in a subtropical karst conduit system during dry and wet seasons. *Hydrogeology Journal*, 26(7), 2241-2255. https://doi.org/10.1007/s10040-018-1746-x

Enemark, T., Peeters, L. J., Mallants, D., & Batelaan, O. (2019). Hydrogeological conceptual model building and testing: A review. *Journal of Hydrology*, 569, 310-329. https://doi.org/10.1016/j.jhydrol.2018.12.007

Ewers, R. O., White, K. A., & Fuller, J. F. (2012). Contaminant plumes and pseudoplumes in karst aquifers. *Carbonates and Evaporites*, 27(2), 153-159. https://doi.org/10.1007/s13146-012-0099-0

Fandel, C., Miville, F., Ferré, T., Goldscheider, N., & Renard, P. (2022). The stochastic simulation of karst conduit network structure using anisotropic fast marching, and its application to a geologically complex alpine karst system. *Hydrogeology Journal*, 30(3), 927-946. https://doi.org/10.1007/s10040-022-02464-x

Faulkner, J., Hu, B. X., Kish, S., & Hua, F. (2009). Laboratory analog and numerical study of groundwater flow and solute transport in a karst aquifer with conduit and matrix domains. *Journal of Contaminant Hydrology*, 110(1-2), 34-44. https://doi.org/10.1016/j.jconhyd.2009.08.004

Feehley, C. E., Zheng, C., & Molz, F. J. (2000). A dual-domain mass transfer approach for modeling solute transport in heterogeneous aquifers: Application to the Macrodispersion Experiment (MADE) site. *Water Resources Research*, 36(9), 2501-2515. https://doi.org/10.1029/2000wr900148
Fenicia, F., Wrede, S., Kavetski, D., Pfister, L., Hoffmann, L., Savenije, H. H., & McDonnell, J. J. (2010). Assessing the impact of mixing assumptions on the estimation of streamwater mean residence time. *Hydrological Processes*, 24(12), 1730–1741. https://doi.org/10.1002/hyp.7595

Field, M. S. (2002). *The QTRACER2 program for tracer-breakthrough curve analysis for tracer tests in karstic aquifers and other hydrologic systems*. National Center for Environmental Assessment–Washington Office, Office of Research and Development, US Environmental Protection Agency.

Field, M. S. (2018). Investigating and remediating contaminated karst aquifers. In *Karst groundwater contamination and public health: Beyond case studies* (pp. 101–115). Springer International Publishing.

Field, M. S. (2020). Groundwater sampling in karst terranes: Passive sampling in comparison to event-driven sampling strategy. *Hydrogeology Journal*, 29(1), 53–65. https://doi.org/10.1007/s10040-020-02240-9

Field, M. S., & Leij, F. J. (2012). Solute transport in solution conduits exhibiting multi-peaked breakthrough curves. *Journal of Hydrology*, 440, 26–35. https://doi.org/10.1016/j.jhydrol.2012.03.018

Field, M. S., & Leij, F. J. (2014). Combined physical and chemical nonequilibrium transport model for solution conduits. *Journal of Contaminant Hydrology*, 157, 37–46. https://doi.org/10.1016/j.jconhyd.2013.11.001

Field, M. S., & Li, G. (2011). Inversion for the input history of a dye tracing experiment. *Journal of Cave and Karst Studies*, 73(1), 16–20. https://doi.org/10.4311/jcks2010es0143

Field, M. S., & Nash, S. G. (1997). Risk assessment methodology for karst aquifers: (1) Estimating karst conduit-flow parameters. *Environmental Monitoring and Assessment*, 47, 1–21. https://doi.org/10.1023/a:1005753919403

Field, M. S., & Pinsky, P. F. (2000). A two-region nonequilibrium model for solute transport in solution conduits in karstic aquifers. *Journal of Contaminant Hydrology*, 44(3–4), 329–351. https://doi.org/10.1016/s0169-7722(00)00009-1

Finsterle, S., Kowalsky, M. B., & Pruess, K. (2012). TOUGH: Model use, calibration, and validation. *Transactions of the ASABE*, 55(4), 1275–1290. https://doi.org/10.13031/2013.42240

Fischer, H. B., List, J. E., Koh, C. R., Imberger, J., & Brooks, N. H. (1979). *Mixing in inland and coastal waters*. Elsevier.

Fischer, P., Jourde, H., Brunet, P., & Leonardi, V. (2024). Large-scale spatial reconstitution of groundwater fluxes in a karst aquifer on the basis of hydrodynamic and hydrodispersive measurements. *Journal of Hydrology*, 632, 130866. https://doi.org/10.1016/j.jhydrol.2024.130866

Florea, L. J., & Wicks, C. M. (2001). Solute transport through laboratory-scale karstic aquifers. *Journal of Cave and Karst Studies*, 63(2), 59–66.

Flynn, R. M., & Sinreich, M. (2010). Characterisation of virus transport and attenuation in epikarst using short pulse and prolonged injection multi-tracer testing. *Water Research*, 44(4), 1138–1149. https://doi.org/10.1016/j.watres.2009.11.032

Ford, D., & Williams, P. D. (2007). *Karst hydrogeology and geomorphology*. John Wiley & Sons.

Frank, S., Goeppert, N., & Goldscheider, N. (2021). Field tracer tests to evaluate transport properties of tryptophan and humic acid in karst. *Groundwater*, 59(1), 59–70. https://doi.org/10.1111/gwat.13015

Gelhar, L. W., Welty, C., & Rehfeldt, K. R. (1992). A critical review of data on field-scale dispersion in aquifers. *Water Resources Research*, 28(7), 1955–1974. https://doi.org/10.1029/92wr00607

Gerke, H. H., & Van Genuchten, M. T. (1993). A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media. *Water Resources Research*, 29(2), 305–319. https://doi.org/10.1029/92wr02339

Geyer, T. (2008). *Process-based characterisation of flow and transport in karst aquifers at catchment scale*. (Doctoral dissertation). Göttingen: Georg-August-Universität. Retrieved from https://ediss.uni-goettingen.de/handle/11858/00-1735-0000-0006-B276-E?locale-attribute=de

Geyer, T., Birk, S., Licha, T., Liedl, R., & Sauter, M. (2007). Multitracer test approach to characterize reactive transport in karst aquifers. *Groundwater*, 45(1), 36–45. https://doi.org/10.1111/j.1745-6584.2006.00261.x

Geyer, T., Birk, S., Liedl, R., & Sauter, M. (2008). Quantification of temporal distribution of recharge in karst systems from spring hydrographs. *Journal of Hydrology*, 348(3–4), 452–463. https://doi.org/10.1016/j.jhydrol.2007.10.015

Geyer, T., Birk, S., Reimann, T., Dörfliger, N., & Sauter, M. (2013). Differentiated characterization of karst aquifers: Some contributions. *Carbonates and Evaporites*, 28(1–2), 41–46. https://doi.org/10.1007/s13146-013-0150-9

Ghanbarian, B., Hunt, A. G., Ewing, R. P., & Sahimi, M. (2013). Tortuosity in porous media: A critical review. *Soil Science Society of America Journal*, 77(5), 1461–1477. https://doi.org/10.2136/sssj2012.0435

Gharari, S., Gupta, H. V., Clark, M. P., Hrachowitz, M., Fenicia, F., Matgen, P., & Savenije, H. H. (2021). Understanding the information content in the hierarchy of model development decisions: Learning from data. *Water Resources Research*, 57(6), e2020WR027948. https://doi.org/10.1029/2020wr027948

Ghasemizadeh, R., Hellweger, F., Butscher, C., Padilla, I., Vesper, D., Field, M., & Alshawabkeh, A. (2012). Groundwater flow and transport modeling of karst aquifers, with particular reference to the North Coast Limestone aquifer system of Puerto Rico. *Hydrogeology Journal*, 20(8), 1441–1461. https://doi.org/10.1007/s10040-012-0897-4

Ghasemizadeh, R., Yu, X., Butscher, C., Padilla, I. Y., & Alshawabkeh, A. (2016). Improved regional groundwater flow modeling using drainage features: A case study of the central northern karst aquifer system of Puerto Rico (USA). *Hydrogeology Journal*, 24(6), 1463–1478. https://doi.org/10.1007/s10040-016-1419-6

Giese, M., Reimann, T., Bailly-Comte, V., Maréchal, J. C., Sauter, M., & Geyer, T. (2018). Turbulent and laminar flow in karst conduits under unsteady flow conditions: Interpretation of pumping tests by discrete conduit-continuum modeling. *Water Resources Research*, 54(3), 1918–1933. https://doi.org/10.1002/2017wr020658

Gill, L. W., Naughton, O., & Johnston, P. M. (2013). Modeling a network of turloughs in lowland karst. *Water Resources Research*, 49(6), 3487–3503. https://doi.org/10.1002/wrcr.20299

Gill, L. W., Schuler, P., Duran, L., Morrissey, P., & Johnston, P. M. (2021). An evaluation of semidistributed-pipe-network and distributed-finite-difference models to simulate karst systems. *Hydrogeology Journal*, 29(1), 259–279. https://doi.org/10.1007/s10040-020-02241-8

Ginn, T. R. (1999). On the distribution of multicomponent mixtures over generalized exposure time in subsurface flow and reactive transport: Foundations, and formulations for groundwater age, chemical heterogeneity, and biodegradation. *Water Resources Research*, 35(5), 1395–1407. https://doi.org/10.1029/1999wr900013

Ginn, T. R., Haeri, H., Massoudieh, A., & Foglia, L. (2009). Notes on groundwater age in forward and inverse modeling. *Transport in Porous Media*, 79(1), 117–134. https://doi.org/10.1007/s11242-009-9406-1

Goeppert, N., Goldscheider, N., & Berkowitz, B. (2020). Experimental and modeling evidence of kilometer-scale anomalous tracer transport in an alpine karst aquifer. *Water Research*, 178, 115755. https://doi.org/10.1016/j.watres.2020.115755

Goldscheider, N. (2005). Karst groundwater vulnerability mapping: Application of a new method in the Swabian Alb, Germany. *Hydrogeology Journal*, 13(4), 555–564. https://doi.org/10.1007/s10040-003-0291-3

Goldscheider, N. (2019). A holistic approach to groundwater protection and ecosystem services in karst terrains. *Carbonates and Evaporites*, 34(4), 1241–1249. https://doi.org/10.1007/s13146-019-00492-5