# A Review of the Application of the Soil and Water Assessment Tool (SWAT) in Karst Watersheds

Ibrahim Al Khoury, Laurie Boithias, David Labat

► To cite this version:

Ibrahim Al Khoury, Laurie Boithias, David Labat. A Review of the Application of the Soil and Water Assessment Tool (SWAT) in Karst Watersheds. Water, 2023, 15 (5), pp.954. 10.3390/w15050954. hal-04584215

HAL Id: hal-04584215

https://ut3-toulouseinp.hal.science/hal-04584215

Submitted on 23 May 2024

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
# A Review of the Application of the Soil and Water Assessment Tool (SWAT) in Karst Watersheds

Ibrahim Al Khoury, Laurie Boithias, David Labat

► To cite this version:

Ibrahim Al Khoury, Laurie Boithias, David Labat. A Review of the Application of the Soil and Water Assessment Tool (SWAT) in Karst Watersheds. Water, 2023, 15 (5), pp.954. 10.3390/w15050954. hal-04584215

HAL Id: hal-04584215

https://ut3-toulouseinp.hal.science/hal-04584215v1

Submitted on 23 May 2024

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
<figure><img src="image_2.png" /></figure>

Review

# A Review of the Application of the Soil and Water Assessment Tool (SWAT) in Karst Watersheds

Ibrahim Al Khoury *, Laurie Boithias and David Labat

GET, Université de Toulouse, CNRS, IRD, UPS, 31400 Toulouse, France

* Correspondence: ibrahim.al-khoury@get.omp.eu

**Abstract:** Karst water resources represent a primary source of freshwater supply, accounting for nearly 25% of the global population water needs. Karst aquifers have complex recharge characteristics, storage patterns, and flow dynamics. They also face a looming stress of depletion and quality degradation due to natural and anthropogenic pressures. This prompted hydrogeologists to apply innovative numerical approaches to better understand the functioning of karst watersheds and support karst water resources management. The Soil and Water Assessment Tool (SWAT) is a semi-distributed hydrological model that has been used to simulate flow and water pollutant transport, among other applications, in basins including karst watersheds. Its source code has also been modified by adding distinctive karst features and subsurface hydrology models to more accurately represent the karst aquifer discharge components. This review summarizes and discusses the findings of 75 SWAT-based studies in watersheds that are at least partially characterized by karst geology, with a primary focus on the hydrological assessment in modified SWAT models. Different karst processes were successfully implemented in SWAT, including the recharge in the epikarst, flows of the conduit and matrix systems, interbasin groundwater flow, and allogenic recharge from sinkholes and sinking streams. Nonetheless, additional improvements to the existing SWAT codes are still needed to better reproduce the heterogeneity and non-linearity of karst flow and storage mechanisms in future research.

**Keywords:** karst hydrology; SWAT model; modified SWAT model; hydrological modeling

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

check for
updates

**Citation:** Al Khoury, I.; Boithias, L.; Labat, D. A Review of the Application of the Soil and Water Assessment Tool (SWAT) in Karst Watersheds. *Water* **2023**, *15*, 954. https://doi.org/10.3390/w15050954

Academic Editor: Renato Morbidelli

Received: 2 February 2023

Revised: 25 February 2023

Accepted: 27 February 2023

Published: 1 March 2023

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

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Karst aquifers are an abundant source of water in many regions across the globe, providing freshwater supply to 20–25% of the world population [1] and upwards of 50% of the total drinking water supply in some countries [2]. They cover nearly 15.2% of Earth’s continental surface [3] and form by chemical dissolution of soluble carbonate rocks (i.e., limestone, dolomite, marble or evaporates) exerted by water enriched with carbon dioxide (CO₂) from the atmosphere or soil zone [4]. Depending on the degree of karstification, distinctive karst features can develop, including sinkholes and dolines, losing streams, springs, and vast networks of subsurface and hydrologically connected cracks, fissures, conduits, and caves [5].

### 1.1. Characteristics of Karst Systems

A karst system is generally composed of four main water-bearing mediums with distinct geomorphology, hydrodynamic properties, storage, and flow patterns: (1) the soil and non-karstic zone, (2) the epikarst, (3) the transmission zone—the latter three forming the unsaturated zone, and (4) the saturated zone [6]. These contrasting layers, which are interactively connected by water flow and solute transport, form the karstic critical zone [7,8].

Figure 1 shows a schematic model of a typical karst aquifer, including the surface hydrological processes and flow mechanisms of the underground karst subsystems. The
<figure><img src="image_2.png" /></figure>

Review

# A Review of the Application of the Soil and Water Assessment Tool (SWAT) in Karst Watersheds

Ibrahim Al Khoury *, Laurie Boithias and David Labat

GET, Université de Toulouse, CNRS, IRD, UPS, 31400 Toulouse, France

* Correspondence: ibrahim.al-khoury@get omp.eu

**Abstract:** Karst water resources represent a primary source of freshwater supply, accounting for nearly 25% of the global population water needs. Karst aquifers have complex recharge characteristics, storage patterns, and flow dynamics. They also face a looming stress of depletion and quality degradation due to natural and anthropogenic pressures. This prompted hydrogeologists to apply innovative numerical approaches to better understand the functioning of karst watersheds and support karst water resources management. The Soil and Water Assessment Tool (SWAT) is a semi- distributed hydrological model that has been used to simulate flow and water pollutant transport, among other applications, in basins including karst watersheds. Its source code has also been modified by adding distinctive karst features and subsurface hydrology models to more accurately represent the karst aquifer discharge components. This review summarizes and discusses the findings of 75 SWAT-based studies in watersheds that are at least partially characterized by karst geology, with a primary focus on the hydrological assessment in modified SWAT models. Different karst processes were successfully implemented in SWAT, including the recharge in the epikarst, flows of the conduit and matrix systems, interbasin groundwater flow, and allogenic recharge from sinkholes and sinking streams. Nonetheless, additional improvements to the existing SWAT codes are still needed to better reproduce the heterogeneity and non-linearity of karst flow and storage mechanisms in future research.

**Keywords:** karst hydrology; SWAT model; modified SWAT model; hydrological modeling

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

check for
updates

**Citation:** Al Khoury, I.; Boithias, L.; Labat, D. A Review of the Application of the Soil and Water Assessment Tool (SWAT) in Karst Watersheds. *Water* **2023**, *15*, 954. https://doi.org/10.3390/w15050954

Academic Editor: Renato Morbidelli

Received: 2 February 2023

Revised: 25 February 2023

Accepted: 27 February 2023

Published: 1 March 2023

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

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Karst aquifers are an abundant source of water in many regions across the globe, providing freshwater supply to 20–25% of the world population [1] and upwards of 50% of the total drinking water supply in some countries [2]. They cover nearly 15.2% of Earth’s continental surface [3] and form by chemical dissolution of soluble carbonate rocks (i.e., limestone, dolomite, marble or evaporates) exerted by water enriched with carbon dioxide (CO₂) from the atmosphere or soil zone [4]. Depending on the degree of karstification, distinctive karst features can develop, including sinkholes and dolines, losing streams, springs, and vast networks of subsurface and hydrologically connected cracks, fissures, conduits, and caves [5].

### 1.1. Characteristics of Karst Systems

A karst system is generally composed of four main water-bearing mediums with distinct geomorphology, hydrodynamic properties, storage, and flow patterns: (1) the soil and non-karstic zone, (2) the epikarst, (3) the transmission zone—the latter three forming the unsaturated zone, and (4) the saturated zone [6]. These contrasting layers, which are interactively connected by water flow and solute transport, form the karstic critical zone [7,8].

Figure 1 shows a schematic model of a typical karst aquifer, including the surface hydrological processes and flow mechanisms of the underground karst subsystems. The
epikarst represents a weathered horizon of a few meters above the vadose zone, with high permeability and porosity driven by the large supply of CO₂ that increases dissolution of carbonate rocks near the land surface. The epikarst, together with the soil cover, controls water infiltration, storage, temporal delay of recharge, and mixing processes. On the other hand, seasonal changes in surface temperature and vegetation cover substantially alter evapotranspiration and recharge, thus intensifying the variability of spring discharge and affecting the quality of the underground water resources [9–11]. Two recharge mechanisms are generally observed in a karst system: (1) diffuse recharge by slow percolation of infiltrated water from the epikarst to the saturated zone through low permeability small fissures in the vadose zone and (2) concentrated recharge via highly conductive karst features (enlarged fractures, sinkholes), allowing a fast transit of flow through the vadose zone to the saturated zone [12]. The transmission zone connects and transfers recharge water from the epikarst to the saturated zone where highly permeable karst conduits drain the fissured rock matrix, generating a flow to the groundwater discharge. Karst systems thus exhibit a duality of the storage and subsurface flow fields with (1) prolonged groundwater storage and low flow velocity/laminar flow in the matrix, and (2) low groundwater storage with rapid flow velocity/non-linear (turbulent) flow in the conduits. Dual discharge patterns to the aquifer outlet are also observed with (1) slow and continuous flow from the matrix during dry periods and (2) fast flow from the conduits during heavy rainfall events [6,13–15]. There is also limited karst hydrogeological research on the flow exchange mechanism between the conduit and the matrix, which primarily depends on the conduit properties and hydraulic head differences between the two mediums [16]. Internal runoff and infiltration provide autogenic recharge, whereas external runoff and sinking streams represent allogenic recharge from the neighboring areas [13]. Moreover, many karst watersheds are non-conservative (losing or gaining watersheds) due to inter-basin groundwater flow (IGF) through their topographic divides. IGF may represent a substantial component of the streamflow where rivers traverse karst areas, thus affecting the catchment annual water balance. IGF is also not explicitly measurable, requiring the application of hydrogeochemical approaches based on major dissolved elements, isotopes, electrical conductivity, and water temperature monitoring or hydrogeological studies of groundwater flow paths [17].

Due to their intrinsic properties and complex hydrodynamic behavior, karst aquifers are vulnerable to contamination, overexploitation, and climate change [18,19]. In well-developed karst systems, natural processes such as absorption, degradation, and filtration are inefficient due to low storage capacity, fast water movement, short residence time, and limited interaction with the material of the aquifer. Contaminants can rapidly reach the groundwater table by concentrated recharge and propagate easily through karst conduits over large distances [20,21]. Moreover, climate change could result in more frequent and extended periods of high/transit or low groundwater levels [22]. Therefore, anticipating the impacts of climate change and anthropogenic hazards, and understanding the functioning of the karst aquifer water bearing components are compulsory tasks to safeguard these dwindling water supplies and set effective management schemes for karst water resources [23,24].
epikarst represents a weathered horizon of a few meters above the vadose zone, with high permeability and porosity driven by the large supply of CO₂ that increases dissolution of carbonate rocks near the land surface. The epikarst, together with the soil cover, controls water infiltration, storage, temporal delay of recharge, and mixing processes. On the other hand, seasonal changes in surface temperature and vegetation cover substantially alter evapotranspiration and recharge, thus intensifying the variability of spring discharge and affecting the quality of the underground water resources [9–11]. Two recharge mechanisms are generally observed in a karst system: (1) diffuse recharge by slow percolation of infiltrated water from the epikarst to the saturated zone through low permeability small fissures in the vadose zone and (2) concentrated recharge via highly conductive karst features (enlarged fractures, sinkholes), allowing a fast transit of flow through the vadose zone to the saturated zone [12]. The transmission zone connects and transfers recharge water from the epikarst to the saturated zone where highly permeable karst conduits drain the fissured rock matrix, generating a flow to the groundwater discharge. Karst systems thus exhibit a duality of the storage and subsurface flow fields with (1) prolonged groundwater storage and low flow velocity/laminar flow in the matrix, and (2) low groundwater storage with rapid flow velocity/non-linear (turbulent) flow in the conduits. Dual discharge patterns to the aquifer outlet are also observed with (1) slow and continuous flow from the matrix during dry periods and (2) fast flow from the conduits during heavy rainfall events [6,13–15]. There is also limited karst hydrogeological research on the flow exchange mechanism between the conduit and the matrix, which primarily depends on the conduit properties and hydraulic head differences between the two mediums [16]. Internal runoff and infiltration provide autogenic recharge, whereas external runoff and sinking streams represent allogenic recharge from the neighboring areas [13]. Moreover, many karst watersheds are non-conservative (losing or gaining watersheds) due to inter-basin groundwater flow (IGF) through their topographic divides. IGF may represent a substantial component of the streamflow where rivers traverse karst areas, thus affecting the catchment annual water balance. IGF is also not explicitly measurable, requiring the application of hydrogeochemical approaches based on major dissolved elements, isotopes, electrical conductivity, and water temperature monitoring or hydrogeological studies of groundwater flow paths [17].

Due to their intrinsic properties and complex hydrodynamic behavior, karst aquifers are vulnerable to contamination, overexploitation, and climate change [18,19]. In well-developed karst systems, natural processes such as absorption, degradation, and filtration are inefficient due to low storage capacity, fast water movement, short residence time, and limited interaction with the material of the aquifer. Contaminants can rapidly reach the groundwater table by concentrated recharge and propagate easily through karst conduits over large distances [20,21]. Moreover, climate change could result in more frequent and extended periods of high/transit or low groundwater levels [22]. Therefore, anticipating the impacts of climate change and anthropogenic hazards, and understanding the functioning of the karst aquifer water bearing components are compulsory tasks to safeguard these dwindling water supplies and set effective management schemes for karst water resources [23,24].
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Conceptual schematic of a karst aquifer illustrating the heterogeneous hydrological behavior of a karst system (epikarst, matrix, and conduits), with dual infiltration and recharge processes, dual subsurface flow fields, and dual discharge characteristics. Karst aquifers are a primary source of freshwater supply for residential, industrial, and agricultural uses and are highly vulnerable to climate change and anthropogenic hazards.</figcaption></figure>

## 1.2. *Numerical Modeling Approaches in Karst Hydrology*

Hydrological models have become robust tools to simulate karst aquifer processes for a wide range of applications, which include integrated hydrodynamic analysis, modeling of karst hydrology, and forecasting of karst water resources availability under climate change [25]. These models can be classified into the following categories: the black box models, lumped models, and semi-distributed to distributed models.

Black box models use mathematical transfer functions or neural networks to relate the input rainfall signal to the output spring discharge without spatial information or an explicit representation of the watershed physical processes [26]. Thus, they may not be adequate to estimate future karst water resources over long time periods [23].

Lumped models conceptualize the physical processes at the scale of the entire hydrological system. They are generally used by modelers facing data scarcity problems [27]. Most lumped groundwater models adopt a series of linear or non-linear reservoirs to simulate the storage and flow components of the karst aquifer mediums, with parameters that represent the spatially averaged characteristics of the system [23].

In comparison to lumped models, the semi-distributed and distributed models explicitly represent the spatial variability of the watershed land and subsurface characteristics, boundary conditions, inputs, and hydrological processes [27]. Semi-distributed models could divide the catchment into hydrological response units (HRUs) and simulate the various hydrological processes in each HRU, use conceptual reservoirs to model areal recharge processes that lack spatial resolution, or represent the internal structure of a karst aquifer using pipe networks as conduit domains. On the other hand, fully-distributed models represent processes by discretizing the system in two- or three-dimensional grids and assigning parameters to each grid cell. In terms of karst hydrological modeling, the distributed karst models are subdivided into three categories: (1) the fully equivalent porous media approach, which uses average hydraulic properties over the aquifer area
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Conceptual schematic of a karst aquifer illustrating the heterogeneous hydrological behavior of a karst system (epikarst, matrix, and conduits), with dual infiltration and recharge processes, dual subsurface flow fields, and dual discharge characteristics. Karst aquifers are a primary source of freshwater supply for residential, industrial, and agricultural uses and are highly vulnerable to climate change and anthropogenic hazards.</figcaption></figure>

## 1.2. *Numerical Modeling Approaches in Karst Hydrology*

Hydrological models have become robust tools to simulate karst aquifer processes for a wide range of applications, which include integrated hydrodynamic analysis, modeling of karst hydrology, and forecasting of karst water resources availability under climate change [25]. These models can be classified into the following categories: the black box models, lumped models, and semi-distributed to distributed models.

Black box models use mathematical transfer functions or neural networks to relate the input rainfall signal to the output spring discharge without spatial information or an explicit representation of the watershed physical processes [26]. Thus, they may not be adequate to estimate future karst water resources over long time periods [23].

Lumped models conceptualize the physical processes at the scale of the entire hydrological system. They are generally used by modelers facing data scarcity problems [27]. Most lumped groundwater models adopt a series of linear or non-linear reservoirs to simulate the storage and flow components of the karst aquifer mediums, with parameters that represent the spatially averaged characteristics of the system [23].

In comparison to lumped models, the semi-distributed and distributed models explicitly represent the spatial variability of the watershed land and subsurface characteristics, boundary conditions, inputs, and hydrological processes [27]. Semi-distributed models could divide the catchment into hydrological response units (HRUs) and simulate the various hydrological processes in each HRU, use conceptual reservoirs to model areal recharge processes that lack spatial resolution, or represent the internal structure of a karst aquifer using pipe networks as conduit domains. On the other hand, fully-distributed models represent processes by discretizing the system in two- or three-dimensional grids and assigning parameters to each grid cell. In terms of karst hydrological modeling, the distributed karst models are subdivided into three categories: (1) the fully equivalent porous media approach, which uses average hydraulic properties over the aquifer area
(concentrated fast flows and diffuse slow flows are not explicitly simulated), (2) the double continuum approach, which represents the matrix and karst conduits as two interacting continua with their hydraulic attributes, and (3) the combined discrete-continuum approach in which the karst conduits are embedded as discrete elements inside the matrix [26,27]. The use of distributed parameters models remains a challenge due to the complexity of karst aquifer mechanisms and the need for extensive hydrological and hydrogeological investigations to define the characteristics of the karst system (i.e., aquifer geometry, conduits network geometry and location, hydraulic properties and interactions between the matrix and conduits) [23,26].

### 1.3. Rationale of This Review of Karst Hydrological Modeling with SWAT

Improving the management and sustainability of karst groundwater resources remains a challenge. This is due to the limited understanding of the critical zone processes in karst watersheds across space and time, as well as the lack of research that characterizes the influence of vegetation cover, climate change, and anthropogenic activities on these processes [8,9]. Jeannin et al. [28] recently tested 13 karst numerical models in a karst watershed. The models consisted of lumped neural networks, reservoir models, and semi- to fully distributed models, and were compared for their performance efficiency in simulating groundwater recharge and karst spring hydrographs. The impact of the spatial distribution of recharge (land use, vegetation, precipitation) on the discharge was found to be low, as the semi- and fully distributed models had a comparable performance to the lumped reservoir models. The modelers stated, however, that the relative significance of the spatial distribution of the recharge function depends on the watershed characteristics. Other studies showed the substantial impact of vegetation and soil parameters (e.g., leaf area index, root depth, soil hydraulic conductivity and moisture content at saturation) on the evapotranspiration/infiltration and recharge functions in karst watersheds [29]. Sarrazin et al. [30] have found that karst recharge is not only sensitive to climatic factors but also to changes in land cover, using the large scale semi-distributed karst recharge model V2Karst V1.1. Bittner et al. [31] successfully simulated spring discharge in a karst watershed using the semi-distributed model Land use change modeling in KARSt systems (LuKARS). They validated the impact of land-use change on the spring water supply. Different techniques have also been developed to estimate evapotranspiration from remote sensing data and estimate groundwater recharge based on the reconstruction of distributed models of precipitation and evapotranspiration. These approaches were classified into four categories, namely: (1) the empirical direct methods, (2) the residual methods of the energy budget, (3) the deterministic methods, and (4) the vegetation index methods [32]. Ollivier et al. [33] established that incorporating the remote sensing-driven evapotranspiration model Simple Crop coefficient for Evapotranspiration (SimpKcET) into the grid-based distributed Karst Recharge and discharge Model (KaRaMel) improved the karst spring discharge simulations in low, intermediate, and high flow seasons, as well as the correlation between recharge events and recharge volumes. Yang et al. [34] also concluded that accounting for the heterogeneous spatial distribution of land cover and karst geological properties in a conceptually-based distributed karst hydrological model, referred to as the distributed karst Xinanjiang (DK-XAJ) model, improved the runoff simulation and the separation of groundwater recharge into rapid conduit flow and slow matrix flow.

As the need to predict future karst water resources under climate change projections and scenarios of land-use change increased, the use of the Soil and Water Assessment Tool (SWAT) in karst hydrological studies has gradually gained popularity. SWAT is a time-continuous, semi-distributed, process-based model that is widely used to simulate the spatial and temporal evolution of a catchment hydrological cycle, soil erosion and water quality, as well as the effects of land-use change and climate variability on catchment processes by dividing the watershed into subbasins and further into HRUs based on the land use, soil characteristics, and slope [35]. Numerous studies have directly applied the standard SWAT in karst basins, while others have modified its source code to improve the
(concentrated fast flows and diffuse slow flows are not explicitly simulated), (2) the double continuum approach, which represents the matrix and karst conduits as two interacting continua with their hydraulic attributes, and (3) the combined discrete-continuum approach in which the karst conduits are embedded as discrete elements inside the matrix [26,27]. The use of distributed parameters models remains a challenge due to the complexity of karst aquifer mechanisms and the need for extensive hydrological and hydrogeological investigations to define the characteristics of the karst system (i.e., aquifer geometry, conduits network geometry and location, hydraulic properties and interactions between the matrix and conduits) [23,26].

## 1.3. *Rationale of This Review of Karst Hydrological Modeling with SWAT*

Improving the management and sustainability of karst groundwater resources remains a challenge. This is due to the limited understanding of the critical zone processes in karst watersheds across space and time, as well as the lack of research that characterizes the influence of vegetation cover, climate change, and anthropogenic activities on these processes [8,9]. Jeannin et al. [28] recently tested 13 karst numerical models in a karst watershed. The models consisted of lumped neural networks, reservoir models, and semi- to fully distributed models, and were compared for their performance efficiency in simulating groundwater recharge and karst spring hydrographs. The impact of the spatial distribution of recharge (land use, vegetation, precipitation) on the discharge was found to be low, as the semi- and fully distributed models had a comparable performance to the lumped reservoir models. The modelers stated, however, that the relative significance of the spatial distribution of the recharge function depends on the watershed characteristics. Other studies showed the substantial impact of vegetation and soil parameters (e.g., leaf area index, root depth, soil hydraulic conductivity and moisture content at saturation) on the evapotranspiration/infiltration and recharge functions in karst watersheds [29]. Sarrazin et al. [30] have found that karst recharge is not only sensitive to climatic factors but also to changes in land cover, using the large scale semi-distributed karst recharge model V2Karst V1.1. Bittner et al. [31] successfully simulated spring discharge in a karst watershed using the semi-distributed model Land use change modeling in KARSt systems (LuKARS). They validated the impact of land-use change on the spring water supply. Different techniques have also been developed to estimate evapotranspiration from remote sensing data and estimate groundwater recharge based on the reconstruction of distributed models of precipitation and evapotranspiration. These approaches were classified into four categories, namely: (1) the empirical direct methods, (2) the residual methods of the energy budget, (3) the deterministic methods, and (4) the vegetation index methods [32]. Ollivier et al. [33] established that incorporating the remote sensing-driven evapotranspiration model Simple Crop coefficient for Evapotranspiration (SimpKcET) into the grid-based distributed Karst Recharge and discharge Model (KaRaMel) improved the karst spring discharge simulations in low, intermediate, and high flow seasons, as well as the correlation between recharge events and recharge volumes. Yang et al. [34] also concluded that accounting for the heterogeneous spatial distribution of land cover and karst geological properties in a conceptually-based distributed karst hydrological model, referred to as the distributed karst Xinanjiang (DK-XAJ) model, improved the runoff simulation and the separation of groundwater recharge into rapid conduit flow and slow matrix flow.

As the need to predict future karst water resources under climate change projections and scenarios of land-use change increased, the use of the Soil and Water Assessment Tool (SWAT) in karst hydrological studies has gradually gained popularity. SWAT is a time-continuous, semi-distributed, process-based model that is widely used to simulate the spatial and temporal evolution of a catchment hydrological cycle, soil erosion and water quality, as well as the effects of land-use change and climate variability on catchment processes by dividing the watershed into subbasins and further into HRUs based on the land use, soil characteristics, and slope [35]. Numerous studies have directly applied the standard SWAT in karst basins, while others have modified its source code to improve the
representation of karst hydrology by considering different karst flow regimes and features (e.g., sinkholes, springs, and IGF) [36]. Therefore, there exists a need for a comprehensive review on the applications of SWAT in karst watersheds in order to: (a) identify the areas of research in the water resource field in which SWAT was implemented in karstified watersheds, (b) investigate the different approaches in which karst processes were incorporated into SWAT, and (c) evaluate the modified SWAT models performance and adequacy in representing the heterogeneous and non-linear karst flow mechanisms.

## 1.4. Overview of Previous SWAT Applications and Modifications for Karst Modeling

The earliest SWAT studies in karst watersheds (a total of 4 articles) have been reported by Gassman et al. [37] as part of a full range review of research findings and methods for different application categories with SWAT (e.g., discharge, hydrological analyses, sensitivity analyses and calibration techniques, climate change impacts on hydrology, pollutant transport and fate). Their review was based on more than 250 SWAT articles identified in the literature up to the year 2007. Since then, the use of SWAT has seen a tremendous growth globally for a wide range of scales and complex environmental studies, with more than 5000 articles currently published in peer-reviewed journals [38]. The number of SWAT review studies has also expanded to cover a variety of applications, such as: SWAT developments in landscape representation, stream routing, and soil phosphorus dynamics [39], SWAT improvements in addressing environmental issues [40], quantification of ecosystem services [41], runoff simulation, hydrological impacts under changing environment, and non-point source pollution [42], SWAT limitations in simulating subdaily processes [43], methods used to develop a SWAT model at field-scale [44], SWAT simulations of hydroclimatic extremes [45], and SWAT applications in Mediterranean catchments [46], to name a few.

Despite these advancements, the numerical simulation of karst watersheds and their processes in SWAT is still underway. In fact, a recent research study by Eini et al. [36] cited only 36 articles describing SWAT-based applications in partially karstified and karst dominated watersheds, with just 13 studies featuring a modified SWAT code. To note, Eini et al. [36] did not provide a detailed overview of the karst modeling approaches adopted in articles that they cited but rather an introductory synopsis prior to presenting two modified SWAT codes that they developed and applied in a karst watershed. Therefore, our paper is the first—to our best knowledge—to present an in-depth review of the studies conducted with SWAT in karst watersheds, building on the selected list of publications by Eini et al. [36] and extending to the full range of studies between the years 2000–2022. The objectives of our present review are to: (a) describe the SWAT subroutines that correspond to the different processes driving the flow of water in the critical zone (i.e., surface runoff, evapotranspiration, infiltration, interflow, recharge, baseflow), (b) summarize and discuss the research methods and findings for the standard and modified SWAT models in karst influenced watersheds, (c) identify potential constraints of the existing SWAT modeling approaches in representing the heterogeneous and non-linear flow mechanisms in karst aquifers, and (d) propose future research directions in order to enhance the applicability of SWAT in karst watersheds and the reliability assessment of karst water resources for future management and planning.

This review will present the different applications of SWAT (i.e., water quantity and quality, land-use and climate change, erosion processes, ecohydrological assessment, and water resources management) in karst influenced and karst dominated watersheds. However, the primary focus of the discussion will be the hydrological assessment in the SWAT applications that featured SWAT coupling with other hydrological models or modifications to the SWAT recharge and groundwater flow equations. These studies aimed to improve the representation of karst features, baseflow, and peak flows in SWAT prior to simulating other watershed processes, such as sediment or pollutant transport.
representation of karst hydrology by considering different karst flow regimes and features (e.g., sinkholes, springs, and IGF) [36]. Therefore, there exists a need for a comprehensive review on the applications of SWAT in karst watersheds in order to: (a) identify the areas of research in the water resource field in which SWAT was implemented in karstified watersheds, (b) investigate the different approaches in which karst processes were incorporated into SWAT, and (c) evaluate the modified SWAT models performance and adequacy in representing the heterogeneous and non-linear karst flow mechanisms.

## 1.4. Overview of Previous SWAT Applications and Modifications for Karst Modeling

The earliest SWAT studies in karst watersheds (a total of 4 articles) have been reported by Gassman et al. [37] as part of a full range review of research findings and methods for different application categories with SWAT (e.g., discharge, hydrological analyses, sensitivity analyses and calibration techniques, climate change impacts on hydrology, pollutant transport and fate). Their review was based on more than 250 SWAT articles identified in the literature up to the year 2007. Since then, the use of SWAT has seen a tremendous growth globally for a wide range of scales and complex environmental studies, with more than 5000 articles currently published in peer-reviewed journals [38]. The number of SWAT review studies has also expanded to cover a variety of applications, such as: SWAT developments in landscape representation, stream routing, and soil phosphorus dynamics [39], SWAT improvements in addressing environmental issues [40], quantification of ecosystem services [41], runoff simulation, hydrological impacts under changing environment, and non-point source pollution [42], SWAT limitations in simulating subdaily processes [43], methods used to develop a SWAT model at field-scale [44], SWAT simulations of hydroclimatic extremes [45], and SWAT applications in Mediterranean catchments [46], to name a few.

Despite these advancements, the numerical simulation of karst watersheds and their processes in SWAT is still underway. In fact, a recent research study by Eini et al. [36] cited only 36 articles describing SWAT-based applications in partially karstified and karst dominated watersheds, with just 13 studies featuring a modified SWAT code. To note, Eini et al. [36] did not provide a detailed overview of the karst modeling approaches adopted in articles that they cited but rather an introductory synopsis prior to presenting two modified SWAT codes that they developed and applied in a karst watershed. Therefore, our paper is the first—to our best knowledge—to present an in-depth review of the studies conducted with SWAT in karst watersheds, building on the selected list of publications by Eini et al. [36] and extending to the full range of studies between the years 2000–2022. The objectives of our present review are to: (a) describe the SWAT subroutines that correspond to the different processes driving the flow of water in the critical zone (i.e., surface runoff, evapotranspiration, infiltration, interflow, recharge, baseflow), (b) summarize and discuss the research methods and findings for the standard and modified SWAT models in karst influenced watersheds, (c) identify potential constraints of the existing SWAT modeling approaches in representing the heterogeneous and non-linear flow mechanisms in karst aquifers, and (d) propose future research directions in order to enhance the applicability of SWAT in karst watersheds and the reliability assessment of karst water resources for future management and planning.

This review will present the different applications of SWAT (i.e., water quantity and quality, land-use and climate change, erosion processes, ecohydrological assessment, and water resources management) in karst influenced and karst dominated watersheds. However, the primary focus of the discussion will be the hydrological assessment in the SWAT applications that featured SWAT coupling with other hydrological models or modifications to the SWAT recharge and groundwater flow equations. These studies aimed to improve the representation of karst features, baseflow, and peak flows in SWAT prior to simulating other watershed processes, such as sediment or pollutant transport.
## 2. Equations in SWAT for Hydrological Simulation

SWAT is a continuous-time semi-distributed agro-eco-hydrological model developed by the United States Department of Agriculture (USDA) [47]. The model has been successfully applied to monitor and predict the impacts of environmental and anthropogenic changes on the physical processes in watersheds at small, regional, and subcontinental scales (e.g., [48–57]). SWAT uses meteorological data, i.e., precipitation, air temperature, relative humidity, wind speed, and solar radiation, in addition to topography, soil properties, and land-use data, to simulate the watershed water balance components at different time steps (subdaily to annual). It can also model water quality and soil erosion [43,58]. The watershed is first disaggregated into subbasins connected through a stream channel and further into HRUs that represent areas of homogenous land use, soil, and slope properties [59]. The definition of HRUs is performed using a geographic information system (GIS), such as the ArcSWAT interface of ArcGIS or the QSWAT plugin of QGIS, coupled to the SWAT model to integrate the topographic, soil, and land-use inputs [60]. The simulated catchment processes in SWAT include surface runoff, infiltration, evapotranspiration, lateral flow, tile drainage, percolation, water stored in the soil profile, return flow from unconfined aquifers, consumptive water use through pumping (if any), recharge from surface water bodies, and in-stream processes, such as channel routing (main and tributary) and transformation of nutrients and pesticides [61]. These components are represented in each HRU by five storage volumes, namely the canopy interception, snow pack, soil profile, shallow aquifer, and deep aquifer [62].

Watershed hydrology in SWAT is represented by a land phase and a routing phase, whereby runoff, sediments, and agricultural chemical yields from all subbasin HRUs are aggregated to the main reach of the subbasin and routed through the channel network to the outlet(s) of the main catchment [63]. The fundamental daily water balance equation used in SWAT to represent the land phase of the hydrological cycle is given as follows [64]:

$$SW_t - SW_0 = \sum_{i=1}^{t} (P_{day} - Q_{surf} - ET_a - W_{seep} - Q_{gw}) \quad (1)$$

where $SW_0$ and $SW_t$ are the initial and final soil water content of the entire soil profile for the simulation period, respectively, $P_{day}$, $Q_{surf}$, $ET_a$, $W_{seep}$, and $Q_{gw}$ are precipitation, surface runoff, actual evapotranspiration, percolation and bypass flow exiting the soil bottom to the vadose zone, and return flow, respectively (all variables are expressed in mm H$_2$O.day$^{-1}$).

SWAT offers different options to simulate scheduled irrigation and auto-irrigation of crops. The auto-irrigation approach is generally used when irrigation scheduling data are lacking. Auto-irrigation is triggered by two stress identifiers: (1) plant water stress, whereby irrigation is applied to meet the plant water demand if the ratio of actual transpiration to potential transpiration falls below a user-specified threshold, and (2) soil water deficit, whereby irrigation is applied if the water content in the soil profile drops below field capacity by more than a user-defined soil water depletion threshold [65]. Sources of irrigation include river reaches, reservoirs, shallow and deep aquifers, or a source from outside the watershed, and irrigation demand is met based on the source water availability [66]. When irrigation is applied, the SWAT water balance is adjusted as follows:

$$SW_t - SW_0 = \sum_{i=1}^{t} (P + Irr) - \sum_{i=1}^{t} (Q_{surf} - ET_a - W_{seep} - Q_{gw}) \quad (2)$$

where $SW_0$ and $SW_t$ are the initial and the final soil water content of the entire soil profile for the simulation period, $P$, $Irr$, $Q_{surf}$, $ET_a$, $W_{seep}$, and $Q_{gw}$ represent precipitation, irrigation, surface runoff, actual evapotranspiration, percolation and bypass flow exiting the soil bottom to the vadose zone, and return flow, respectively (all variables are expressed in mm H$_2$O.day$^{-1}$).

Water routed through channels to the main watershed outlets is generated from direct surface runoff, lateral soil flow, baseflow from groundwater storage, and tile flow [67].
## 2. Equations in SWAT for Hydrological Simulation

SWAT is a continuous-time semi-distributed agro-eco-hydrological model developed by the United States Department of Agriculture (USDA) [47]. The model has been successfully applied to monitor and predict the impacts of environmental and anthropogenic changes on the physical processes in watersheds at small, regional, and subcontinental scales (e.g., [48–57]). SWAT uses meteorological data, i.e., precipitation, air temperature, relative humidity, wind speed, and solar radiation, in addition to topography, soil properties, and land-use data, to simulate the watershed water balance components at different time steps (subdaily to annual). It can also model water quality and soil erosion [43,58]. The watershed is first disaggregated into subbasins connected through a stream channel and further into HRUs that represent areas of homogenous land use, soil, and slope properties [59]. The definition of HRUs is performed using a geographic information system (GIS), such as the ArcSWAT interface of ArcGIS or the QSWAT plugin of QGIS, coupled to the SWAT model to integrate the topographic, soil, and land-use inputs [60]. The simulated catchment processes in SWAT include surface runoff, infiltration, evapotranspiration, lateral flow, tile drainage, percolation, water stored in the soil profile, return flow from unconfined aquifers, consumptive water use through pumping (if any), recharge from surface water bodies, and in-stream processes, such as channel routing (main and tributary) and transformation of nutrients and pesticides [61]. These components are represented in each HRU by five storage volumes, namely the canopy interception, snow pack, soil profile, shallow aquifer, and deep aquifer [62].

Watershed hydrology in SWAT is represented by a land phase and a routing phase, whereby runoff, sediments, and agricultural chemical yields from all subbasin HRUs are aggregated to the main reach of the subbasin and routed through the channel network to the outlet(s) of the main catchment [63]. The fundamental daily water balance equation used in SWAT to represent the land phase of the hydrological cycle is given as follows [64]:

$$SW_t - SW_0 = \sum_{i=1}^{t} (P_{day} - Q_{surf} - ET_a - W_{seep} - Q_{gw}) \quad (1)$$

where $SW_0$ and $SW_t$ are the initial and final soil water content of the entire soil profile for the simulation period, respectively, $P_{day}$, $Q_{surf}$, $ET_a$, $W_{seep}$, and $Q_{gw}$ are precipitation, surface runoff, actual evapotranspiration, percolation and bypass flow exiting the soil bottom to the vadose zone, and return flow, respectively (all variables are expressed in mm H$_2$O.day$^{-1}$).

SWAT offers different options to simulate scheduled irrigation and auto-irrigation of crops. The auto-irrigation approach is generally used when irrigation scheduling data are lacking. Auto-irrigation is triggered by two stress identifiers: (1) plant water stress, whereby irrigation is applied to meet the plant water demand if the ratio of actual transpiration to potential transpiration falls below a user-specified threshold, and (2) soil water deficit, whereby irrigation is applied if the water content in the soil profile drops below field capacity by more than a user-defined soil water depletion threshold [65]. Sources of irrigation include river reaches, reservoirs, shallow and deep aquifers, or a source from outside the watershed, and irrigation demand is met based on the source water availability [66]. When irrigation is applied, the SWAT water balance is adjusted as follows:

$$SW_t - SW_0 = \sum_{i=1}^{t} (P + Irr) - \sum_{i=1}^{t} (Q_{surf} - ET_a - W_{seep} - Q_{gw}) \quad (2)$$

where $SW_0$ and $SW_t$ are the initial and the final soil water content of the entire soil profile for the simulation period, $P$, $Irr$, $Q_{surf}$, $ET_a$, $W_{seep}$, and $Q_{gw}$ represent precipitation, irrigation, surface runoff, actual evapotranspiration, percolation and bypass flow exiting the soil bottom to the vadose zone, and return flow, respectively (all variables are expressed in mm H$_2$O.day$^{-1}$).

Water routed through channels to the main watershed outlets is generated from direct surface runoff, lateral soil flow, baseflow from groundwater storage, and tile flow [67].
These flow components contribute to the catchment water yield, which is considered a critical parameter in the sustainable management of water resources [68]. Hence, water yield is defined as the net water volume leaving the HRU and entering a reach at the subbasin level into the main channel, as follows:

$$WYLD = Q_{surf} + Q_{lat} + Q_{gw} + Q_{tile} - T_{loss} \quad (3)$$

where WYLD is the water yield, $Q_{surf}$, $Q_{lat}$, and $Q_{gw}$ are the surface runoff, soil lateral flow, and return flow from the shallow aquifer to the main channel, respectively, $Q_{tile}$ is the tile flow, and $T_{loss}$ represents the losses from the tributary in the HRU via transmission through the riverbed (all variables are expressed in mm H₂O.day⁻¹).

The governing equations of the watershed hydrological components are presented thoroughly in the theoretical documentation of SWAT [69]. Hence, the primary focus in Sections 2.1 and 2.2 below was devoted to flow processes in the critical zone that directly impact streamflow simulation in standard SWAT. These processes were grouped under: (1) surface water hydrology (Section 2.1) and (2) subsurface water hydrology (Section 2.2). A set of equations that are fundamental in hydrological modeling with SWAT were provided in each subsection, with a corresponding list of SWAT variables in Appendix A (Table A1). These equations serve to help the reader understand the methods used in SWAT to simulate surface and groundwater flows, prior to discussing the modifications made to the SWAT source code for the different applications in the realm of karst hydrology (Section 4 of the review).

## 2.1. Surface Water Hydrology

### 2.1.1. Evapotranspiration

SWAT provides three methods to simulate daily potential evapotranspiration (PET) at the HRU scale, namely the Penman-Monteith [70], the Priestley-Taylor [71], and the Hargreaves methods [72]. Between the three approaches, the Penman-Monteith equation is considered the most suited to estimate PET, as it explicitly separates the effects of climate and land cover properties on each of the evapotranspiration components [30,33]. This method is represented with Equation (4), as follows [69]:

$$\lambda E = \frac{\Delta(H_{net} - G) + \rho_{air} \times c_p (e_z^0 - e_z) / r_a}{\Delta + \gamma(1 + r_c / r_a)} \quad (4)$$

where $\lambda E$ is the latent heat flux density (MJ.m⁻².day⁻¹), $\lambda$ is the latent heat of vaporization (MJ.kg⁻¹), $E$ is the depth rate evaporation (mm.day⁻¹), $\Delta$ is the slope of the saturation vapor pressure-temperature curve (de/dT) (kPa.°C⁻¹), $H_{net}$ is the net radiation (MJ.m⁻².day⁻¹), $G$ is the heat flux density to the ground (MJ.m⁻².day⁻¹), $\rho_{air}$ is the air density (kg.m⁻³), $c_p$ is the specific heat at constant pressure (MJ.kg⁻¹.C⁻¹), $e_z^0$ is the saturation vapor pressure of air at height z (kPa), $e_z$ is the water vapor pressure of air at height z (kPa), $\gamma$ is the psychometric constant (kPa.°C⁻¹), $r_c$ is the plant canopy resistance (s.m⁻¹), and $r_a$ is the diffusion resistance of the air layer (aerodynamic resistance) (s.m⁻¹).

PET in SWAT depends on plant growth, which considers canopy resistance expressed as a function of the minimum effective stomatal resistance for a single leaf and the leaf area index (LAI). The LAI, defined as one half the total leaf area per unit ground area, reflects the structural characteristics of the plant canopy and defines the size of the interface for energy and mass exchanges between the vegetation surface and the atmosphere [73]. Evapotranspiration is also related to the canopy height required to determine the aerodynamic resistance parameter [69].

SWAT uses the LAI in conjunction with a simplified version of the Environmental Policy Integrated Climate (EPIC) plant growth model to simulate the phenological development of plants and estimate evapotranspiration [74]. In addition to the LAI development, the plant growth module of SWAT includes the simulation of the light interception and the conversion of intercepted light into biomass, assuming a plant species-specific radiation-
These flow components contribute to the catchment water yield, which is considered a critical parameter in the sustainable management of water resources [68]. Hence, water yield is defined as the net water volume leaving the HRU and entering a reach at the subbasin level into the main channel, as follows:

$$WYLD = Q_{surf} + Q_{lat} + Q_{gw} + Q_{tile} - T_{loss} \quad (3)$$

where WYLD is the water yield, $Q_{surf}$, $Q_{lat}$, and $Q_{gw}$ are the surface runoff, soil lateral flow, and return flow from the shallow aquifer to the main channel, respectively, $Q_{tile}$ is the tile flow, and $T_{loss}$ represents the losses from the tributary in the HRU via transmission through the riverbed (all variables are expressed in mm H$_{2}$O.day$^{-1}$).

The governing equations of the watershed hydrological components are presented thoroughly in the theoretical documentation of SWAT [69]. Hence, the primary focus in Sections 2.1 and 2.2 below was devoted to flow processes in the critical zone that directly impact streamflow simulation in standard SWAT. These processes were grouped under: (1) surface water hydrology (Section 2.1) and (2) subsurface water hydrology (Section 2.2). A set of equations that are fundamental in hydrological modeling with SWAT were provided in each subsection, with a corresponding list of SWAT variables in Appendix A (Table A1). These equations serve to help the reader understand the methods used in SWAT to simulate surface and groundwater flows, prior to discussing the modifications made to the SWAT source code for the different applications in the realm of karst hydrology (Section 4 of the review).

## 2.1. Surface Water Hydrology

### 2.1.1. Evapotranspiration

SWAT provides three methods to simulate daily potential evapotranspiration (PET) at the HRU scale, namely the Penman-Monteith [70], the Priestley-Taylor [71], and the Hargreaves methods [72]. Between the three approaches, the Penman-Monteith equation is considered the most suited to estimate PET, as it explicitly separates the effects of climate and land cover properties on each of the evapotranspiration components [30,33]. This method is represented with Equation (4), as follows [69]:

$$\lambda E = \frac{\Delta(H_{net} - G) + \rho_{air} \times c_p (e_z^0 - e_z) / r_a}{\Delta + \gamma(1 + r_c / r_a)} \quad (4)$$

where $\lambda E$ is the latent heat flux density (MJ.m$^{-2}$.day$^{-1}$), $\lambda$ is the latent heat of vaporization (MJ.kg$^{-1}$), $E$ is the depth rate evaporation (mm.day$^{-1}$), $\Delta$ is the slope of the saturation vapor pressure-temperature curve (de/dT) (kPa.$^{\circ}$C$^{-1}$), $H_{net}$ is the net radiation (MJ.m$^{-2}$.day$^{-1}$), $G$ is the heat flux density to the ground (MJ.m$^{-2}$.day$^{-1}$), $\rho_{air}$ is the air density (kg.m$^{-3}$), $c_p$ is the specific heat at constant pressure (MJ.kg$^{-1}$.C$^{-1}$), $e_z^0$ is the saturation vapor pressure of air at height $z$ (kPa), $e_z$ is the water vapor pressure of air at height $z$ (kPa), $\gamma$ is the psychometric constant (kPa.$^{\circ}$C$^{-1}$), $r_c$ is the plant canopy resistance (s.m$^{-1}$), and $r_a$ is the diffusion resistance of the air layer (aerodynamic resistance) (s.m$^{-1}$).

PET in SWAT depends on plant growth, which considers canopy resistance expressed as a function of the minimum effective stomatal resistance for a single leaf and the leaf area index (LAI). The LAI, defined as one half the total leaf area per unit ground area, reflects the structural characteristics of the plant canopy and defines the size of the interface for energy and mass exchanges between the vegetation surface and the atmosphere [73]. Evapotranspiration is also related to the canopy height required to determine the aerodynamic resistance parameter [69].

SWAT uses the LAI in conjunction with a simplified version of the Environmental Policy Integrated Climate (EPIC) plant growth model to simulate the phenological development of plants and estimate evapotranspiration [74]. In addition to the LAI development, the plant growth module of SWAT includes the simulation of the light interception and the conversion of intercepted light into biomass, assuming a plant species-specific radiation-
use efficiency [75]. Plant development is primarily dependent on the base temperature for growth derived from minimum, maximum, and optimum temperature requirements. The plants heat unit requirements are quantified and related to the time of planting and maturity [76]. The LAI is incremented daily based on the accumulated potential heat units. It first increases to a crop-specific maximum value, remains constant until the senescence stage, then decreases linearly to zero at harvest. Similarly, the canopy height increases until a crop-specific maximum is achieved and stays at this height through the remainder of the growing season [77]. The potential crop leaf growth and biomass are first computed under optimal conditions and further adjusted for actual growth under stress factors such as water, temperature, and nutrients [78]. SWAT also uses dormancy in function of day length and latitude to repeat the annual growth cycle for trees and perennials [73].

After estimating potential evapotranspiration, SWAT calculates actual evapotranspiration (ETa), which includes four components: the canopy evaporation, the plant transpiration, the sublimation and soil surface evaporation, and the groundwater evapotranspiration [79]. The model first evaporates any precipitation intercepted by the plant canopy. Then, actual plant transpiration is estimated as a function of the potential transpiration adjusted for the wet canopy storage, root depth, soil water content, and the leaf area index, which depends on the plant developmental stage. Soil evaporation is modeled as a function of potential evapotranspiration adjusted for canopy evaporation and the rate of shading. If snow is present in the HRUs, sublimation takes place until evaporation from soil could occur after snow melting. Subsequently, SWAT proceeds to adjust the maximum possible soil evaporation for plant water use and partitions the evaporative demand between the different soil layers, in order to estimate the actual evaporation at each layer based on the soil water content [74,79,80].

### 2.1.2. Surface Runoff and Infiltration

In SWAT, soil surface runoff and infiltration are estimated from precipitation by one of the following two approaches: (1) the modified Soil Conservation Service Curve Number (SCS-CN) procedure and (2) the Green and Ampt Mein Larson (GAML) excess rainfall method. The SCS-CN approach simulates cumulative surface runoff based on cumulative precipitation and soil retention properties for daily time step, whereas the GAML approach simulates surface runoff for subdaily time step applications using subdaily precipitation input data [45,81,82].

Surface runoff is estimated with the SCN-CN procedure as follows [66]:

$$Q_{surf} = \frac{(P - I_a)^2}{(P - I_a) + S} \quad \text{when } P > I_a; \text{ else } Q_{surf} = 0 \quad (5)$$

where $Q_{surf}$ is the accumulated runoff, $P$ is the total precipitation, $I_a$ is the initial water abstraction prior to runoff due to surface storage interception and infiltration (generally approximated as 0.2S, but can vary with the soil type), and $S$ is the soil moisture retention parameter (all variables are expressed in mm H$_2$O.day$^{-1}$).

The soil retention parameter (S) varies temporally with the changes in moisture content, and spatially in function of the soil type, land use, and management practices. It can also be assumed to vary with the accumulated plant evapotranspiration. The retention parameter is expressed as a function of the daily curve number (CN) of the Antecedent Moisture Condition-II (AMC-II) for a given land use/cover and hydrological soil group as follows [66]:

$$S = \frac{25,400}{CN} - 254 \quad (6)$$

The SCS approach defines three antecedent moisture conditions, namely AMC-I for dry/wilting point condition, AMC-II for average moisture, and AMC-III wet/field capacity,
use efficiency [75]. Plant development is primarily dependent on the base temperature for growth derived from minimum, maximum, and optimum temperature requirements. The plants heat unit requirements are quantified and related to the time of planting and maturity [76]. The LAI is incremented daily based on the accumulated potential heat units. It first increases to a crop-specific maximum value, remains constant until the senescence stage, then decreases linearly to zero at harvest. Similarly, the canopy height increases until a crop-specific maximum is achieved and stays at this height through the remainder of the growing season [77]. The potential crop leaf growth and biomass are first computed under optimal conditions and further adjusted for actual growth under stress factors such as water, temperature, and nutrients [78]. SWAT also uses dormancy in function of day length and latitude to repeat the annual growth cycle for trees and perennials [73].

After estimating potential evapotranspiration, SWAT calculates actual evapotranspiration (ETa), which includes four components: the canopy evaporation, the plant transpiration, the sublimation and soil surface evaporation, and the groundwater evapotranspiration [79]. The model first evaporates any precipitation intercepted by the plant canopy. Then, actual plant transpiration is estimated as a function of the potential transpiration adjusted for the wet canopy storage, root depth, soil water content, and the leaf area index, which depends on the plant developmental stage. Soil evaporation is modeled as a function of potential evapotranspiration adjusted for canopy evaporation and the rate of shading. If snow is present in the HRUs, sublimation takes place until evaporation from soil could occur after snow melting. Subsequently, SWAT proceeds to adjust the maximum possible soil evaporation for plant water use and partitions the evaporative demand between the different soil layers, in order to estimate the actual evaporation at each layer based on the soil water content [74,79,80].

### 2.1.2. Surface Runoff and Infiltration

In SWAT, soil surface runoff and infiltration are estimated from precipitation by one of the following two approaches: (1) the modified Soil Conservation Service Curve Number (SCS-CN) procedure and (2) the Green and Ampt Mein Larson (GAML) excess rainfall method. The SCS-CN approach simulates cumulative surface runoff based on cumulative precipitation and soil retention properties for daily time step, whereas the GAML approach simulates surface runoff for subdaily time step applications using subdaily precipitation input data [45,81,82].

Surface runoff is estimated with the SCN-CN procedure as follows [66]:

$$Q_{surf} = \frac{(P - I_a)^2}{(P - I_a) + S} \quad \text{when } P > I_a; \text{ else } Q_{surf} = 0 \quad (5)$$

where $Q_{surf}$ is the accumulated runoff, $P$ is the total precipitation, $I_a$ is the initial water abstraction prior to runoff due to surface storage interception and infiltration (generally approximated as 0.2S, but can vary with the soil type), and $S$ is the soil moisture retention parameter (all variables are expressed in mm H$_2$O.day$^{-1}$).

The soil retention parameter (S) varies temporally with the changes in moisture content, and spatially in function of the soil type, land use, and management practices. It can also be assumed to vary with the accumulated plant evapotranspiration. The retention parameter is expressed as a function of the daily curve number (CN) of the Antecedent Moisture Condition-II (AMC-II) for a given land use/cover and hydrological soil group as follows [66]:

$$S = \frac{25,400}{CN} - 254 \quad (6)$$

The SCS approach defines three antecedent moisture conditions, namely AMC-I for dry/wilting point condition, AMC-II for average moisture, and AMC-III wet/field capacity,
represented by curve numbers CN1, CN2 and CN3, respectively. CN1 and CN3 are computed as a function of CN2 as follows:

$$CN1 = CN2 - \frac{20(100 - CN2)}{(100 - CN2 + e^{2.533-0.0636(100-CN2)})} \quad (7)$$

$$CN3 = CN2 \times e^{0.000673(100-CN2)} \quad (8)$$

Infiltration rate is calculated using the GAML equation as follows [83]:

$$I(t) = K_e \left( 1 + \frac{\psi \times \Delta\theta}{I_{acc}(t)} \right) \quad (9)$$

where $I(t)$ is the infiltration rate (mm H$_2$O) at the simulation time step $t$ (subdaily), $K_e$ is the effective hydraulic conductivity, which considers soil water content and land-use impact as a function of CN (mm.h$^{-1}$), $\psi$ is the wetting front matric potential (mm), $\Delta\theta$ is the change in soil moisture content (mm.mm$^{-1}$), and $I_{acc}(t)$ is the cumulative infiltration after ponding (mm H$_2$O hour$^{-1}$). The cumulative depth of water infiltration $I_{acc}(t)$ is computed as follows:

$$I_{acc}(t) = I_{acc}(t-1) + K_e \times \Delta t - \psi \times \Delta\theta \times \ln \left[ \frac{I_{acc}(t) + \psi \times \Delta\theta}{I_{acc}(t-1) + \psi \times \Delta\theta} \right] \quad (10)$$

where $t-1$ is the previous simulation time step. Equation (10) is solved using a successive substitution technique. Subsequently, the infiltration rate is calculated using Equation (9) for each time step. Surface runoff is generated when the rainfall intensity exceeds infiltration rate. Otherwise, the total rainfall volume during the time step infiltrates into the soil.

### 2.1.3. Channel Flow and Flow Routing

For stream channel routing, Manning's equation is used to calculate the rate and velocity of flow in the reach of each subbasin when the streamflow is less than the bankfull discharge rate, computed as a function of the bankfull channel width and depth. SWAT incorporates floodplain inundation geometry into the channel routing simulation if the streamflow is greater than bankfull flow [84].

The peak runoff rate, reached when all the subbasins are contributing to flow at the outlet, is estimated using the modified rational method, as follows [85]:

$$q_{peak} = \frac{\alpha_{tc} \times Q_{surf} \times A}{3.6 \times t_{conc}} \quad (11)$$

where $q_{peak}$ is the peak runoff rate (m$^3$.s$^{-1}$), $\alpha_{tc}$ is the fraction of daily rainfall that occurs during the time of concentration, $Q_{surf}$ is the surface runoff (mm H$_2$O.day$^{-1}$), A is the sub-basin area (km$^2$), and $t_{conc}$ is the time of concentration for the subbasin (hours), calculated as the sum of the overland flow time and channel flow time. Water is routed through the channel network using either the Muskingum routing method (based on the continuity and empirical linear storage equations) [86] or the variable storage routing method (based on the continuity equation) [87,88].

Water transmission losses can occur through the side and bottom of the river channels and enter the bank storage or the deep aquifer. Transmission losses are estimated as follows [89]:

$$T_{loss} = K_{ch} \times L_{ch} \times P_{ch} \times TT \quad (12)$$

where $T_{loss}$ represents the channel transmission losses (m$^3$ H$_2$O), $K_{ch}$ is the effective hydraulic conductivity of the channel alluvium (mm.h$^{-1}$), $L_{ch}$ is the channel length (km), $P_{ch}$ is the wetted perimeter in the channel (m), and $TT$ is the flow travel time (hours).
represented by curve numbers CN1, CN2 and CN3, respectively. CN1 and CN3 are computed as a function of CN2 as follows:

$$CN1 = CN2 - \frac{20(100 - CN2)}{(100 - CN2 + e^{2.533-0.0636(100-CN2)})} \quad (7)$$

$$CN3 = CN2 \times e^{0.000673(100-CN2)} \quad (8)$$

Infiltration rate is calculated using the GAML equation as follows [83]:

$$I(t) = K_e \left( 1 + \frac{\psi \times \Delta\theta}{I_{acc}(t)} \right) \quad (9)$$

where $I(t)$ is the infiltration rate (mm H₂O) at the simulation time step $t$ (subdaily), $K_e$ is the effective hydraulic conductivity, which considers soil water content and land-use impact as a function of CN (mm.h⁻¹), $\psi$ is the wetting front matric potential (mm), $\Delta\theta$ is the change in soil moisture content (mm.mm⁻¹), and $I_{acc}(t)$ is the cumulative infiltration after ponding (mm H₂O hour⁻¹). The cumulative depth of water infiltration $I_{acc}(t)$ is computed as follows:

$$I_{acc}(t) = I_{acc}(t-1) + K_e \times \Delta t - \psi \times \Delta\theta \times \ln \left[ \frac{I_{acc}(t) + \psi \times \Delta\theta}{I_{acc}(t-1) + \psi \times \Delta\theta} \right] \quad (10)$$

where $t-1$ is the previous simulation time step. Equation (10) is solved using a successive substitution technique. Subsequently, the infiltration rate is calculated using Equation (9) for each time step. Surface runoff is generated when the rainfall intensity exceeds infiltration rate. Otherwise, the total rainfall volume during the time step infiltrates into the soil.

### 2.1.3. Channel Flow and Flow Routing

For stream channel routing, Manning's equation is used to calculate the rate and velocity of flow in the reach of each subbasin when the streamflow is less than the bankfull discharge rate, computed as a function of the bankfull channel width and depth. SWAT incorporates floodplain inundation geometry into the channel routing simulation if the streamflow is greater than bankfull flow [84].

The peak runoff rate, reached when all the subbasins are contributing to flow at the outlet, is estimated using the modified rational method, as follows [85]:

$$q_{peak} = \frac{\alpha_{tc} \times Q_{surf} \times A}{3.6 \times t_{conc}} \quad (11)$$

where $q_{peak}$ is the peak runoff rate (m³.s⁻¹), $\alpha_{tc}$ is the fraction of daily rainfall that occurs during the time of concentration, $Q_{surf}$ is the surface runoff (mm H₂O.day⁻¹), A is the sub-basin area (km²), and $t_{conc}$ is the time of concentration for the subbasin (hours), calculated as the sum of the overland flow time and channel flow time. Water is routed through the channel network using either the Muskingum routing method (based on the continuity and empirical linear storage equations) [86] or the variable storage routing method (based on the continuity equation) [87,88].

Water transmission losses can occur through the side and bottom of the river channels and enter the bank storage or the deep aquifer. Transmission losses are estimated as follows [89]:

$$T_{loss} = K_{ch} \times L_{ch} \times P_{ch} \times TT \quad (12)$$

where $T_{loss}$ represents the channel transmission losses (m³ H₂O), $K_{ch}$ is the effective hydraulic conductivity of the channel alluvium (mm.h⁻¹), $L_{ch}$ is the channel length (km), $P_{ch}$ is the wetted perimeter in the channel (m), and $TT$ is the flow travel time (hours).
## 2.2. Subsurface Water Hydrology

### 2.2.1. Soil Water Percolation and Lateral Flow

The water percolation component in SWAT redistributes infiltrated water in the soil profile using a storage routing method combined with an optional crack-flow routine. Percolation is simulated when the water content of a soil layer exceeds its field capacity defined as the sum of the available soil water content and permanent wilting point. Percolated water moves to the subsequent layer unless it is saturated, frozen, or impervious [90–92]. Water percolation is estimated as follows:

$$W_{perc,ly} = SW_{ly,excess} \left( 1 - e^{\frac{-\Delta t}{TT_{perc,ly}}} \right) \quad (13)$$

where $W_{perc,ly}$ is the water percolating from soil layer (ly) to the underlying soil layer (mm H₂O.day⁻¹), $SW_{ly,excess}$ is the drainable volume of water in the soil layer (ly) on a given day (computed as the difference between the water content of the soil layer and field capacity, in mm H₂O.day⁻¹), $\Delta t$ is the length of the time step (hours), and $TT_{perc,ly}$ is the travel time through the soil layer (hours), calculated as follows:

$$TT_{perc,ly} = \frac{SAT_{ly} - FC_{ly}}{K_{sat,ly}} \quad (14)$$

where $K_{sat,ly}$ (mm.h⁻¹), $SAT_{ly}$ (mm H₂O), and $FC_{ly}$ (mm H₂O) represent the saturated hydraulic conductivity, saturation water content, and field capacity water content of the soil layer (ly), respectively.

SWAT incorporates a crack flow module that can be used to simulate bypass (crack) or preferential flow in the soil. The use of the crack flow approach to increase infiltration rates from the surface is optional and requires the activation of a crack flow code by the user [36]. Crack volume for each soil layer is modeled in the dry seasons, which allows infiltrated rainwater to move rapidly through the soil profile along vertical cracks, and disappears in wet conditions [93]. Bypass flow from the bottom of the soil profile to the saturated zone is computed using Equation (15), and excess water that leaves the bottom of the soil profile through the vadose zone is calculated by combining percolation and bypass flow, as shown in Equation (16) [69]:

$$W_{crk,btm} = 0.5 \times crk \left( \frac{crk_{ly=nn}}{depth_{ly=nn}} \right) \quad (15)$$

$$W_{seep} = W_{perc,ly=n} + W_{crk,btm} \quad (16)$$

where $W_{perc,ly=n}$ is the water percolating out of the lowest soil layer, $W_{crk,btm}$ is the crack flow past the lower boundary of the soil profile (mm H₂O.day⁻¹), $crk$ is the total crack volume for the soil profile on a given day (mm), $crk_{ly=nn}$ is the crack volume for the deepest soil layer on a given day (mm), $depth_{ly=nn}$ is the depth of the deepest soil layer (mm), and $W_{seep}$ is the total volume of water drained from the bottom of the soil profile (mm H₂O.day⁻¹).

Lateral flow (soil interflow) along a steep hillslope is computed simultaneously with percolation when the soil water content exceeds its field capacity. It is simulated using a kinematic storage routing method (Equation (17)) that is based on the on slope, slope length, and saturated conductivity of each soil layer [63,92], as follows:

$$Q_{lat} = 0.024 \left( \frac{2 \times SW_{ly,excess} \times K_{sat} \times slp}{\emptyset_d \times L_{hill}} \right) \quad (17)$$

where $Q_{lat}$ is the daily water flux from the hillslope outlet (mm H₂O.day⁻¹), $SW_{ly,excess}$ is the drainable volume of water stored in the saturated zone of the hillslope per unit area (mm H₂O), $K_{sat}$ is the saturated hydraulic conductivity of the soil (mm.h⁻¹), $slp$ is the
## 2.2. Subsurface Water Hydrology

### 2.2.1. Soil Water Percolation and Lateral Flow

The water percolation component in SWAT redistributes infiltrated water in the soil profile using a storage routing method combined with an optional crack-flow routine. Percolation is simulated when the water content of a soil layer exceeds its field capacity defined as the sum of the available soil water content and permanent wilting point. Percolated water moves to the subsequent layer unless it is saturated, frozen, or impervious [90–92]. Water percolation is estimated as follows:

$$W_{perc,ly} = SW_{ly,excess} \left( 1 - e^{\frac{-\Delta t}{TT_{perc,ly}}} \right) \qquad (13)$$

where $W_{perc,ly}$ is the water percolating from soil layer (ly) to the underlying soil layer (mm H₂O.day⁻¹), $SW_{ly,excess}$ is the drainable volume of water in the soil layer (ly) on a given day (computed as the difference between the water content of the soil layer and field capacity, in mm H₂O.day⁻¹), $\Delta t$ is the length of the time step (hours), and $TT_{perc,ly}$ is the travel time through the soil layer (hours), calculated as follows:

$$TT_{perc,ly} = \frac{SAT_{ly} - FC_{ly}}{K_{sat,ly}} \qquad (14)$$

where $K_{sat,ly}$ (mm.h⁻¹), $SAT_{ly}$ (mm H₂O), and $FC_{ly}$ (mm H₂O) represent the saturated hydraulic conductivity, saturation water content, and field capacity water content of the soil layer (ly), respectively.

SWAT incorporates a crack flow module that can be used to simulate bypass (crack) or preferential flow in the soil. The use of the crack flow approach to increase infiltration rates from the surface is optional and requires the activation of a crack flow code by the user [36]. Crack volume for each soil layer is modeled in the dry seasons, which allows infiltrated rainwater to move rapidly through the soil profile along vertical cracks, and disappears in wet conditions [93]. Bypass flow from the bottom of the soil profile to the saturated zone is computed using Equation (15), and excess water that leaves the bottom of the soil profile through the vadose zone is calculated by combining percolation and bypass flow, as shown in Equation (16) [69]:

$$W_{crk,btm} = 0.5 \times crk \left( \frac{crk_{ly=nn}}{depth_{ly=nn}} \right) \qquad (15)$$

$$W_{seep} = W_{perc,ly=n} + W_{crk,btm} \qquad (16)$$

where $W_{perc,ly=n}$ is the water percolating out of the lowest soil layer, $W_{crk,btm}$ is the crack flow past the lower boundary of the soil profile (mm H₂O.day⁻¹), $crk$ is the total crack volume for the soil profile on a given day (mm), $crk_{ly=nn}$ is the crack volume for the deepest soil layer on a given day (mm), $depth_{ly=nn}$ is the depth of the deepest soil layer (mm), and $W_{seep}$ is the total volume of water drained from the bottom of the soil profile (mm H₂O.day⁻¹).

Lateral flow (soil interflow) along a steep hillslope is computed simultaneously with percolation when the soil water content exceeds its field capacity. It is simulated using a kinematic storage routing method (Equation (17)) that is based on the on slope, slope length, and saturated conductivity of each soil layer [63,92], as follows:

$$Q_{lat} = 0.024 \left( \frac{2 \times SW_{ly,excess} \times K_{sat} \times slp}{\emptyset_d \times L_{hill}} \right) \qquad (17)$$

where $Q_{lat}$ is the daily water flux from the hillslope outlet (mm H₂O.day⁻¹), $SW_{ly,excess}$ is the drainable volume of water stored in the saturated zone of the hillslope per unit area (mm H₂O), $K_{sat}$ is the saturated hydraulic conductivity of the soil (mm.h⁻¹), $slp$ is the
increase in elevation per unit distance, $\emptyset_d$ is the drainable (residual) porosity of the soil layer (mm/mm), and $L_{hill}$ is the hillslope length (m).

The daily water balance for each soil layer is expressed using Equation (18), as follows [94]:

$$ \Delta SW_{ly} = W_{perc,ly-1} - W_{perc,ly} - Q_{lat,ly} - E_{e,ly} - E_{t,ly} \quad (18) $$

where $\Delta SW_{ly}$ is the change of soil water content at soil layer (ly), $W_{perc,ly-1}$ is the percolation received from layer (ly-1), $W_{perc,ly}$ and $Q_{lat,ly}$ are the percolation and lateral flow generated from soil layer (ly), respectively, and $E_{e,ly}$ and $E_{t,ly}$ are the evaporation and transpiration drawn from the soil layer (ly), respectively (all variables are expressed in mm H₂O.day⁻¹).

### 2.2.2. Groundwater Flow and Baseflow to the Stream

The groundwater module of SWAT comprises a system of two aquifers in each sub-basin: (1) a shallow unconfined aquifer that generates baseflow into the stream and (2) a deep confined aquifer contributing to streamflow outside of the watershed (flow lost from the system) [95]. Recharge from the unsaturated soil profile to the aquifers on a given day is calculated using an exponential decay weighting function that accounts for the time delay of the recharge mechanism, as follows [67]:

$$ W_{rchrg,i} = \left(1 - e^{-\frac{1}{\delta_{gw}}}\right) W_{seep,i} + \left(e^{-\frac{1}{\delta_{gw}}}\right) W_{rchrg,i-1} \quad (19) $$

where $W_{rchrg,i}$ and $W_{rchrg,i-1}$ represent the recharge to the aquifers (shallow and deep) at days i and i-1 (mm H₂O.day⁻¹), respectively, $W_{seep,i}$ is the water drained from the bottom of the soil profile (mm H₂O.day⁻¹), and $\delta_{gw}$ is the delay time required for recharge to reach the aquifers (days).

Recharge components routed to the shallow (unconfined) aquifer and the deep (confined) aquifer are computed using Equations (20) and (21), respectively, as follows:

$$ W_{rchrg,sh,i} = (1 - \beta_{dp}) W_{rchrg,i} \quad (20) $$

$$ W_{rchrg,dp,i} = \beta_{dp} W_{rchrg,i} \quad (21) $$

where $W_{rchrg,sh,i}$ and $W_{rchrg,dp,i}$ represent the water diverted to the shallow and deep aquifers (mm H₂O.day⁻¹), respectively, and $\beta_{dp}$ is a coefficient of percolation to the deep aquifer.

The shallow aquifer contributes to the streamflow if water stored in the aquifer exceeds a user-specified threshold. Otherwise, return flow is set to zero. The daily groundwater flow to the main river channel is computed using an exponential storage-discharge relationship, which incorporates the recharge from the shallow aquifer and a baseflow recession constant, as follows:

$$ Q_{gw,sh,i} = \begin{cases} W_{rchrg,sh,i} (1 - e^{-\alpha_{gw,sh} \times \Delta t}) + Q_{gw,sh,i-1} (e^{-\alpha_{gw,sh} \times \Delta t}), & aq_{sh} > aq_{shthr,q} \\ 0, & aq_{sh} \le aq_{shthr,q} \end{cases} \quad (22) $$

where $Q_{gw,sh,i}$ is the baseflow from the shallow aquifer to the main stream channel (mm H₂O.day⁻¹), $\alpha_{gw,sh}$ is the groundwater recession constant of shallow aquifer (days⁻¹), $aq_{sh}$ is the amount of water stored in the shallow aquifer (mm H₂O.day⁻¹), $\Delta t$ is the time step (1 day), and $aq_{shthr,q}$ is the threshold water level in the shallow aquifer for return flow to occur (mm H₂O).

The groundwater flow from the deep aquifer is represented by Equation (23), as follows:

$$ Q_{gw,dp,i} = W_{rchrg,dp,i} (1 - e^{-\alpha_{gw,dp} \times \Delta t}) + Q_{gw,dp,i-1} (e^{-\alpha_{gw,dp} \times \Delta t}) \quad (23) $$

where $Q_{gw,sh,i}$ is the groundwater flow from confined aquifer (mm H₂O.day⁻¹), $\Delta t$ is the time step (1 day), and $\alpha_{gw,dp}$ is the groundwater recession constant of the deep aquifer (days⁻¹).
increase in elevation per unit distance, $\emptyset_d$ is the drainable (residual) porosity of the soil layer (mm/mm), and $L_{hill}$ is the hillslope length (m).

The daily water balance for each soil layer is expressed using Equation (18), as follows [94]:

$$ \Delta SW_{ly} = W_{perc,ly-1} - W_{perc,ly} - Q_{lat,ly} - E_{e,ly} - E_{t,ly} \quad (18) $$

where $\Delta SW_{ly}$ is the change of soil water content at soil layer (ly), $W_{perc,ly-1}$ is the percolation received from layer (ly-1), $W_{perc,ly}$ and $Q_{lat,ly}$ are the percolation and lateral flow generated from soil layer (ly), respectively, and $E_{e,ly}$ and $E_{t,ly}$ are the evaporation and transpiration drawn from the soil layer (ly), respectively (all variables are expressed in mm H₂O.day⁻¹).

### 2.2.2. Groundwater Flow and Baseflow to the Stream

The groundwater module of SWAT comprises a system of two aquifers in each sub-basin: (1) a shallow unconfined aquifer that generates baseflow into the stream and (2) a deep confined aquifer contributing to streamflow outside of the watershed (flow lost from the system) [95]. Recharge from the unsaturated soil profile to the aquifers on a given day is calculated using an exponential decay weighting function that accounts for the time delay of the recharge mechanism, as follows [67]:

$$ W_{rchrg,i} = \left(1 - e^{-\frac{1}{\delta_{gw}}}\right) W_{seep,i} + \left(e^{-\frac{1}{\delta_{gw}}}\right) W_{rchrg,i-1} \quad (19) $$

where $W_{rchrg,i}$ and $W_{rchrg,i-1}$ represent the recharge to the aquifers (shallow and deep) at days i and i-1 (mm H₂O.day⁻¹), respectively, $W_{seep,i}$ is the water drained from the bottom of the soil profile (mm H₂O.day⁻¹), and $\delta_{gw}$ is the delay time required for recharge to reach the aquifers (days).

Recharge components routed to the shallow (unconfined) aquifer and the deep (confined) aquifer are computed using Equations (20) and (21), respectively, as follows:

$$ W_{rchrg,sh,i} = (1 - \beta_{dp}) W_{rchrg,i} \quad (20) $$

$$ W_{rchrg,dp,i} = \beta_{dp} W_{rchrg,i} \quad (21) $$

where $W_{rchrg,sh,i}$ and $W_{rchrg,dp,i}$ represent the water diverted to the shallow and deep aquifers (mm H₂O.day⁻¹), respectively, and $\beta_{dp}$ is a coefficient of percolation to the deep aquifer.

The shallow aquifer contributes to the streamflow if water stored in the aquifer exceeds a user-specified threshold. Otherwise, return flow is set to zero. The daily groundwater flow to the main river channel is computed using an exponential storage-discharge relationship, which incorporates the recharge from the shallow aquifer and a baseflow recession constant, as follows:

$$ Q_{gw,sh,i} = \begin{cases} W_{rchrg,sh,i} (1 - e^{-\alpha_{gw,sh} \times \Delta t}) + Q_{gw,sh,i-1} (e^{-\alpha_{gw,sh} \times \Delta t}), & aq_{sh} > aq_{shthr,q} \\ 0, & aq_{sh} \le aq_{shthr,q} \end{cases} \quad (22) $$

where $Q_{gw,sh,i}$ is the baseflow from the shallow aquifer to the main stream channel (mm H₂O.day⁻¹), $\alpha_{gw,sh}$ is the groundwater recession constant of shallow aquifer (days⁻¹), $aq_{sh}$ is the amount of water stored in the shallow aquifer (mm H₂O.day⁻¹), $\Delta t$ is the time step (1 day), and $aq_{shthr,q}$ is the threshold water level in the shallow aquifer for return flow to occur (mm H₂O).

The groundwater flow from the deep aquifer is represented by Equation (23), as follows:

$$ Q_{gw,dp,i} = W_{rchrg,dp,i} (1 - e^{-\alpha_{gw,dp} \times \Delta t}) + Q_{gw,dp,i-1} (e^{-\alpha_{gw,dp} \times \Delta t}) \quad (23) $$

where $Q_{gw,sh,i}$ is the groundwater flow from confined aquifer (mm H₂O.day⁻¹), $\Delta t$ is the time step (1 day), and $\alpha_{gw,dp}$ is the groundwater recession constant of the deep aquifer (days⁻¹).
In dry periods, water in the shallow aquifers may be removed by evaporation to the partially saturated overlaying soil through the capillary fringe that separates the saturated and vadose zones. Water can also be directly absorbed by deep rooted plants through transpiration [96]. SWAT accounts for this phenomenon via a process defined as revap, which occurs when water storage in the shallow aquifer exceeds a user-defined threshold. The amount of water that can be potentially consumed by revap is calculated as follows [97]:

$$ET_{rvp-max} = \beta_{rvp} \times PET \qquad (24)$$

where $ET_{rvp,max}$ is the maximum amount of water that can be removed from the shallow aquifer (mm H₂O.day⁻¹), $\beta_{rvp}$ is the groundwater evaporation coefficient, and PET is the potential evapotranspiration (mm H₂O.day⁻¹). The actual groundwater evapotranspiration $ET_{rvp}$ is subsequently calculated based on water availability in the shallow aquifer, considering the following cases [69]:

$$ET_{rvp} = \begin{cases} ET_{rvp,max} - aq_{shthr,rvp}, & 0, \quad aq_{sh} \le aq_{shthr,rvp} \\ ET_{rvp,max} - aq_{shthr,rvp}, & aq_{shthr,rvp} < aq_{sh} < (aq_{shthr,rvp} + ET_{rvp,max}) \\ ET_{rvp,max}, & aq_{sh} \ge (aq_{shthr,rvp} + ET_{rvp,max}) \end{cases} \qquad (25)$$

where $aq_{sh}$ is the water stored in the shallow aquifer at the beginning of day i (mm H₂O.day⁻¹) and $aq_{shthr,rvp}$ is the threshold water level in the shallow aquifer for groundwater evaporation to occur.

The volumetric water balance for the shallow aquifer is represented as follows [98]:

$$aq_{sh,i} = aq_{sh,i-1} + W_{rchrg,sh,i} - Q_{gw,sh,i} - ET_{rvp,i} - Q_{pump,sh,i} \qquad (26)$$

where $aq_{sh,i}$ and $aq_{sh,i-1}$ represent water stored in the shallow aquifer on days i and i − 1, respectively, $ET_{rvp,i}$ is the volume of water that moves upward by capillary rise, and $Q_{pump,sh,i}$ is the water withdrawn by pumping from the shallow aquifer (all variables are expressed in mm H₂O.day⁻¹).

SWAT also simulates other types of water bodies, including wetlands, ponds, and depressions or potholes. These water bodies are modeled within the subbasins of the main stream channel and are fed by runoff originating from the subbasin in which they are located [99]. They can also contribute to seepage and groundwater recharge, adding to the recharge from soil water percolation [100].

The downward daily seepage from the pond or wetland $V_{seep}$ (m³ H₂O.day⁻¹) is estimated using Equation (27) [69]:

$$V_{seep} = 240 \times K_{sat} \times A_{wet} \qquad (27)$$

where $K_{sat}$ is the saturated hydraulic conductivity of the pond or wetland bottom (mm.h⁻¹) and $A_{wet}$ is the water surface area of the pond or wetland (hectares).

Daily seepage from the pothole/depression is computed as a function of soil water content, as follows [101]:

$$V_{seep,pot} = \begin{cases} 240K_{sat} \times SA, & SW < 0.5FC \\ 240\left(1 - \frac{SW}{FC}\right) \times K_{sat} \times SA, & 0.5FC \le SW < FC \\ 0, & SW \ge FC \end{cases} \qquad (28)$$

where $V_{seep,pot}$ is the seepage from a pothole (m³ H₂O.day⁻¹), $K_{sat}$ is the saturated hydraulic conductivity of the top soil layer (mm.h⁻¹), $SA$ is the pothole surface area (hectares), $SW$ is the daily soil water content of the profile (mm H₂O), and $FC$ is the field capacity moisture content (mm H₂O).
In dry periods, water in the shallow aquifers may be removed by evaporation to the partially saturated overlaying soil through the capillary fringe that separates the saturated and vadose zones. Water can also be directly absorbed by deep rooted plants through transpiration [96]. SWAT accounts for this phenomenon via a process defined as revap, which occurs when water storage in the shallow aquifer exceeds a user-defined threshold. The amount of water that can be potentially consumed by revap is calculated as follows [97]:

$$ET_{rvp-max} = \beta_{rvp} \times PET \quad (24)$$

where $ET_{rvp,max}$ is the maximum amount of water that can be removed from the shallow aquifer (mm H₂O.day⁻¹), $\beta_{rvp}$ is the groundwater evaporation coefficient, and PET is the potential evapotranspiration (mm H₂O.day⁻¹). The actual groundwater evapotranspiration $ET_{rvp}$ is subsequently calculated based on water availability in the shallow aquifer, considering the following cases [69]:

$$ET_{rvp} = \begin{cases} ET_{rvp,max} - aq_{shthr,rvp}, & 0, \quad aq_{sh} \le aq_{shthr,rvp} \\ ET_{rvp,max} - aq_{shthr,rvp}, & aq_{shthr,rvp} < aq_{sh} < (aq_{shthr,rvp} + ET_{rvp,max}) \\ ET_{rvp,max}, & aq_{sh} \ge (aq_{shthr,rvp} + ET_{rvp,max}) \end{cases} \quad (25)$$

where $aq_{sh}$ is the water stored in the shallow aquifer at the beginning of day i (mm H₂O.day⁻¹) and $aq_{shthr,rvp}$ is the threshold water level in the shallow aquifer for groundwater evaporation to occur.

The volumetric water balance for the shallow aquifer is represented as follows [98]:

$$aq_{sh,i} = aq_{sh,i-1} + W_{rchrg,sh,i} - Q_{gw,sh,i} - ET_{rvp,i} - Q_{pump,sh,i} \quad (26)$$

where $aq_{sh,i}$ and $aq_{sh,i-1}$ represent water stored in the shallow aquifer on days i and i − 1, respectively, $ET_{rvp,i}$ is the volume of water that moves upward by capillary rise, and $Q_{pump,sh,i}$ is the water withdrawn by pumping from the shallow aquifer (all variables are expressed in mm H₂O.day⁻¹).

SWAT also simulates other types of water bodies, including wetlands, ponds, and depressions or potholes. These water bodies are modeled within the subbasins of the main stream channel and are fed by runoff originating from the subbasin in which they are located [99]. They can also contribute to seepage and groundwater recharge, adding to the recharge from soil water percolation [100].

The downward daily seepage from the pond or wetland $V_{seep}$ (m³ H₂O.day⁻¹) is estimated using Equation (27) [69]:

$$V_{seep} = 240 \times K_{sat} \times A_{wet} \quad (27)$$

where $K_{sat}$ is the saturated hydraulic conductivity of the pond or wetland bottom (mm.h⁻¹) and $A_{wet}$ is the water surface area of the pond or wetland (hectares).

Daily seepage from the pothole/depression is computed as a function of soil water content, as follows [101]:

$$V_{seep,pot} = \begin{cases} 240K_{sat} \times SA, & SW < 0.5FC \\ 240\left(1 - \frac{SW}{FC}\right) \times K_{sat} \times SA, & 0.5FC \le SW < FC \\ 0, & SW \ge FC \end{cases} \quad (28)$$

where $V_{seep,pot}$ is the seepage from a pothole (m³ H₂O.day⁻¹), $K_{sat}$ is the saturated hydraulic conductivity of the top soil layer (mm.h⁻¹), $SA$ is the pothole surface area (hectares), $SW$ is the daily soil water content of the profile (mm H₂O), and $FC$ is the field capacity moisture content (mm H₂O).
### 3. SWAT Studies in Karst Watersheds: Selection and Classification Methods

We used the SWAT Literature Database (CARD) [38] and Google Scholar engine to identify SWAT research studies in karst watersheds, published between the years 2000 (the year that the first SWAT study in a karst watershed was published) and 2022. Searching priority was initially accorded to the 5400+ articles available in CARD and grouped by specific application categories. All SWAT code iterations (standard and modified) were included in the search and selection process of the articles, based on the keywords “hydro-logic”, “hydrologic and pollutants”, and “karst”. Consequently, 17 articles were identified in CARD. Then, multiple searches were performed using Google Scholar to identify the studies that have not been included in CARD, considering the above-mentioned criteria terms in combination with the term “SWAT”. Only peer-reviewed articles and published thesis reports in Google Scholar were selected for further assessment, whereas technical reports, abstracts/conference papers, and non-English articles were excluded. Combining both literature databases, a total of 75 studies related to SWAT simulations in karstic and partially karstified watersheds were identified. We classified these studies into two main categories: (1) the standard SWAT model applications (category I) and (2) the coupled/modified SWAT model applications (category II). Subsequently, 25 studies reporting an application of a modified SWAT or SWAT coupled with a karstic flow model fell under category II, while the remaining 50 studies fell under the first category I.

In this paper, we grouped the articles under category I by region (North and Latin America, Europe, Asia, and Africa) and study scope (i.e., hydrological or water quality modeling, climate or land-use change impacts) (Table 1). For the sake of paper length, we discussed the studies under category I that presented a novel simulation approach or a complex application of the standard SWAT in karst watersheds. Next, we subdivided the articles under category II based upon: (1) the conceptual models/algorithms coupled with SWAT or used to modify the SWAT source code, (2) the studied karst processes/features (e.g., matrix, conduits, springs, sinkholes), and (3) the simulation scope (e.g., hydrological or water quality modeling, climate or land-use change impacts) (Table 2). Then, we thoroughly presented the core methodology and major findings of the SWAT studies of category II, which focused primarily on hydrological simulation. Appendices A–K summarize the equations of the karstic models coupled with SWAT and used in the different modified variants of the code. Finally, we identified potential constraints of the modified SWAT models so that they can so that they can be considered in developing future SWAT models adapted to karst hydrology.

The accuracy of the SWAT models’ outputs was reported in their respective studies using different statistical indicators, such as the Nash-Sutcliffe Efficiency (NSE), the coefficient of determination (R²), the percent of bias (PBIAS%), the root mean square error observations standard deviation ratio (RSR), and the Kling-Gupta Efficiency (KGE) [102,103]. In this review, the overall trends of the hydrological models’ performance were examined using NSE, being the most commonly applied statistical indicator across all the reported studies. NSE is a measure of the relative magnitude of the residual variance against the observed data variance. It is used to assess the goodness of fit of the plot of observed versus simulated data, and is computed as follows [102]:

$$NSE = 1 - \frac{\sum_{i=1}^{n} (O_i - S_i)^2}{\sum_{i=1}^{n} (O_i - \bar{O})^2} \quad (29)$$

where $O_i$ and $S_i$ represent the ith value of the observed and simulated data, respectively, $\bar{O}$ is the mean of the observed and simulated data, and n is the total number of observations. NSE values can vary between -∞ and 1. In particular, watershed streamflow simulation at the daily, monthly, and annual scales is judged as satisfactory if 0.5 < NSE ≤ 0.7, good if 0.7 < NSE ≤ 0.8, and very good for NSE ≥ 0.8. Conversely it is unsatisfactory if NSE ≤ 0.5, while negative NSE values indicate an unacceptable model performance [102].
### 3. SWAT Studies in Karst Watersheds: Selection and Classification Methods

We used the SWAT Literature Database (CARD) [38] and Google Scholar engine to identify SWAT research studies in karst watersheds, published between the years 2000 (the year that the first SWAT study in a karst watershed was published) and 2022. Searching priority was initially accorded to the 5400+ articles available in CARD and grouped by specific application categories. All SWAT code iterations (standard and modified) were included in the search and selection process of the articles, based on the keywords “hydro-logic”, “hydrologic and pollutants”, and “karst”. Consequently, 17 articles were identified in CARD. Then, multiple searches were performed using Google Scholar to identify the studies that have not been included in CARD, considering the above-mentioned criteria terms in combination with the term “SWAT”. Only peer-reviewed articles and published thesis reports in Google Scholar were selected for further assessment, whereas technical reports, abstracts/conference papers, and non-English articles were excluded. Combining both literature databases, a total of 75 studies related to SWAT simulations in karstic and partially karstified watersheds were identified. We classified these studies into two main categories: (1) the standard SWAT model applications (category I) and (2) the coupled/modified SWAT model applications (category II). Subsequently, 25 studies reporting an application of a modified SWAT or SWAT coupled with a karstic flow model fell under category II, while the remaining 50 studies fell under the first category I.

In this paper, we grouped the articles under category I by region (North and Latin America, Europe, Asia, and Africa) and study scope (i.e., hydrological or water quality modeling, climate or land-use change impacts) (Table 1). For the sake of paper length, we discussed the studies under category I that presented a novel simulation approach or a complex application of the standard SWAT in karst watersheds. Next, we subdivided the articles under category II based upon: (1) the conceptual models/algorithms coupled with SWAT or used to modify the SWAT source code, (2) the studied karst processes/features (e.g., matrix, conduits, springs, sinkholes), and (3) the simulation scope (e.g., hydrological or water quality modeling, climate or land-use change impacts) (Table 2). Then, we thoroughly presented the core methodology and major findings of the SWAT studies of category II, which focused primarily on hydrological simulation. Appendices A–K summarize the equations of the karstic models coupled with SWAT and used in the different modified variants of the code. Finally, we identified potential constraints of the modified SWAT models so that they can so that they can be considered in developing future SWAT models adapted to karst hydrology.

The accuracy of the SWAT models’ outputs was reported in their respective studies using different statistical indicators, such as the Nash–Sutcliffe Efficiency (NSE), the coefficient of determination (R²), the percent of bias (PBIAS%), the root mean square error observations standard deviation ratio (RSR), and the Kling–Gupta Efficiency (KGE) [102,103]. In this review, the overall trends of the hydrological models’ performance were examined using NSE, being the most commonly applied statistical indicator across all the reported studies. NSE is a measure of the relative magnitude of the residual variance against the observed data variance. It is used to assess the goodness of fit of the plot of observed versus simulated data, and is computed as follows [102]:

$$NSE = 1 - \frac{\sum_{i=1}^{n} (O_i - S_i)^2}{\sum_{i=1}^{n} (O_i - \bar{O})^2} \qquad (29)$$

where $O_i$ and $S_i$ represent the ith value of the observed and simulated data, respectively, $\bar{O}$ is the mean of the observed and simulated data, and n is the total number of observations. NSE values can vary between -∞ and 1. In particular, watershed streamflow simulation at the daily, monthly, and annual scales is judged as satisfactory if 0.5 < NSE ≤ 0.7, good if 0.7 < NSE ≤ 0.8, and very good for NSE ≥ 0.8. Conversely it is unsatisfactory if NSE ≤ 0.5, while negative NSE values indicate an unacceptable model performance [102].
<table><caption>**Table 1.** Reference, basin description, and application of the standard SWAT studies in karst watersheds (category I).</caption><thead><tr><th>Region</th><th>Reference</th><th>Basin Name (Country, Size in km²)</th><th>Application</th></tr></thead><tbody><tr><td rowspan="15">North and Latin America</td><td>Spruill et al. [104]</td><td>University of KY Research Site (USA; 5.5)</td><td>Simulation of streamflow</td></tr><tr><td>Coffey et al. [105]</td><td>(University of KY Research Site (USA; 5.5)</td><td>Simulation of streamflow</td></tr><tr><td>Benham et al. [106]</td><td>Shoal Creek (USA; 367)</td><td>Simulation of streamflow and bacteria fate and transport</td></tr><tr><td>Amatya and Jha [107]</td><td>Chapel Branch Creek (USA; 15.55)</td><td>Simulation of streamflow in a watershed with a flooded embayment outlet draining to a lake</td></tr><tr><td>Amatya and Jha [108]</td><td>Chapel Branch Creek (USA; 15.55)</td><td>Simulation of streamflow and phosphorus loads and concentrations in karst watershed tributaries and downstream a reservoir-like embayment outlet</td></tr><tr><td>Williams et al. [109]</td><td>Chapel Branch Creek (USA; 15.55)</td><td>Simulation of streamflow, nitrogen loads, and phosphorus loads in a karst watershed draining to a lake via a reservoir-like embayment</td></tr><tr><td>Wilson et al. [110]</td><td>South Branch, Root River (USA; 301.8)</td><td>Impacts of traditional and alternative conservation management practices on water quality (sediments and phosphorus)</td></tr><tr><td>Jain et al. [111]</td><td>Nueces River Headwaters (USA; 2126)</td><td>Impacts of land-use/cover change on watershed hydrology</td></tr><tr><td>Sunde et al. [112]</td><td>Hinkson Creek (USA; 231)</td><td>Impacts of future urban development on watershed hydrology</td></tr><tr><td>Sunde et al. [113]</td><td>Hinkson Creek (USA; 231)</td><td>Impacts of climate change on watershed hydrological processes</td></tr><tr><td>Sunde et al. [114]</td><td>Hinkson Creek (USA; 231)</td><td>Impacts of future urbanization and climate change on watershed hydrology</td></tr><tr><td>Sarkar et al. [115]</td><td>Conestoga River (USA; 1230)</td><td>Simulation of flow, sediment loads from upland watershed sources, flow routing, and sediment processes using a coupled SWAT-HSPF model</td></tr><tr><td>Merriman et al. [116]</td><td>Upper East River (USA; 375.3)</td><td>Impacts of agricultural best management practices on flow, sediment loads, and nutrient loads</td></tr><tr><td>Sullivan et al. [117]</td><td>Edwards aquifer overlain by Cibolo Creek watershed (USA; 707) and Dry Comal Creek watershed (USA; 337)</td><td>Simulation of nitrate concentration inputs to MODFLOW CFPv2 and CMT3D models used to assess nitrate transport in an aquifer</td></tr><tr><td>Chen et al. [118]</td><td>Blanco River (N/A)</td><td>Multi-model projections of hydrological drought characteristics under climate change</td></tr><tr><td>Zeiger et al. [5]</td><td>James River (USA; 3770)</td><td>Impacts of climate and land use on streamflow, sediment, and nutrient loads, and identification of critical source areas of non-point source pollution</td></tr><tr><td>Al Aamery et al. [119]</td><td>Cane Run-Royal Spring (58)</td><td>Simulation of surface runoff, surface routing, and soil water percolation inputs for a fluviokarst-specific combined discrete continuum numerical model</td></tr><tr><td>Karki et al. [120]</td><td>Apalachicola-Chattahoochee-Flint River (USA; 12,000)</td><td>Simulation of groundwater areal recharge input for a MODFLOW-NWT aquifer model</td></tr></tbody></table>
<table><caption>**Table 1.** Reference, basin description, and application of the standard SWAT studies in karst watersheds (category I).</caption><thead><tr><th>Region</th><th>Reference</th><th>Basin Name (Country, Size in km²)</th><th>Application</th></tr></thead><tbody><tr><td rowspan="15">North and Latin America</td><td>Spruill et al. [104]</td><td>University of KY Research Site (USA; 5.5)</td><td>Simulation of streamflow</td></tr><tr><td>Coffey et al. [105]</td><td>(University of KY Research Site (USA; 5.5)</td><td>Simulation of streamflow</td></tr><tr><td>Benham et al. [106]</td><td>Shoal Creek (USA; 367)</td><td>Simulation of streamflow and bacteria fate and transport</td></tr><tr><td>Amatya and Jha [107]</td><td>Chapel Branch Creek (USA; 15.55)</td><td>Simulation of streamflow in a watershed with a flooded embayment outlet draining to a lake</td></tr><tr><td>Amatya and Jha [108]</td><td>Chapel Branch Creek (USA; 15.55)</td><td>Simulation of streamflow and phosphorus loads and concentrations in karst watershed tributaries and downstream a reservoir-like embayment outlet</td></tr><tr><td>Williams et al. [109]</td><td>Chapel Branch Creek (USA; 15.55)</td><td>Simulation of streamflow, nitrogen loads, and phosphorus loads in a karst watershed draining to a lake via a reservoir-like embayment</td></tr><tr><td>Wilson et al. [110]</td><td>South Branch, Root River (USA; 301.8)</td><td>Impacts of traditional and alternative conservation management practices on water quality (sediments and phosphorus)</td></tr><tr><td>Jain et al. [111]</td><td>Nueces River Headwaters (USA; 2126)</td><td>Impacts of land-use/cover change on watershed hydrology</td></tr><tr><td>Sunde et al. [112]</td><td>Hinkson Creek (USA; 231)</td><td>Impacts of future urban development on watershed hydrology</td></tr><tr><td>Sunde et al. [113]</td><td>Hinkson Creek (USA; 231)</td><td>Impacts of climate change on watershed hydrological processes</td></tr><tr><td>Sunde et al. [114]</td><td>Hinkson Creek (USA; 231)</td><td>Impacts of future urbanization and climate change on watershed hydrology</td></tr><tr><td>Sarkar et al. [115]</td><td>Conestoga River (USA; 1230)</td><td>Simulation of flow, sediment loads from upland watershed sources, flow routing, and sediment processes using a coupled SWAT-HSPF model</td></tr><tr><td>Merriman et al. [116]</td><td>Upper East River (USA; 375.3)</td><td>Impacts of agricultural best management practices on flow, sediment loads, and nutrient loads</td></tr><tr><td>Sullivan et al. [117]</td><td>Edwards aquifer overlain by Cibolo Creek watershed (USA; 707) and Dry Comal Creek watershed (USA; 337)</td><td>Simulation of nitrate concentration inputs to MODFLOW CFPv2 and CMT3D models used to assess nitrate transport in an aquifer</td></tr><tr><td>Chen et al. [118]</td><td>Blanco River (N/A)</td><td>Multi-model projections of hydrological drought characteristics under climate change</td></tr><tr><td>Zeiger et al. [5]</td><td>James River (USA; 3770)</td><td>Impacts of climate and land use on streamflow, sediment, and nutrient loads, and identification of critical source areas of non-point source pollution</td></tr><tr><td>Al Aamery et al. [119]</td><td>Cane Run-Royal Spring (58)</td><td>Simulation of surface runoff, surface routing, and soil water percolation inputs for a fluviokarst-specific combined discrete continuum numerical model</td></tr><tr><td>Karki et al. [120]</td><td>Apalachicola-Chattahoochee-Flint River (USA; 12,000)</td><td>Simulation of groundwater areal recharge input for a MODFLOW-NWT aquifer model</td></tr></tbody></table>
<table><caption>Table 1. *Cont.*</caption><thead><tr><th>Region</th><th>Reference</th><th>Basin Name (Country, Size in km²)</th><th>Application</th></tr></thead><tbody><tr><td rowspan="10">Europe</td><td>Salerno and Tartari [121]</td><td>Subbasin of the Lake Pusiano watershed (Italy, 52.5)</td><td>Simulation of discharge using SWAT supported by wavelet analysis to assess the contribution of external flow component to streamflow</td></tr><tr><td>Vale and Holman [122]</td><td>Bosherston Lakes (UK; N/A)</td><td>Quantitative assessment of the hydrological processes controlling water levels and groundwater-surface water interactions in a lake system</td></tr><tr><td>Tzoraki et al. [123]</td><td>Evrotas (Greece; 2050)</td><td>Simulation and analysis of flood events characteristics</td></tr><tr><td>Palazón and Navas [124]</td><td>Linsoles River (Spain; 284)</td><td>Simulation of surface runoff and sediment yield</td></tr><tr><td>Palazón and Navas [125]</td><td>The Barasona reservoir catchment (Spain; 1509)</td><td>Simulation of erosion and sediment yield</td></tr><tr><td>Sellami et al. [126]</td><td>Thau catchment (France; 280)</td><td>Assessment of SWAT model accuracy in predicting discharge at gauged and ungauged catchments within an uncertainty framework</td></tr><tr><td>Gamvroudis et al. [127]</td><td>Evrotas River (Greece; 1348)</td><td>Simulation of watershed water budget and spatial distribution of runoff and sediment transport</td></tr><tr><td>Malagò et al. [128]</td><td>Scandanavian Peninsula (10⁶); Iberian Peninsula (556,000)</td><td>Hydrological simulation, sensitivity analysis, multi-variable calibration, and regionalization of the calibrated parameters for the identification of dominant hydrological processes in each region</td></tr><tr><td>Mehdi et al. [62]</td><td>Altmühl River (Germany; 980)</td><td>Impacts of climate and land-use changes on streamflow and nutrients loads</td></tr><tr><td>Sellami et al. [129]</td><td>Thau catchment (France; 280)</td><td>Impacts of climate change on watershed hydrology</td></tr><tr><td>Palazón and Navas [130]</td><td>The Barasona reservoir catchment (Spain; 1509)</td><td>Simulation of streamflow under different precipitation characterization scenarios</td></tr><tr><td>Vigiak et al. [131]</td><td>Danube River (800,000)</td><td>Simulation of sediment fluxes under soil conservation measures and identification of sediment budget knowledge gaps</td></tr><tr><td>Efthimiou [132]</td><td>Kalamas River (Greece; 1899.25)</td><td>Simulation of watershed hydrological budget</td></tr><tr><td>Martínez-Salvador and Conesa-García [133]</td><td>Upper Argos River (Spain; 510)</td><td>Simulation of streamflow and sediment load</td></tr><tr><td>Senent-Aparicio et al. [134]</td><td>Castril River (Spain; 120)</td><td>Simulation of streamflow using SWAT supported by chloride mass balance to estimate IGF contribution to streamflow</td></tr><tr><td>Busico et al. [135]</td><td>Anthemountas (Greece; 374)</td><td>Assessment of groundwater recharge variations and their relationship with other hydrological parameters under climate change</td></tr><tr><td>Sánchez-Gómez et al. [136]</td><td>Henares River (Spain; 4070)</td><td>Optimization of SWAT streamflow simulation by incorporating watershed geological properties in model calibration</td></tr></tbody></table>
<table><caption>Table 1. *Cont.*</caption><thead><tr><th>Region</th><th>Reference</th><th>Basin Name (Country, Size in km²)</th><th>Application</th></tr></thead><tbody><tr><td rowspan="13">Europe</td><td>Salerno and Tartari [121]</td><td>Subbasin of the Lake Pusiano watershed (Italy, 52.5)</td><td>Simulation of discharge using SWAT supported by wavelet analysis to assess the contribution of external flow component to streamflow</td></tr><tr><td>Vale and Holman [122]</td><td>Bosherston Lakes (UK; N/A)</td><td>Quantitative assessment of the hydrological processes controlling water levels and groundwater-surface water interactions in a lake system</td></tr><tr><td>Tzoraki et al. [123]</td><td>Evrotas (Greece; 2050)</td><td>Simulation and analysis of flood events characteristics</td></tr><tr><td>Palazón and Navas [124]</td><td>Linsoles River (Spain; 284)</td><td>Simulation of surface runoff and sediment yield</td></tr><tr><td>Palazón and Navas [125]</td><td>The Barasona reservoir catchment (Spain; 1509)</td><td>Simulation of erosion and sediment yield</td></tr><tr><td>Sellami et al. [126]</td><td>Thau catchment (France; 280)</td><td>Assessment of SWAT model accuracy in predicting discharge at gauged and ungauged catchments within an uncertainty framework</td></tr><tr><td>Gamvroudis et al. [127]</td><td>Evrotas River (Greece; 1348)</td><td>Simulation of watershed water budget and spatial distribution of runoff and sediment transport</td></tr><tr><td>Malagò et al. [128]</td><td>Scandanavian Peninsula (10⁶); Iberian Peninsula (556,000)</td><td>Hydrological simulation, sensitivity analysis, multi-variable calibration, and regionalization of the calibrated parameters for the identification of dominant hydrological processes in each region</td></tr><tr><td>Mehdi et al. [62]</td><td>Altmühl River (Germany; 980)</td><td>Impacts of climate and land-use changes on streamflow and nutrients loads</td></tr><tr><td>Sellami et al. [129]</td><td>Thau catchment (France; 280)</td><td>Impacts of climate change on watershed hydrology</td></tr><tr><td>Palazón and Navas [130]</td><td>The Barasona reservoir catchment (Spain; 1509)</td><td>Simulation of streamflow under different precipitation characterization scenarios</td></tr><tr><td>Vigiak et al. [131]</td><td>Danube River (800,000)</td><td>Simulation of sediment fluxes under soil conservation measures and identification of sediment budget knowledge gaps</td></tr><tr><td>Efthimiou [132]</td><td>Kalamas River (Greece; 1899.25)</td><td>Simulation of watershed hydrological budget</td></tr><tr><td>Martínez-Salvador and Conesa-García [133]</td><td>Upper Argos River (Spain; 510)</td><td>Simulation of streamflow and sediment load</td></tr><tr><td>Senent-Aparicio et al. [134]</td><td>Castril River (Spain; 120)</td><td>Simulation of streamflow using SWAT supported by chloride mass balance to estimate IGF contribution to streamflow</td></tr><tr><td>Busico et al. [135]</td><td>Anthemountas (Greece; 374)</td><td>Assessment of groundwater recharge variations and their relationship with other hydrological parameters under climate change</td></tr><tr><td>Sánchez-Gómez et al. [136]</td><td>Henares River (Spain; 4070)</td><td>Optimization of SWAT streamflow simulation by incorporating watershed geological properties in model calibration</td></tr></tbody></table>
<table><caption>Table 1. *Cont.*</caption><thead><tr><th>Region</th><th>Reference</th><th>Basin Name (Country, Size in km²)</th><th>Application</th></tr></thead><tbody><tr><td rowspan="11">Asia</td><td>Jiang et al. [137]</td><td>Shibetsu River (Japan; 672)</td><td>Simulation of streamflow and external flow contribution to discharge from the water balance equation, using measured data</td></tr><tr><td>Tian et al. [138]</td><td>Shibantang River (China; 2248)</td><td>Assessment of trade-offs and synergic relationships between ecosystem services (water yield, sediment yield, and net primary productivity)</td></tr><tr><td>Bucak et al. [139]</td><td>Lake Beyşehir catchment (Turkey; 4704)</td><td>Impacts of climate and land-use changes on the hydrological balance of a lake catchment and water levels</td></tr><tr><td>Hou and Gao [140]</td><td>Sancha River (China, 4068)<fnref n="1" /></td><td>Simulation of the spatial variability of streamflow, surface runoff, and groundwater runoff, and analysis of their spatial correlation with environmental factors</td></tr><tr><td>Jakada and Chen [141]</td><td>Miaogou subbasin of Gaolan River Basin (China; 45)</td><td>Simulation of watershed hydrology using SWAT supported by a geological survey and a tracer test</td></tr><tr><td>Mo et al. [142]</td><td>Xiajia River (China; 799.2)</td><td>Simulation of watershed runoff under different precipitation input data</td></tr><tr><td>Hou et al. [143]</td><td>Guizhou Province (China 4681)</td><td>Analysis of the factors affecting streamflow, surface runoff, and groundwater, and their interactions for different geomorphic types</td></tr><tr><td>Gao et al. [144]</td><td>Sancha River (China, 7061)<fnref n="1" /></td><td>Assessment of trade-offs and synergic relationships between ecosystem services (sediment yield and surface/slope runoff, water yield, and slope runoff) and main factors affecting their relationships, for different geomorphic types</td></tr><tr><td>Jiang et al. [145]</td><td>Sancha River (China, 7061)<fnref n="1" /></td><td>Simulation of the spatial distributions of rainfall erosivity and runoff erosivity, and identification of the dominant factors and their interactions affecting the spatial distributions of rainfall/runoff erosivity, for different geomorphic types</td></tr><tr><td>Chang et al. [146]</td><td>Nanpan River (China; 43,200)</td><td>Simulation of soil moisture using SWAT and development of a methodology for a comprehensive drought index based on the watershed hydrological processes (precipitation, runoff, and soil moisture)</td></tr><tr><td>Zhang et al. [147]</td><td>Lijiang River (China, 5444)</td><td>Simulation of streamflow and water quality using SWAT and HSPF models driven by different precipitation input data, and impacts of best management practices on non-point-source pollution reduction</td></tr><tr><td>Mo et al. [148]</td><td>Chengbi River (China; 2087)</td><td>Simulation of runoff under different calibration methods and precipitation input data</td></tr><tr><td>Yuan et al. [149]</td><td>Gaoche catchment area of the Dabang River basin (China; 1877.20)</td><td>Assessment of trade-offs and synergic relationships between ecosystem services (surface/underground runoff and surface sediment yield) and driving factors affecting their variation</td></tr><tr><td rowspan="2">Africa</td><td>Zettam et al. [150]</td><td>Tafna watershed (Algeria; 7245)</td><td>Simulation of watershed hydrological processes and assessment of the impacts of dam construction on water balance and sediment flux</td></tr><tr><td>Zaibak and Meddi [151]</td><td>Cheliff basin (Algeria; 43,750)</td><td>Simulation of streamflow at watershed dam-feeding subbasins and outlet</td></tr></tbody></table>

<footnote n="1">¹ There is a variation in the area of the Sancha River basin reported by Hou and Gao [140] compared to Gao et al. [144] and Jiang et al. [145].</footnote>
<table><caption>Table 1. *Cont.*</caption><thead><tr><th>Region</th><th>Reference</th><th>Basin Name (Country, Size in km²)</th><th>Application</th></tr></thead><tbody><tr><td rowspan="10">Asia</td><td>Jiang et al. [137]</td><td>Shibetsu River (Japan; 672)</td><td>Simulation of streamflow and external flow contribution to discharge from the water balance equation, using measured data</td></tr><tr><td>Tian et al. [138]</td><td>Shibantang River (China; 2248)</td><td>Assessment of trade-offs and synergic relationships between ecosystem services (water yield, sediment yield, and net primary productivity)</td></tr><tr><td>Bucak et al. [139]</td><td>Lake Beyşehir catchment (Turkey; 4704)</td><td>Impacts of climate and land-use changes on the hydrological balance of a lake catchment and water levels</td></tr><tr><td>Hou and Gao [140]</td><td>Sancha River (China, 4068) <sup>1</sup></td><td>Simulation of the spatial variability of streamflow, surface runoff, and groundwater runoff, and analysis of their spatial correlation with environmental factors</td></tr><tr><td>Jakada and Chen [141]</td><td>Miaogou subbasin of Gaolan River Basin (China; 45)</td><td>Simulation of watershed hydrology using SWAT supported by a geological survey and a tracer test</td></tr><tr><td>Mo et al. [142]</td><td>Xiajia River (China; 799.2)</td><td>Simulation of watershed runoff under different precipitation input data</td></tr><tr><td>Hou et al. [143]</td><td>Guizhou Province (China 4681)</td><td>Analysis of the factors affecting streamflow, surface runoff, and groundwater, and their interactions for different geomorphic types</td></tr><tr><td>Gao et al. [144]</td><td>Sancha River (China, 7061) <sup>1</sup></td><td>Assessment of trade-offs and synergic relationships between ecosystem services (sediment yield and surface/slope runoff, water yield, and slope runoff) and main factors affecting their relationships, for different geomorphic types</td></tr><tr><td>Jiang et al. [145]</td><td>Sancha River (China, 7061) <sup>1</sup></td><td>Simulation of the spatial distributions of rainfall erosivity and runoff erosivity, and identification of the dominant factors and their interactions affecting the spatial distributions of rainfall/runoff erosivity, for different geomorphic types</td></tr><tr><td>Chang et al. [146]</td><td>Nanpan River (China; 43,200)</td><td>Simulation of soil moisture using SWAT and development of a methodology for a comprehensive drought index based on the watershed hydrological processes (precipitation, runoff, and soil moisture)</td></tr><tr><td>Zhang et al. [147]</td><td>Lijiang River (China, 5444)</td><td>Simulation of streamflow and water quality using SWAT and HSPF models driven by different precipitation input data, and impacts of best management practices on non-point-source pollution reduction</td></tr><tr><td>Zhang et al. [147]</td><td>Chengbi River (China; 2087)</td><td>Simulation of runoff under different calibration methods and precipitation input data</td></tr><tr><td>Zhang et al. [147]</td><td>Gaoche catchment area of the Dabang River basin (China; 1877.20)</td><td>Assessment of trade-offs and synergic relationships between ecosystem services (surface/underground runoff and surface sediment yield) and driving factors affecting their variation</td></tr><tr><td rowspan="2">Africa</td><td>Zettam et al. [150]</td><td>Tafna watershed (Algeria; 7245)</td><td>Simulation of watershed hydrological processes and assessment of the impacts of dam construction on water balance and sediment flux</td></tr><tr><td>Zaibak and Meddi [151]</td><td>Cheliff basin (Algeria; 43,750)</td><td>Simulation of streamflow at watershed dam-feeding subbasins and outlet</td></tr></tbody></table>

<footnote n="1">¹ There is a variation in the area of the Sancha River basin reported by Hou and Gao [140] compared to Gao et al. [144] and Jiang et al. [145].</footnote>
<table><caption>**Table 2.** Groundwater modeling approach, reference, basin description and application of the modified SWAT codes in karst studies (category II).</caption><thead><tr><th>Groundwater Modeling Approach</th><th>Reference</th><th>Basin Name (Region; Size in km²)</th><th>Application—Modified SWAT Name (When Applicable)</th></tr></thead><tbody><tr><td rowspan="6">Conceptual linear one-reservoir groundwater model</td><td>Afinowicz et al. [152]</td><td>North Fork, Upper Guadalupe River (Texas-USA; 360)</td><td>Simulation of streamflow and water budget, and assessment crop management impacts on water budget</td></tr><tr><td>Baffaut and Benson [153]¹</td><td>James River (Missouri, USA; 3600)</td><td>Simulation of streamflow and pollutant transport (in-stream phosphorous loads and fecal coliform concentrations)—Adapted SWAT/SWAT-B&amp;B</td></tr><tr><td>Yactayo [154]¹</td><td>Opequon Creek (Virginia, USA; 890.2)</td><td>Simulation of streamflow and nitrate transport through the sinkholes in a karstic watershed—SWAT-karst</td></tr><tr><td>Palanisamy and Workman [155]¹</td><td>Cane Run Creek (Kentucky, USA; 115.6)</td><td>Simulation of streamflow through sinkholes in the streambed—KarstSWAT</td></tr><tr><td>Zhou et al. [156]¹</td><td>South and North Panjiang River (China; 2762)</td><td>Simulation of streamflow through sinkholes in the watershed</td></tr><tr><td>Nikolaidis et al. [157]</td><td>Koiliaris River (Crete, Greece; 132)</td><td>Simulation of water budget and in-stream nitrate concentrations, and assessment of climate change impacts on hydrology and water quality—Karst-SWAT</td></tr><tr><td rowspan="10">Conceptual linear two-reservoir groundwater model</td><td>Nerantzaki et al. [158]²</td><td>Koiliaris River (Crete, Greece; 130)</td><td>Simulation of suspended sediment transport, and assessment of climate change impacts on flow, soil erosion, and sediment transport</td></tr><tr><td>Tapoglou et al. [159]²</td><td>Crete Island (Greece; 8337)</td><td>Assessment of climate change impacts on extreme hydrometeorological events</td></tr><tr><td>Demetropoulou et al. [160]²</td><td>Geropotamos (Crete, Greece; 525 km²)</td><td>Methodology for the prioritization of a Program of Measures for water quantity and quality protection</td></tr><tr><td>Nerantzaki et al. [161]²</td><td>Crete Island (Greece; 8265)</td><td>Assessment of climate change impacts on hydrology</td></tr><tr><td>Lilli et al. [162]²</td><td>Koiliaris River (Crete, Greece; 130)</td><td>Development of erosion and flood protection nature-based solutions</td></tr><tr><td>Nerantzaki et al. [163]²</td><td>Koiliaris River (Crete, Greece; 130)</td><td>Uncertainty analysis of flow simulation due to the parameter uncertainty of the SWAT and Karst-SWAT models and internal variability of climate scenarios</td></tr><tr><td>Lilli et al. [164]²</td><td>Koiliaris River (Crete, Greece; 132)</td><td>Analysis of hydrological and geochemical processes</td></tr><tr><td>Malagò et al. [165]</td><td>Crete Island (Greece; 8336)</td><td>Simulation of hydrological water balance—KSWAT</td></tr><tr><td>Nguyen et al. [166]</td><td>Area in southwest Harz Mountains and southern Harz rim (Lower Saxony; Germany; 384)</td><td>Streamflow simulation, including IGF—SWAT_IGF</td></tr></tbody></table>
<table><caption>**Table 2.** Groundwater modeling approach, reference, basin description and application of the modified SWAT codes in karst studies (category II).</caption><thead><tr><th>Groundwater Modeling Approach</th><th>Reference</th><th>Basin Name (Region; Size in km²)</th><th>Application—Modified SWAT Name (When Applicable)</th></tr></thead><tbody><tr><td rowspan="6">Conceptual linear one-reservoir groundwater model</td><td>Afinowicz et al. [152]</td><td>North Fork, Upper Guadalupe River (Texas-USA; 360)</td><td>Simulation of streamflow and water budget, and assessment crop management impacts on water budget</td></tr><tr><td>Baffaut and Benson [153]¹</td><td>James River (Missouri, USA; 3600)</td><td>Simulation of streamflow and pollutant transport (in-stream phosphorous loads and fecal coliform concentrations)—Adapted SWAT/SWAT-B&amp;B</td></tr><tr><td>Yactayo [154]¹</td><td>Opequon Creek (Virginia, USA; 890.2)</td><td>Simulation of streamflow and nitrate transport through the sinkholes in a karstic watershed—SWAT-karst</td></tr><tr><td>Palanisamy and Workman [155]¹</td><td>Cane Run Creek (Kentucky, USA; 115.6)</td><td>Simulation of streamflow through sinkholes in the streambed—KarstSWAT</td></tr><tr><td>Zhou et al. [156]¹</td><td>South and North Panjiang River (China; 2762)</td><td>Simulation of streamflow through sinkholes in the watershed</td></tr><tr><td>Nikolaidis et al. [157]</td><td>Koiliaris River (Crete, Greece; 132)</td><td>Simulation of water budget and in-stream nitrate concentrations, and assessment of climate change impacts on hydrology and water quality—Karst-SWAT</td></tr><tr><td rowspan="10">Conceptual linear two-reservoir groundwater model</td><td>Nerantzaki et al. [158]²</td><td>Koiliaris River (Crete, Greece; 130)</td><td>Simulation of suspended sediment transport, and assessment of climate change impacts on flow, soil erosion, and sediment transport</td></tr><tr><td>Tapoglou et al. [159]²</td><td>Crete Island (Greece; 8337)</td><td>Assessment of climate change impacts on extreme hydrometeorological events</td></tr><tr><td>Demetropoulou et al. [160]²</td><td>Geropotamos (Crete, Greece; 525 km²)</td><td>Methodology for the prioritization of a Program of Measures for water quantity and quality protection</td></tr><tr><td>Nerantzaki et al. [161]²</td><td>Crete Island (Greece; 8265)</td><td>Assessment of climate change impacts on hydrology</td></tr><tr><td>Lilli et al. [162]²</td><td>Koiliaris River (Crete, Greece; 130)</td><td>Development of erosion and flood protection nature-based solutions</td></tr><tr><td>Nerantzaki et al. [163]²</td><td>Koiliaris River (Crete, Greece; 130)</td><td>Uncertainty analysis of flow simulation due to the parameter uncertainty of the SWAT and Karst-SWAT models and internal variability of climate scenarios</td></tr><tr><td>Lilli et al. [164]²</td><td>Koiliaris River (Crete, Greece; 132)</td><td>Analysis of hydrological and geochemical processes</td></tr><tr><td>Malagò et al. [165]</td><td>Crete Island (Greece; 8336)</td><td>Simulation of hydrological water balance—KSWAT</td></tr><tr><td>Nguyen et al. [166]</td><td>Area in southwest Harz Mountains and southern Harz rim (Lower Saxony; Germany; 384)</td><td>Streamflow simulation, including IGF—SWAT_IGF</td></tr></tbody></table>
<table><caption>Table 2. *Cont.*</caption><thead><tr><th>Groundwater Modeling Approach</th><th>Reference</th><th>Basin Name (Region; Size in km²)</th><th>Application—Modified SWAT Name (When Applicable)</th></tr></thead><tbody><tr><td>Conceptual linear three-reservoir groundwater model</td><td>Wang et al. [167]<fnref n="1" /><br>Geng et al. [168]<fnref n="1" /></td><td>Xianghualing River (Hunan, China; 26.8)<br>Daotian River (Guizhou, China; 99.21)</td><td>Streamflow simulation<br>Simulation of flow (including IGF) and water budget</td></tr><tr><td>Conceptual non-linear one-reservoir groundwater model</td><td>Wang and Brubaker [169]<fnref n="1" /></td><td>Shenandoah River of the Potomac River Basin (USA; 7607)</td><td>Streamflow simulation</td></tr><tr><td>Modified crack flow module; conceptual linear one-reservoir groundwater model</td><td>Eini et al. [36]<fnref n="1" /><br>Amin et al. [170]<fnref n="1" /></td><td>Maharlu Lake (Province of Fars, Iran; 4270)<br>Spring Creek (Pennsylvania, USA; 370)</td><td>Simulation of crack/preferential flow, discharge, and water budget—SWAT-ML and SWAT-CF<br>Simulation of streamflow, nutrient loads, and sediment loads for different agricultural management practices—Topo-SWAT</td></tr><tr><td>Variable source area hydrology; conceptual linear one-reservoir groundwater model</td><td>Amin et al. [171]³<br>Amin et al. [172]³<br>Gunn et al. [173]³</td><td>Spring Creek (Pennsylvania, USA; 370)<br>Spring Creek (Pennsylvania, USA; 370)<br>Spring Creek (Pennsylvania, USA; 370)</td><td>Impact of dairy cropping practices on nutrient and sediment loads<br>Impact of agricultural best management practices on nutrient and sediment loads<br>Impact of climate change with increasing atmospheric CO₂ on watershed hydrology—SWAT-VSA_CO₂ and SWAT-VSA_CO₂+Plant</td></tr><tr><td>SWAT + Water Accounting Plus (WA+) framework</td><td>Delavar et al. [174]<br>Delavar et al. [175]</td><td>Tashk-Bakhtegan (Iran; 27,520)<br>Karkheh River (Iran; 42,267)</td><td>Assessment of water consumption and supply trends under different water management strategies—SWAT-FARS<br>Assessment of water supply and demand conditions in wet and dry periods, based on the water resources, consumption, and withdrawal indicators of the WA+ framework—kSWAT-Karkheh</td></tr></tbody></table>

<footnote n="1">¹These studies reported applications of both a standard SWAT model and a modified SWAT model. ² These studies used the Karst-SWAT version of SWAT developed by Nikolaidis et al. [157] without making any additional modifications to the model. ³ These studies used the Topo-SWAT version of SWAT developed by Amin et al. [170].</footnote>
<table><caption>Table 2. *Cont.*</caption><thead><tr><th>Groundwater Modeling Approach</th><th>Reference</th><th>Basin Name (Region; Size in km²)</th><th>Application—Modified SWAT Name (When Applicable)</th></tr></thead><tbody><tr><td>Conceptual linear three-reservoir groundwater model</td><td>Wang et al. [167]<sup>1</sup><br>Geng et al. [168]<sup>1</sup></td><td>Xianghualing River (Hunan, China; 26.8)<br>Daotian River (Guizhou, China; 99.21)</td><td>Streamflow simulation<br>Simulation of flow (including IGF) and water budget</td></tr><tr><td>Conceptual non-linear one-reservoir groundwater model</td><td>Wang and Brubaker [169]<sup>1</sup></td><td>Shenandoah River of the Potomac River Basin (USA; 7607)</td><td>Streamflow simulation</td></tr><tr><td>Modified crack flow module; conceptual linear one-reservoir groundwater model</td><td>Eini et al. [36]<sup>1</sup><br>Amin et al. [170]<sup>1</sup></td><td>Maharlu Lake (Province of Fars, Iran; 4270)<br>Spring Creek (Pennsylvania, USA; 370)</td><td>Simulation of crack/preferential flow, discharge, and water budget—SWAT-ML and SWAT-CF<br>Simulation of streamflow, nutrient loads, and sediment loads for different agricultural management practices—Topo-SWAT</td></tr><tr><td>Variable source area hydrology; conceptual linear one-reservoir groundwater model</td><td>Amin et al. [171]<sup>3</sup><br>Amin et al. [172]<sup>3</sup><br>Gunn et al. [173]<sup>3</sup></td><td>Spring Creek (Pennsylvania, USA; 370)<br>Spring Creek (Pennsylvania, USA; 370)<br>Spring Creek (Pennsylvania, USA; 370)</td><td>Impact of dairy cropping practices on nutrient and sediment loads<br>Impact of agricultural best management practices on nutrient and sediment loads<br>Impact of climate change with increasing atmospheric CO₂ on watershed hydrology—SWAT-VSA_CO₂ and SWAT-VSA_CO₂+Plant</td></tr><tr><td>SWAT + Water Accounting Plus (WA+) framework</td><td>Delavar et al. [174]<br>Delavar et al. [175]</td><td>Tashk-Bakhtegan (Iran; 27,520)<br>Karkheh River (Iran; 42,267)</td><td>Assessment of water consumption and supply trends under different water management strategies—SWAT-FARS<br>Assessment of water supply and demand conditions in wet and dry periods, based on the water resources, consumption, and withdrawal indicators of the WA+ framework—kSWAT-Karkheh</td></tr></tbody></table>

<footnote n="1">¹These studies reported applications of both a standard SWAT model and a modified SWAT model. ² These studies used the Karst-SWAT version of SWAT developed by Nikolaidis et al. [157] without making any additional modifications to the model. ³ These studies used the Topo-SWAT version of SWAT developed by Amin et al. [170].</footnote>
# 4. Results and Discussion

## 4.1. Applications of Standard SWAT in Karst Watersheds

### 4.1.1. Research Areas Covered in Standard SWAT Applications

The studies that reported an application of standard SWAT in a karst watershed (Table 1) were mainly conducted in the USA (18 articles; 36%), followed by Europe (17 articles; 34%), and Asia (13 articles; 26%). In contrast, only two studies (4%) were identified in Africa.

Different versions of SWAT have been developed over the years to meet the growing need for water resources modeling and management tools, the latest being SWAT+. SWAT+ is a completely restructured version of SWAT that offers an enhanced flexibility in watershed configuration and spatial representation of landscape processes [176,177]. The identified studies in this review were conducted using the previous SWAT versions, including SWAT v2000, v2005, v2009, and v2012. Noticeably, SWAT+ has not yet been implemented in karst regions.

The standard SWAT model has been applied to a wide range of karst dominated and karst influenced watershed scales to assess the hydrological cycle and simulate streamflow [104,105,107,132,150], flood events [123], erosion processes and sediment yield [124,125,133], as well as pollutant (nutrients and pathogens) transport [108,109,127]. The model was also used to compare water quality impacts between scenarios of different crop types and agricultural management practices [110,116].

Several studies evaluated climate change impacts on watershed hydrology based on historical climate patterns and climate projections [113,118,129,135], as well as the effects of land-use change on the water budget [111]. Other studies assessed the combined impacts of land-use and climatic changes on watershed hydrology and or water quality [62,139], including the influence of future urbanization and impervious surface growth [114], and other anthropogenic factors, such as wastewater treatment [5].

In other applications, SWAT was used to simulate the spatial and temporal evolution of runoff, groundwater, erosivity, and surface sediment yield in karst watersheds, considering various climatic and land features. These studies identified the driving factors affecting the variation of ecosystem services and analyzed the trade-offs and synergic relationships between them for rocky desertification containment and ecological protection [138,140,143–145,149].

Additionally, SWAT has been coupled with other models to expand the assessment of flow and water quality. For instance, Sarkar et al. [115] linked SWAT with the Hydrological Simulation Program-FORTRAN (HSPF) to simulate flow and sediment loading from upland agricultural areas in a karstified watershed using SWAT, followed by in-stream sediment processes in HSPF. Sullivan et al. [117] applied SWAT to model recharge nitrate concentrations from natural and anthropogenic sources in a karst watershed. Then, the recharge output from SWAT was incorporated into the Modular Three-Dimensional Finite-Difference Groundwater Flow Model Conduit Flow Process version 2 (MODFLOW CFPv2) and the Conduit Modular 3-Dimensional Transport (CMT3D) model to predict groundwater flow and nitrate transport and levels in the aquifer. Similarly, Karki et al. [120] estimated groundwater recharge in a karst watershed using SWAT, then integrated the recharge output from SWAT into a MODFLOW model with Newton-Raphson formulation (MODFLOW-NWT) to evaluate the impacts of irrigation withdrawals on groundwater levels and the stream-aquifer fluxes. Al Aamery et al. [119] also simulated surface runoff, surface routing, and soil water percolation in SWAT as inputs for a combined discrete-continuum fluviokarst numerical model.

Moreover, the performance of SWAT for karst watersheds hydrological and water quality simulations was evaluated under different precipitation input data [130,142,147] and with respect to various calibration approaches, such as multi-site calibration [148] and zonal calibration that incorporates the basin geological properties [136].
# 4. Results and Discussion

## 4.1. Applications of Standard SWAT in Karst Watersheds

### 4.1.1. Research Areas Covered in Standard SWAT Applications

The studies that reported an application of standard SWAT in a karst watershed (Table 1) were mainly conducted in the USA (18 articles; 36%), followed by Europe (17 articles; 34%), and Asia (13 articles; 26%). In contrast, only two studies (4%) were identified in Africa.

Different versions of SWAT have been developed over the years to meet the growing need for water resources modeling and management tools, the latest being SWAT+. SWAT+ is a completely restructured version of SWAT that offers an enhanced flexibility in watershed configuration and spatial representation of landscape processes [176,177]. The identified studies in this review were conducted using the previous SWAT versions, including SWAT v2000, v2005, v2009, and v2012. Noticeably, SWAT+ has not yet been implemented in karst regions.

The standard SWAT model has been applied to a wide range of karst dominated and karst influenced watershed scales to assess the hydrological cycle and simulate streamflow [104,105,107,132,150], flood events [123], erosion processes and sediment yield [124,125,133], as well as pollutant (nutrients and pathogens) transport [108,109,127]. The model was also used to compare water quality impacts between scenarios of different crop types and agricultural management practices [110,116].

Several studies evaluated climate change impacts on watershed hydrology based on historical climate patterns and climate projections [113,118,129,135], as well as the effects of land-use change on the water budget [111]. Other studies assessed the combined impacts of land-use and climatic changes on watershed hydrology and or water quality [62,139], including the influence of future urbanization and impervious surface growth [114], and other anthropogenic factors, such as wastewater treatment [5].

In other applications, SWAT was used to simulate the spatial and temporal evolution of runoff, groundwater, erosivity, and surface sediment yield in karst watersheds, considering various climatic and land features. These studies identified the driving factors affecting the variation of ecosystem services and analyzed the trade-offs and synergic relationships between them for rocky desertification containment and ecological protection [138,140,143–145,149].

Additionally, SWAT has been coupled with other models to expand the assessment of flow and water quality. For instance, Sarkar et al. [115] linked SWAT with the Hydrological Simulation Program-FORTRAN (HSPF) to simulate flow and sediment loading from upland agricultural areas in a karstified watershed using SWAT, followed by in-stream sediment processes in HSPF. Sullivan et al. [117] applied SWAT to model recharge nitrate concentrations from natural and anthropogenic sources in a karst watershed. Then, the recharge output from SWAT was incorporated into the Modular Three-Dimensional Finite-Difference Groundwater Flow Model Conduit Flow Process version 2 (MODFLOW CFPv2) and the Conduit Modular 3-Dimensional Transport (CMT3D) model to predict groundwater flow and nitrate transport and levels in the aquifer. Similarly, Karki et al. [120] estimated groundwater recharge in a karst watershed using SWAT, then integrated the recharge output from SWAT into a MODFLOW model with Newton-Raphson formulation (MODFLOW-NWT) to evaluate the impacts of irrigation withdrawals on groundwater levels and the stream-aquifer fluxes. Al Aamery et al. [119] also simulated surface runoff, surface routing, and soil water percolation in SWAT as inputs for a combined discrete-continuum fluviokarst numerical model.

Moreover, the performance of SWAT for karst watersheds hydrological and water quality simulations was evaluated under different precipitation input data [130,142,147] and with respect to various calibration approaches, such as multi-site calibration [148] and zonal calibration that incorporates the basin geological properties [136].
4.1.2. Overall Performance of Standard SWAT Models

More than 70% of the studies that used NSE to assess the performance of SWAT hydrological models with daily time series calibration scored NSE values greater than 0.5, and over 90% reported NSE values greater than 0.5 for the daily time series validation. In comparison, more than 90% of the studies scored NSE values higher than 0.5 with monthly calibrated models, while upwards of 80% reported NSE values higher than 0.5 with monthly validation (Figure 2). These results indicate a satisfactory performance, with numerous applications meeting the criteria of a “good” flow simulation, as proposed by Moriasi et al. [102].

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Number of standard SWAT-based studies in karst watersheds under the NSE performance ratings recommended by Moriasi et al. [102] for daily and monthly discharge simulation.</figcaption></figure>

However, some applications conducted in complex karst watersheds scored poor NSE statistics. The studies conducted by Spruill et al. [104] and Coffey et al. [105] in the small experimental watershed of Kentucky revealed that SWAT failed to accurately reproduce peak and low flows. The observed and simulated daily hydrographs were asynchronous, with SWAT often underestimating the peak discharge rates and generating recessions that are faster than the observed data curves. The monthly runoff volumes at the watershed outlet were also underpredicted, which was attributed to the lack of explicit representation of karst geology in SWAT. A similar finding was reached by Benham et al. [106] who concluded that SWAT inability to reproduce the flows sustained by karst features reduced the prediction efficiency of streamflow in their study watershed. At a larger scale, the studies undertaken in the Scandinavian and Iberian peninsulas of 10⁶ km² and 556,000 km² [128], respectively, and in the Danube River basin of 800,000 km² [131] revealed that the performance of SWAT was lower in karst dominated regions in comparison to non-karst areas, due to the model misrepresentation of baseflow in karst streams. Martinez-Salvador and Conesa-Garcia [133] also emphasized on the need to improve the representation of extreme hydrological events (e.g., low flow and peak flow periods) in SWAT.

Furthermore, several studies underlined the need to account for external and interbasin groundwater flows to improve the discharge simulation in SWAT [104,121,124,127,134,137,141]. The hydrological simulations performed by Spruill et al. [104] confirmed the dye tracing results from sinkholes surrounding the study site that an area larger than the watershed topographic boundaries contributes to streamflow. Amatya et al. [107] underlined the need to couple SWAT with a subsurface hydrology model to accurately characterize the dynamics of the karst groundwater flow contribution to the surface drainage network. Gamvroudis et al. [127] estimated that around 33% of the water balance was lost via deep
4.1.2. Overall Performance of Standard SWAT Models

More than 70% of the studies that used NSE to assess the performance of SWAT hydrological models with daily time series calibration scored NSE values greater than 0.5, and over 90% reported NSE values greater than 0.5 for the daily time series validation. In comparison, more than 90% of the studies scored NSE values higher than 0.5 with monthly calibrated models, while upwards of 80% reported NSE values higher than 0.5 with monthly validation (Figure 2). These results indicate a satisfactory performance, with numerous applications meeting the criteria of a “good” flow simulation, as proposed by Moriasi et al. [102].

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Number of standard SWAT-based studies in karst watersheds under the NSE performance ratings recommended by Moriasi et al. [102] for daily and monthly discharge simulation.</figcaption></figure>

However, some applications conducted in complex karst watersheds scored poor NSE statistics. The studies conducted by Spruill et al. [104] and Coffey et al. [105] in the small experimental watershed of Kentucky revealed that SWAT failed to accurately reproduce peak and low flows. The observed and simulated daily hydrographs were asynchronous, with SWAT often underestimating the peak discharge rates and generating recessions that are faster than the observed data curves. The monthly runoff volumes at the watershed outlet were also underpredicted, which was attributed to the lack of explicit representation of karst geology in SWAT. A similar finding was reached by Benham et al. [106] who concluded that SWAT inability to reproduce the flows sustained by karst features reduced the prediction efficiency of streamflow in their study watershed. At a larger scale, the studies undertaken in the Scandinavian and Iberian peninsulas of 10⁶ km² and 556,000 km² [128], respectively, and in the Danube River basin of 800,000 km² [131] revealed that the performance of SWAT was lower in karst dominated regions in comparison to non-karst areas, due to the model misrepresentation of baseflow in karst streams. Martinez-Salvador and Conesa-Garcia [133] also emphasized on the need to improve the representation of extreme hydrological events (e.g., low flow and peak flow periods) in SWAT.

Furthermore, several studies underlined the need to account for external and interbasin groundwater flows to improve the discharge simulation in SWAT [104,121,124,127,134,137,141]. The hydrological simulations performed by Spruill et al. [104] confirmed the dye tracing results from sinkholes surrounding the study site that an area larger than the watershed topographic boundaries contributes to streamflow. Amatya et al. [107] underlined the need to couple SWAT with a subsurface hydrology model to accurately characterize the dynamics of the karst groundwater flow contribution to the surface drainage network. Gamvroudis et al. [127] estimated that around 33% of the water balance was lost via deep
groundwater flow to areas outside their study watershed due to karst formations, while Palazón and Navas [124] simulated the discharge losses by underground flow through swallow holes in the upper part of the study basin. On the other hand, Jakada and Chen [141] confirmed the absence of runoff losses by subterranean flow diversion from their study watershed prior to conducting a hydrological simulation in SWAT. Their finding was based on the results of tracer tests conducted through the sinkholes in the watershed and monitoring of the springs within and outside the basin.

In more complex applications, Salerno and Tartari [121] coupled wavelet analysis with hydrological modeling in SWAT to identify the streamflow components in a non-conservative karst subbasin. After excluding the possibility of an incorrect assessment of the precipitation data, streamflow measurements, and evapotranspiration estimates, a series of continuous wavelet transform, cross wavelet transform, wavelet coherence, and phase difference analyses were applied to precipitation, groundwater levels, observed streamflow, and the time series constructed by the difference between the observed daily discharge and the streamflow simulated by a calibrated SWAT model of the study site. Based on the ensemble of correlations, it was established that the external water contribution to the river discharge was primarily due to groundwater seepage from a hydrogeological catchment that is larger than the surface watershed. The daily time series of the external water contribution was generated by multiplying the SWAT-simulated groundwater inflow by a yearly coefficient. This coefficient was adjusted to match the external contribution time series with the groundwater fluctuations simulated by SWAT and have the annual simulated flows equal to the observed flows. The additional water component improved the prediction efficiency of daily streamflow at the watershed outlet, with NSE increasing from 0.61–0.56 in the calibration and validation periods to 0.66–0.62, and $R^2$ increasing from 0.71–0.69 to 0.74–0.72. The mean absolute error of streamflow underestimation was also reduced from 47 to 33%.

Jian et al. [137] simulated discharge in a non-conservative karst watershed with an initial average discrepancy of 47% between the observed and measured water balances. After ruling out the possibility of invalid precipitation, evapotranspiration, and discharge measurements, the external contribution of the underground flow to streamflow was added as a point source discharge in SWAT, adopting the mean value of the difference in the annual water budget. The hydrological calibration and validation were carried out in a two-stage process. In the first step, the SWAT model and external flow value were calibrated using discharge data, while surface runoff, baseflow, and evapotranspiration were calibrated in the next step using available observational data. As a result, the baseflow component (excluding the external flow contribution) was calibrated in SWAT, and the inclusion of IGF reduced the underestimation bias of streamflow from nearly 50% to less than 3% at the monthly scale and 15% at the daily scale. NSE and $R^2$ values greater than 0.5 and 0.65, respectively, were also reached both in the calibration and validation periods.

More recently, Senent-Aparicio et al. [134] applied SWAT with the atmospheric Chloride Mass Balance (CMB) method to simulate streamflow of the Castil River basin (Spain). The study site is steep karst watershed fed by IGF from adjacent aquifers under steady conditions (i.e., no groundwater abstraction, evapotranspiration from shallow aquifers, or underflow to deep aquifers). The net aquifer discharge was equated to the baseflow component of streamflow, and the CMB approach was used to estimate the fraction of net aquifer recharge from the upstream areas as a proxy for the IGF contributing to additional baseflow. The corrected baseflow time series with IGF improved the SWAT model performance, reducing the underestimation bias of the streamflow simulations to less than 20% in both calibration and validation.

## 4.2. Applications of Modified SWAT in Karst Watersheds

### 4.2.1. Conceptual Linear One-Reservoir Model

This subsection describes the SWAT models with modified recharge functions and a linear one-reservoir groundwater module to simulate karst flows. These models include:
groundwater flow to areas outside their study watershed due to karst formations, while Palazón and Navas [124] simulated the discharge losses by underground flow through swallow holes in the upper part of the study basin. On the other hand, Jakada and Chen [141] confirmed the absence of runoff losses by subterranean flow diversion from their study watershed prior to conducting a hydrological simulation in SWAT. Their finding was based on the results of tracer tests conducted through the sinkholes in the watershed and monitoring of the springs within and outside the basin.

In more complex applications, Salerno and Tartari [121] coupled wavelet analysis with hydrological modeling in SWAT to identify the streamflow components in a non-conservative karst subbasin. After excluding the possibility of an incorrect assessment of the precipitation data, streamflow measurements, and evapotranspiration estimates, a series of continuous wavelet transform, cross wavelet transform, wavelet coherence, and phase difference analyses were applied to precipitation, groundwater levels, observed streamflow, and the time series constructed by the difference between the observed daily discharge and the streamflow simulated by a calibrated SWAT model of the study site. Based on the ensemble of correlations, it was established that the external water contribution to the river discharge was primarily due to groundwater seepage from a hydrogeological catchment that is larger than the surface watershed. The daily time series of the external water contribution was generated by multiplying the SWAT-simulated groundwater inflow by a yearly coefficient. This coefficient was adjusted to match the external contribution time series with the groundwater fluctuations simulated by SWAT and have the annual simulated flows equal to the observed flows. The additional water component improved the prediction efficiency of daily streamflow at the watershed outlet, with NSE increasing from 0.61–0.56 in the calibration and validation periods to 0.66–0.62, and $R^2$ increasing from 0.71–0.69 to 0.74–0.72. The mean absolute error of streamflow underestimation was also reduced from 47 to 33%.

Jian et al. [137] simulated discharge in a non-conservative karst watershed with an initial average discrepancy of 47% between the observed and measured water balances. After ruling out the possibility of invalid precipitation, evapotranspiration, and discharge measurements, the external contribution of the underground flow to streamflow was added as a point source discharge in SWAT, adopting the mean value of the difference in the annual water budget. The hydrological calibration and validation were carried out in a two-stage process. In the first step, the SWAT model and external flow value were calibrated using discharge data, while surface runoff, baseflow, and evapotranspiration were calibrated in the next step using available observational data. As a result, the baseflow component (excluding the external flow contribution) was calibrated in SWAT, and the inclusion of IGF reduced the underestimation bias of streamflow from nearly 50% to less than 3% at the monthly scale and 15% at the daily scale. NSE and $R^2$ values greater than 0.5 and 0.65, respectively, were also reached both in the calibration and validation periods.

More recently, Senent-Aparicio et al. [134] applied SWAT with the atmospheric Chloride Mass Balance (CMB) method to simulate streamflow of the Castil River basin (Spain). The study site is steep karst watershed fed by IGF from adjacent aquifers under steady conditions (i.e., no groundwater abstraction, evapotranspiration from shallow aquifers, or underflow to deep aquifers). The net aquifer discharge was equated to the baseflow component of streamflow, and the CMB approach was used to estimate the fraction of net aquifer recharge from the upstream areas as a proxy for the IGF contributing to additional baseflow. The corrected baseflow time series with IGF improved the SWAT model performance, reducing the underestimation bias of the streamflow simulations to less than 20% in both calibration and validation.

## 4.2. Applications of Modified SWAT in Karst Watersheds

### 4.2.1. Conceptual Linear One-Reservoir Model

This subsection describes the SWAT models with modified recharge functions and a linear one-reservoir groundwater module to simulate karst flows. These models include:
the modified SWAT by [152], SWAT-B&B [153], SWAT-karst [154], KarstSWAT [155], and the modified SWAT by [156]. The latter four models can also represent karst watersheds dominated by flow through sinkholes.

Referring to Table 2, the first application of a modified SWAT code in a karst watershed was performed by Afinowicz et al. [152] to evaluate the impacts of woody plants management scenarios on the rangeland water cycle of the North Fork of the Upper Guadalupe River, Texas (USA). The watershed has an area of 360 km² and is covered by thin soils that overlie fractured limestone formations.

The return flow (baseflow) function of the groundwater module of SWAT (v2000) was modified to simulate rapid infiltration in karst areas into the deep aquifer. Therefore, the deep aquifer recharge component was deducted from the baseflow component of streamflow to allow a fraction of infiltrated water to bypass the shallow aquifer and enter the deep aquifer instead of flowing into the channel as baseflow, as shown in Equation (A1) of Appendix B.

The hydrological model was adjusted using daily streamflow data at the watershed outlet, with a 5-year warm-up period, a 5-year calibration period, and a 7-year validation period. The model scored monthly NSE values of 0.29 and 0.5 for the calibration and validation periods, respectively. It performed less efficiently at the daily scale, with NSE values of 0.4 and 0.09. It also failed to accurately reproduce all discharge trends at the daily scale, particularly high peak flows. The results of the hydrograph simulations were attributed to the nature of the surface runoff in the watershed, which is characterized by sustained low baseflow and very high flow that brings the soil water capacity to saturation.

Baffaut and Benson [153] modified the groundwater recharge equation of SWAT (v2005) to model fast infiltration from sinkholes and losing streams to the aquifer and groundwater flow contribution to surface water. The improved SWAT, known as SWAT-B&B/Adapted SWAT model, was applied to the 3600 km² James River basin in southwest Missouri (USA), characterized by losing streams, sinkholes, and springs.

In SWAT-B&B, recharge into the aquifer was partitioned to two components: (1) the infiltration from the soil bottom, representing slow flow to the porous matrix, and (2) the recharge from sinkholes and losing streams, representing fast flow to the conduits. Sinkholes in the study basin were modeled as ponds with a small drainage area and high hydraulic conductivity, while losing streams were represented by tributary channels with high streambed hydraulic conductivity. Thus, the soil and karst infiltration components were simulated using two recharge functions, each with a specific groundwater delay coefficient (see Equations (A2) and (A3) of Appendix C). Return flow was then modeled with the standard SWAT function, based on the groundwater flow of the previous day and the total aquifer recharge of that day (see Equation (A4) of Appendix C).

The hydrological model was calibrated for 8 years of daily streamflow records at 5 gauging stations and validated for 7 years. Streamflow biases were all less than 25%, ranging between 4% to 20% during the calibration period and -2% to -21% during validation. The percent of bias in surface runoff simulation were all around 10%, indicating a better representation of the baseflow component to streamflow. Moreover, NSE values of around 0.5 were reached for the calibration and validation periods in the main stem of the stream and at the outlet, but lower values close to 0.3 were obtained in the upstream small tributaries. Although a significant improvement in the NSE values could not be spotted by comparing both the standard and modified SWAT models, SWAT-B&B sustained more flows during the dry periods in comparison to SWAT.

The model was then used to estimate in-stream phosphorus loads and concentrations, and fecal coliform concentrations. Poor water quality simulation results were obtained in almost all observational river reaches of the basin, both in calibration and validation periods.

Yactayo [154] further modified the SWAT-B&B code to simulate fast aquifer recharge through sinkholes at the HRU scale by introducing a new parameter called *sink* to the HRU groundwater input file. This sinkhole partitioning coefficient represented the fraction of the runoff drained by a sinkhole to the unconfined aquifer. With this approach, a fraction
the modified SWAT by [152], SWAT-B&B [153], SWAT-karst [154], KarstSWAT [155], and the modified SWAT by [156]. The latter four models can also represent karst watersheds dominated by flow through sinkholes.

Referring to Table 2, the first application of a modified SWAT code in a karst watershed was performed by Afinowicz et al. [152] to evaluate the impacts of woody plants management scenarios on the rangeland water cycle of the North Fork of the Upper Guadalupe River, Texas (USA). The watershed has an area of 360 km² and is covered by thin soils that overlie fractured limestone formations.

The return flow (baseflow) function of the groundwater module of SWAT (v2000) was modified to simulate rapid infiltration in karst areas into the deep aquifer. Therefore, the deep aquifer recharge component was deducted from the baseflow component of streamflow to allow a fraction of infiltrated water to bypass the shallow aquifer and enter the deep aquifer instead of flowing into the channel as baseflow, as shown in Equation (A1) of Appendix B.

The hydrological model was adjusted using daily streamflow data at the watershed outlet, with a 5-year warm-up period, a 5-year calibration period, and a 7-year validation period. The model scored monthly NSE values of 0.29 and 0.5 for the calibration and validation periods, respectively. It performed less efficiently at the daily scale, with NSE values of 0.4 and 0.09. It also failed to accurately reproduce all discharge trends at the daily scale, particularly high peak flows. The results of the hydrograph simulations were attributed to the nature of the surface runoff in the watershed, which is characterized by sustained low baseflow and very high flow that brings the soil water capacity to saturation.

Baffaut and Benson [153] modified the groundwater recharge equation of SWAT (v2005) to model fast infiltration from sinkholes and losing streams to the aquifer and groundwater flow contribution to surface water. The improved SWAT, known as SWAT-B&B/Adapted SWAT model, was applied to the 3600 km² James River basin in southwest Missouri (USA), characterized by losing streams, sinkholes, and springs.

In SWAT-B&B, recharge into the aquifer was partitioned to two components: (1) the infiltration from the soil bottom, representing slow flow to the porous matrix, and (2) the recharge from sinkholes and losing streams, representing fast flow to the conduits. Sinkholes in the study basin were modeled as ponds with a small drainage area and high hydraulic conductivity, while losing streams were represented by tributary channels with high streambed hydraulic conductivity. Thus, the soil and karst infiltration components were simulated using two recharge functions, each with a specific groundwater delay coefficient (see Equations (A2) and (A3) of Appendix C). Return flow was then modeled with the standard SWAT function, based on the groundwater flow of the previous day and the total aquifer recharge of that day (see Equation (A4) of Appendix C).

The hydrological model was calibrated for 8 years of daily streamflow records at 5 gauging stations and validated for 7 years. Streamflow biases were all less than 25%, ranging between 4% to 20% during the calibration period and -2% to -21% during validation. The percent of bias in surface runoff simulation were all around 10%, indicating a better representation of the baseflow component to streamflow. Moreover, NSE values of around 0.5 were reached for the calibration and validation periods in the main stem of the stream and at the outlet, but lower values close to 0.3 were obtained in the upstream small tributaries. Although a significant improvement in the NSE values could not be spotted by comparing both the standard and modified SWAT models, SWAT-B&B sustained more flows during the dry periods in comparison to SWAT.

The model was then used to estimate in-stream phosphorus loads and concentrations, and fecal coliform concentrations. Poor water quality simulation results were obtained in almost all observational river reaches of the basin, both in calibration and validation periods.

Yactayo [154] further modified the SWAT-B&B code to simulate fast aquifer recharge through sinkholes at the HRU scale by introducing a new parameter called *sink* to the HRU groundwater input file. This sinkhole partitioning coefficient represented the fraction of the runoff drained by a sinkhole to the unconfined aquifer. With this approach, a fraction
of the surface runoff and lateral flow in the karst HRU was no longer included in the calculation of total streamflow in the main channel but allocated to the daily seepage from sinking streams and sinkholes. The transmissions losses from the surface runoff entering the sinkholes were also not simulated. Thus, the unconfined aquifer recharge in non-karst regions was calculated using Equation (A5) of Appendix D, whereas aquifer recharge in karst regions was computed using Equation (A6) of Appendix D.

The modified model known as SWAT-karst was applied in the 890.2 km² Opequon Creek watershed, located in the Potomac and Shenandoah River basin in Virginia. For SWAT-karst, a new land-use category was added to the land-use map so that sinkholes may be represented by HRUs, based on the area of the sinkhole regions and the land use where the sinkholes are located. Similar to SWAT-B&B, sinking streams were represented by tributary channels with high hydraulic conductivities.

SWAT-karst, SWAT-B&B and SWAT were run at the daily time step for a period of 11 years and compared in terms of their performance efficiency in simulating streamflow and other water balance components without any model calibration. All three models overestimated streamflow, and the values of the PBIAS, NSE, and RSR were unsatisfactory at all subbasin outlets and streamflow gages where discharge values were compared. Nonetheless, both SWAT-B&B and SWAT-karst performed better than SWAT in simulating karst discharge, and SWAT-karst had a more significant impact on the distribution of the water balance components, by simulating less runoff and more baseflow in karst regions with sinkholes. The authors noted that aquifer recharge diverted by sinkholes to regions outside the watershed could be a reason behind SWAT-karst overestimating discharge and failing to meet the acceptable performance criteria. However, they maintained that parameter *sink* values could be modified to control the depth of water that recharges the unconfined and confined aquifers (water lost from the watershed, see Equation (A7) of Appendix D). Yactayo [154] also modeled the nitrate loading that recharges the aquifers through the sinkhole as a function of: (1) the volume of surface runoff and lateral flow lost to sinkholes in karst regions, and (2) the nitrate aquifer recharge loading from the soil water percolation. Similar to the flow simulation results, the values of the in-stream nitrate concentrations calculated from aquifer recharge and nitrate in baseflow were unsatisfactory.

On the other hand, Palanisamy and Workman [155] incorporated an orifice flow transfer function and a successive summation routing algorithm (SSRA) into SWAT in order to simulate groundwater flow from sinkholes located in the streambed to a spring. The modified SWAT code, called KarstSWAT, was applied to the Cane Run watershed of 115.6 km² in Kentucky (USA), where numerous sinkholes found along the river streambed divert surface runoff through an underground conduit to the main watershed spring. The karst aquifers to which sinkholes drain the river flow largely overlap the Cane Run surface watershed, and runoff routing into the sinkholes depends on the incoming streamflow volume, the sinkhole size, and the capacity of the underground conduit.

To represent this unique hydrological setting, sinkholes were conceptualized as orifices and were modeled as outlets of the karst subbasins during watershed delineation in SWAT. The discharge capacity of the sinkholes was simulated using a head-discharge relationship (see Equation (A8) of Appendix E) as a function of a diameter range that corresponds to the size of the sinkholes. The discharge from the sinkholes and infiltration from the soil profile bottom were then added to the deep aquifer reservoir in SWAT, aggregated at HRU level, and transferred to the spring outlet using the SSRA algorithm with a maximum travel time of one day. The number of the subbasin in which the sinkholes are located and the diameter of the sinkholes were specified in an input file called *sink.dat*, while groundwater basins that drain the aquifer water to the spring were defined in a file called *gw_flow.dat*.

KarstSWAT was calibrated for 3 years using daily streamflow measurements at the Cane Run River, and validated for another 3 years using runoff data at the Cane Run River and spring outlet. Compared to the original SWAT model, KarstSWAT showed a better representation of the hydrological cycle in the karst watershed. The average annual surface runoff and recharge to shallow aquifer decreased by 65% and 91%, respectively, while deep
of the surface runoff and lateral flow in the karst HRU was no longer included in the calculation of total streamflow in the main channel but allocated to the daily seepage from sinking streams and sinkholes. The transmissions losses from the surface runoff entering the sinkholes were also not simulated. Thus, the unconfined aquifer recharge in non-karst regions was calculated using Equation (A5) of Appendix D, whereas aquifer recharge in karst regions was computed using Equation (A6) of Appendix D.

The modified model known as SWAT-karst was applied in the 890.2 km² Opequon Creek watershed, located in the Potomac and Shenandoah River basin in Virginia. For SWAT-karst, a new land-use category was added to the land-use map so that sinkholes may be represented by HRUs, based on the area of the sinkhole regions and the land use where the sinkholes are located. Similar to SWAT-B&B, sinking streams were represented by tributary channels with high hydraulic conductivities.

SWAT-karst, SWAT-B&B and SWAT were run at the daily time step for a period of 11 years and compared in terms of their performance efficiency in simulating streamflow and other water balance components without any model calibration. All three models overestimated streamflow, and the values of the PBIAS, NSE, and RSR were unsatisfactory at all subbasin outlets and streamflow gages where discharge values were compared. Nonetheless, both SWAT-B&B and SWAT-karst performed better than SWAT in simulating karst discharge, and SWAT-karst had a more significant impact on the distribution of the water balance components, by simulating less runoff and more baseflow in karst regions with sinkholes. The authors noted that aquifer recharge diverted by sinkholes to regions outside the watershed could be a reason behind SWAT-karst overestimating discharge and failing to meet the acceptable performance criteria. However, they maintained that parameter *sink* values could be modified to control the depth of water that recharges the unconfined and confined aquifers (water lost from the watershed, see Equation (A7) of Appendix D). Yactayo [154] also modeled the nitrate loading that recharges the aquifers through the sinkhole as a function of: (1) the volume of surface runoff and lateral flow lost to sinkholes in karst regions, and (2) the nitrate aquifer recharge loading from the soil water percolation. Similar to the flow simulation results, the values of the in-stream nitrate concentrations calculated from aquifer recharge and nitrate in baseflow were unsatisfactory.

On the other hand, Palanisamy and Workman [155] incorporated an orifice flow transfer function and a successive summation routing algorithm (SSRA) into SWAT in order to simulate groundwater flow from sinkholes located in the streambed to a spring. The modified SWAT code, called KarstSWAT, was applied to the Cane Run watershed of 115.6 km² in Kentucky (USA), where numerous sinkholes found along the river streambed divert surface runoff through an underground conduit to the main watershed spring. The karst aquifers to which sinkholes drain the river flow largely overlap the Cane Run surface watershed, and runoff routing into the sinkholes depends on the incoming streamflow volume, the sinkhole size, and the capacity of the underground conduit.

To represent this unique hydrological setting, sinkholes were conceptualized as orifices and were modeled as outlets of the karst subbasins during watershed delineation in SWAT. The discharge capacity of the sinkholes was simulated using a head-discharge relationship (see Equation (A8) of Appendix E) as a function of a diameter range that corresponds to the size of the sinkholes. The discharge from the sinkholes and infiltration from the soil profile bottom were then added to the deep aquifer reservoir in SWAT, aggregated at HRU level, and transferred to the spring outlet using the SSRA algorithm with a maximum travel time of one day. The number of the subbasin in which the sinkholes are located and the diameter of the sinkholes were specified in an input file called *sink.dat*, while groundwater basins that drain the aquifer water to the spring were defined in a file called *gw_flow.dat*.

KarstSWAT was calibrated for 3 years using daily streamflow measurements at the Cane Run River, and validated for another 3 years using runoff data at the Cane Run River and spring outlet. Compared to the original SWAT model, KarstSWAT showed a better representation of the hydrological cycle in the karst watershed. The average annual surface runoff and recharge to shallow aquifer decreased by 65% and 91%, respectively, while deep
aquifer recharge increased many folds as water was partially diverted through the sinkholes rather than the soil. The cumulative observed and simulated streamflow plots, with and without sinkholes, also demonstrated that KarstSWAT reduced channel flow during low flow and high flow periods. The modified model performance was further assessed against the original SWAT model under a multitude of runoff events during which at least 10 mm of peak rainfall was observed. Results showed that KarstSWAT improved the prediction of the peak flows and baseflow, with average the NSE and $R^2$ values increasing from 0.23 to 0.77 and 0.78 to 0.87, respectively. Moreover, the discrepancy between the observed and simulated spring flow was attributed to the capacity of the orifices to transfer flow, whereby the overestimation of streamflow by KarstSWAT resulted in the underestimation of spring discharge and vice versa. Nonetheless, discharge at the watershed spring was continuously simulated, showing a good agreement with the observed spring hydrographs at different time periods.

Zhou et al. [156] also modified SWAT (v2012) to simulate fast infiltration through karst sinkholes in the upper course of the South Panjiang River, Southwest China. The basin extends over an area of 2762 km² that is mainly covered by limestone and under the influence of a subtropical humid monsoon climate. Due to the karst effect, sinkholes have formed across the watershed subbasins as opposed to only in the streambed (the case of the study the Cane Run watershed [155]).

The authors used the pond module of SWAT to represent the sinkhole processes. While infiltration from the bottom of the soil profile in both karst and non-karst areas is modeled using the same delay time variable in the original SWAT, the recharge function was modified to simulate the rapid recharge of groundwater aquifer in sinkholes. Water leaving the ponds to the aquifer was separated from percolation with a delay time variable specific to pond leakage and set to 1/50 of its original value. Hence, recharge was divided into two components: leakage recharge of the soil profile and rapid recharge of the karst sinkholes (Equation (A9) in Appendix F). The pond module was added at the subbasin scale, and sinkholes were represented by one pond in each subbasin whereby a fraction of the subbasin area drains the surface flow into the pond. A high hydraulic conductivity value was set at the bottom of the ponds in order to maximize infiltration and groundwater recharge.

The SWAT model was adjusted using monthly streamflow data, with 2 years of warm-up, calibration, and validation each. The modified SWAT model improved the streamflow simulations: the values of the NSE and $R^2$ indicators increased from 0.35–0.66 (calibration-validation) and 0.7–0.76, respectively, in the original SWAT to 0.61–0.79 and 0.74–0.83 in the modified SWAT, with a higher prediction accuracy of the peak flow and baseflow at the daily time interval. The use of the pond module, with large hydraulic conductivity values and short recharge durations, also reduced the surface runoff and lateral flow in the subbasins with sinkholes and increased baseflow depth by rapidly diverting the surface water to the shallow aquifer.

### 4.2.2. Conceptual Linear Two-Reservoir Model

This subsection describes the modified SWAT models with a linear two-reservoir groundwater module to simulate flows of the matrix and conduits components in a karst system. These models are Karst-SWAT [157], KSWAT [165], and SWAT_IFG [166].

Nikolaidis et al. [157] interfaced SWAT with a spreadsheet version of the linear two-reservoir model proposed by Kourgialas et al. [178] to simulate discharge and nitrate transport in the Koiliaris River basin (132 km²) in Crete, Greece, under climate change. The modified SWAT model, known as Karst-SWAT, comprises an upper reservoir representative of the fast flow in the conduits and a lower reservoir for the slow flow in the matrix and narrow fractures. The model uses two proportionality coefficients to partition karst recharge between the two compartments and models another flow fraction from the upper to the lower reservoir. The sum of outflows from the matrix and conduit reservoirs forms the total discharge of the karstic area (see Equations (A10) to (A14) in Appendix G).
aquifer recharge increased many folds as water was partially diverted through the sinkholes rather than the soil. The cumulative observed and simulated streamflow plots, with and without sinkholes, also demonstrated that KarstSWAT reduced channel flow during low flow and high flow periods. The modified model performance was further assessed against the original SWAT model under a multitude of runoff events during which at least 10 mm of peak rainfall was observed. Results showed that KarstSWAT improved the prediction of the peak flows and baseflow, with average the NSE and $R^2$ values increasing from 0.23 to 0.77 and 0.78 to 0.87, respectively. Moreover, the discrepancy between the observed and simulated spring flow was attributed to the capacity of the orifices to transfer flow, whereby the overestimation of streamflow by KarstSWAT resulted in the underestimation of spring discharge and vice versa. Nonetheless, discharge at the watershed spring was continuously simulated, showing a good agreement with the observed spring hydrographs at different time periods.

Zhou et al. [156] also modified SWAT (v2012) to simulate fast infiltration through karst sinkholes in the upper course of the South Panjiang River, Southwest China. The basin extends over an area of 2762 km² that is mainly covered by limestone and under the influence of a subtropical humid monsoon climate. Due to the karst effect, sinkholes have formed across the watershed subbasins as opposed to only in the streambed (the case of the study the Cane Run watershed [155]).

The authors used the pond module of SWAT to represent the sinkhole processes. While infiltration from the bottom of the soil profile in both karst and non-karst areas is modeled using the same delay time variable in the original SWAT, the recharge function was modified to simulate the rapid recharge of groundwater aquifer in sinkholes. Water leaving the ponds to the aquifer was separated from percolation with a delay time variable specific to pond leakage and set to 1/50 of its original value. Hence, recharge was divided into two components: leakage recharge of the soil profile and rapid recharge of the karst sinkholes (Equation (A9) in Appendix F). The pond module was added at the subbasin scale, and sinkholes were represented by one pond in each subbasin whereby a fraction of the subbasin area drains the surface flow into the pond. A high hydraulic conductivity value was set at the bottom of the ponds in order to maximize infiltration and groundwater recharge.

The SWAT model was adjusted using monthly streamflow data, with 2 years of warm-up, calibration, and validation each. The modified SWAT model improved the streamflow simulations: the values of the NSE and $R^2$ indicators increased from 0.35–0.66 (calibration-validation) and 0.7–0.76, respectively, in the original SWAT to 0.61–0.79 and 0.74–0.83 in the modified SWAT, with a higher prediction accuracy of the peak flow and baseflow at the daily time interval. The use of the pond module, with large hydraulic conductivity values and short recharge durations, also reduced the surface runoff and lateral flow in the subbasins with sinkholes and increased baseflow depth by rapidly diverting the surface water to the shallow aquifer.

### 4.2.2. Conceptual Linear Two-Reservoir Model

This subsection describes the modified SWAT models with a linear two-reservoir groundwater module to simulate flows of the matrix and conduits components in a karst system. These models are Karst-SWAT [157], KSWAT [165], and SWAT_IFG [166].

Nikolaidis et al. [157] interfaced SWAT with a spreadsheet version of the linear two-reservoir model proposed by Kourgialas et al. [178] to simulate discharge and nitrate transport in the Koiliaris River basin (132 km²) in Crete, Greece, under climate change. The modified SWAT model, known as Karst-SWAT, comprises an upper reservoir representative of the fast flow in the conduits and a lower reservoir for the slow flow in the matrix and narrow fractures. The model uses two proportionality coefficients to partition karst recharge between the two compartments and models another flow fraction from the upper to the lower reservoir. The sum of outflows from the matrix and conduit reservoirs forms the total discharge of the karstic area (see Equations (A10) to (A14) in Appendix G).
In the case of the Koiliaris River basin, the recharge area of the springs contributing to the total watershed flow extends at least 50 km² beyond its boundaries. Nikolaidis et al. [157] first used SWAT to model the surface hydrological processes (precipitation, evapotranspiration, infiltration, runoff) and route the percolated water to the deep groundwater aquifer. The extent of karst areas contributing to the emergence of springs from outside the watershed boundaries was established based on the geologic knowledge of the study site and a mass balance modeling approach. Then, the SWAT-simulated deep groundwater flow in karst areas was assigned to the karstic reservoir model in order to estimate the spring flow contribution to discharge. After calibration of the reservoir model parameters, the resultant karst flow time series were input to SWAT as a point source to simulate the overall watershed runoff.

The parameters of the karst flow reservoir and SWAT models were adjusted using high frequency flow measurements at the watershed outlet, surface runoff measurements at a major tributary of the river, and long-term monthly spring flow records. The overall model prediction efficiency of discharge was satisfactory. At the monthly time step, NSE values of 0.77–0.61, PBIAS of –22.1%; –11.8%, and RSR of 0.62–0.63 were reached during calibration and validation, respectively, whereas NSE of 0.62–0.43, PBIAS of (–22.3%; –11.6%), and RSR of 0.48–0.75 were achieved for the daily runoff simulations.

From the water quality perspective, Nikolaidis et al. [157] incorporated a nitrate mass balance model to the upper and lower reservoirs of the karst flow model, assuming that nitrate is conservative in karst. A karst factor was added to the nitrate mass balance equation of the lower reservoir to account for the extra dilution of the incoming nitrate loads by the permanent karst flow volume below the spring level. After calculating the nitrogen inputs in the watershed and the extended karst recharge area based on the local land-use practices, the hydrological and water quality modeling parameters were adjusted using nitrates grab sample data at a river tributary and groundwater wells, coupled with high frequency nitrate data, grab samples at the watershed outlet, and flow measurements. The simulated nitrate concentrations were adequate compared to the nitrate grab sample measurements.

The impact of climate change on the water budget of the Koiliaris River basin was also predicted up to the year 2050, using three climate change scenarios for a combination of general and regional circulation climate models. The results of the climatic projections suggested that precipitation, evapotranspiration, and runoff could decrease by 17%, 8%, and 22%, respectively, for the time horizon 2030–2050 compared to 2010–2029.

Nerantzaki et al. [158] later adopted the Karst-SWAT flow model by Nikolaidis et al. [157] to first simulate the hydrology and suspended sediment transport in the Koiliaris River basin then predict the impacts of climate change on discharge, soil erosion, and sediment transport. The concentration of suspended sediments in the karstic watershed was calculated using the same mass balance equations and deep karst factor adopted by Nikolaidis et al. [157] for nitrates. Four additional years of simulation were added to the validation period of the model previously calibrated by Nikolaidis et al. [157]. Next, climate change scenarios were run up the year 2090 after adjusting the most sensitive flow and water quality parameters. The results of the discharge simulations were adequate, with daily NSE, PBIAS, and RSR of 0.8, 25.3%, and 0.45, respectively, and monthly NSE, PBIAS, and RSR of 0.83, 23.4%, and 0.41, respectively. The suspended sediments calibration results were less adequate, with daily NSE of 0.7, PBIAS of 57%, and RSR of 0.55, suggesting an overestimation bias.

The results of the climate change scenarios showed that surface runoff and spring flow could decrease by nearly 70 to 77% between the time periods of 2010–2049 and 2050–2090. The erosion rate of the watershed main subbasin and surface sediments export were also expected to drop by 48% and 55%, respectively, whereas sediments emerging from the springs were not substantially affected by climate change.

Following an analysis of climate change impacts in the Crete Island using Karst-SWAT, Demetropoulou et al. [160] proposed a program of measures to improve water governance
In the case of the Koiliaris River basin, the recharge area of the springs contributing to the total watershed flow extends at least 50 km² beyond its boundaries. Nikolaidis et al. [157] first used SWAT to model the surface hydrological processes (precipitation, evapotranspiration, infiltration, runoff) and route the percolated water to the deep groundwater aquifer. The extent of karst areas contributing to the emergence of springs from outside the watershed boundaries was established based on the geologic knowledge of the study site and a mass balance modeling approach. Then, the SWAT-simulated deep groundwater flow in karst areas was assigned to the karstic reservoir model in order to estimate the spring flow contribution to discharge. After calibration of the reservoir model parameters, the resultant karst flow time series were input to SWAT as a point source to simulate the overall watershed runoff.

The parameters of the karst flow reservoir and SWAT models were adjusted using high frequency flow measurements at the watershed outlet, surface runoff measurements at a major tributary of the river, and long-term monthly spring flow records. The overall model prediction efficiency of discharge was satisfactory. At the monthly time step, NSE values of 0.77–0.61, PBIAS of –22.1%; –11.8%, and RSR of 0.62–0.63 were reached during calibration and validation, respectively, whereas NSE of 0.62–0.43, PBIAS of (–22.3%; –11.6%), and RSR of 0.48–0.75 were achieved for the daily runoff simulations.

From the water quality perspective, Nikolaidis et al. [157] incorporated a nitrate mass balance model to the upper and lower reservoirs of the karst flow model, assuming that nitrate is conservative in karst. A karst factor was added to the nitrate mass balance equation of the lower reservoir to account for the extra dilution of the incoming nitrate loads by the permanent karst flow volume below the spring level. After calculating the nitrogen inputs in the watershed and the extended karst recharge area based on the local land-use practices, the hydrological and water quality modeling parameters were adjusted using nitrates grab sample data at a river tributary and groundwater wells, coupled with high frequency nitrate data, grab samples at the watershed outlet, and flow measurements. The simulated nitrate concentrations were adequate compared to the nitrate grab sample measurements.

The impact of climate change on the water budget of the Koiliaris River basin was also predicted up to the year 2050, using three climate change scenarios for a combination of general and regional circulation climate models. The results of the climatic projections suggested that precipitation, evapotranspiration, and runoff could decrease by 17%, 8%, and 22%, respectively, for the time horizon 2030–2050 compared to 2010–2029.

Nerantzaki et al. [158] later adopted the Karst-SWAT flow model by Nikolaidis et al. [157] to first simulate the hydrology and suspended sediment transport in the Koiliaris River basin then predict the impacts of climate change on discharge, soil erosion, and sediment transport. The concentration of suspended sediments in the karstic watershed was calculated using the same mass balance equations and deep karst factor adopted by Nikolaidis et al. [157] for nitrates. Four additional years of simulation were added to the validation period of the model previously calibrated by Nikolaidis et al. [157]. Next, climate change scenarios were run up the year 2090 after adjusting the most sensitive flow and water quality parameters. The results of the discharge simulations were adequate, with daily NSE, PBIAS, and RSR of 0.8, 25.3%, and 0.45, respectively, and monthly NSE, PBIAS, and RSR of 0.83, 23.4%, and 0.41, respectively. The suspended sediments calibration results were less adequate, with daily NSE of 0.7, PBIAS of 57%, and RSR of 0.55, suggesting an overestimation bias.

The results of the climate change scenarios showed that surface runoff and spring flow could decrease by nearly 70 to 77% between the time periods of 2010–2049 and 2050–2090. The erosion rate of the watershed main subbasin and surface sediments export were also expected to drop by 48% and 55%, respectively, whereas sediments emerging from the springs were not substantially affected by climate change.

Following an analysis of climate change impacts in the Crete Island using Karst-SWAT, Demetropoulou et al. [160] proposed a program of measures to improve water governance
in the 525 km² Geropotamos basin located in central-southern part of the island. Nerantzaki et al. [161] also used Karst-SWAT to forecast the hydrological response of the Crete region under climate change scenarios up to the year 2098, considering different irrigation sources in SWAT. Moreover, Tapoglou et al. [159] applied Karst-SWAT to predict the impact of climate change on the hydrological cycle and the frequency of extreme hydrological and meteorological events in Crete. Nerantzaki et al. [163] further expanded the research work conducted in the Koiliaris River basin with Karst-SWAT by assessing: (1) the uncertainty of the watershed runoff and karstic flow simulations due to the parameter uncertainty in SWAT and Karst-SWAT, and (2) the impact of internal variability (or stochastic uncertainty) of the meteorological input data on the flow simulations for the reference period and under the climate change scenarios. The uncertainty of the flow models was estimated by combining the Sequential Uncertainty Fitting Version 2 (SUFI-2) in the SWAT Calibration and Uncertainty Program (SWAT-CUP) interface and the @RISK by PALISADE software, while the effect of input internal variability on the flow output was evaluated using Monte Carlo simulations.

Within the framework of studying the hydrological and geochemical processes in the Koiliaris European Critical Zone Observatory, Lilli et al. [164] used Karst-SWAT to simulate the hydrological budget of the Koiliaris River basin and gain insight on the hydrological pathways and response of the karst during extreme events. Additionally, Karst-SWAT was applied to simulate surface flow and spring flow in the Koiliaris River basin, which were required to determine the design flows and flood frequency within the framework of developing nature-based solutions for the riparian forest restoration and flood protection project at the Koiliaris Critical Zone Observatory [162].

Next, Malagò et al. [165] developed a two-reservoir modeling approach by linking SWAT-B&B [153] and Karst-SWAT [157]. The resultant hybrid model, called KSWAT, was used in conjunction with SWAT to simulate the water balance, spring flow, and total discharge in the Island of Crete in Greece. The study area extends over 8336 km² of which 2730 km² is karst.

A SWAT model of the Crete Island was set up and the modified model KSWAT was applied only in the karst subbasins of the region. The daily aquifer recharge from the karst subbasins was simulated using SWAT-B&B. The area of the subbasins contributing to the recharge of a particular spring or group of springs was identified based on the local geological maps and dominant karst soils. Recharge from the soil profile bottom, stream losses, and seepage from other water bodies to the deep aquifer were maximized by: (1) setting the deep aquifer percolation fraction and minimum groundwater delay to 1, (2) adjusting the groundwater coefficient of capillary rise to 0.1 to prevent the upward movement of water to the unsaturated zone, and (3) minimizing the return flow from the shallow and deep aquifers in SWAT. The deep aquifer recharge time series generated by SWAT-B&B were then input to Karst-SWAT in order to simulate and calibrate the discharge of the springs. The parameters of the Karst-SWAT model were adjusted based on daily spring discharge data from 47 gauging stations in Crete.

The hydrologic model in SWAT was adjusted using a step-wise calibration with monthly streamflow data at 15 stream gauging stations. Snow, surface runoff, lateral flow, and baseflow parameters were first calibrated separately in order to adjust the timing of the runoff signal and the discharge values (peak flow and baseflow). Then, the model was recalibrated based on the streamflow-related parameters combined with the adjusted variables of the other water budget components. The final near optimal parameter set of the calibrated subbasins was transferred to the ungauged subbasins using the hydrological similarity approach with the Partial Least Squares Regression, in order to identify similar subbasins based on the correlation between the watershed and discharge characteristics.

Subsequently, the calibrated spring discharge time series from Karst-SWAT were added to the Crete SWAT model as point sources in order to predict the total monthly runoff across the island, and a final calibration was performed to adjust discharge. The results of the performance indicators showed that only 40% of the calibrated gauging