Modified SWAT Model for Agricultural Watershed in Karst Area of Southwest China

Abstract

The Soil and Water Assessment Tool (SWAT) model is extensively used globally for hydrological and water quality assessments but encounters challenges in karst regions due to their complex surface and groundwater hydrological environments. This study aims to refine the delineation of hydrological response units within the SWAT model by combining geomorphological classification and to enhance the model with an epikarst zone hydrological process module, exploring the accuracy improvement of SWAT model simulations in karst regions of Southwest China. Compared with the simulation results of the original SWAT model, we simulated runoff and nutrient concentrations in the Mudong watershed from January 2017 to December 2021 using the improved SWAT model. The simulation results indicated that the modified SWAT model responded more rapidly to precipitation events, particularly in bare karst landform, aligning more closely with the actual hydrological processes in Southwest China’s karst regions. In terms of the predictive accuracy for monthly loads of total nitrogen (TN) and total phosphorus (TP), the coefficient of determination (R²) value of the modified model increased by 10.3% and 9.7%, respectively, and the Nash–Sutcliffe efficiency coefficient (NSE) increased by 11.3% and 9.9%, respectively. The modified SWAT model improves prediction accuracy in karst areas and holds significant practical value for guiding non-point source pollution control in agricultural watersheds.

1. Introduction

With the growth of the global population and the intensification of agriculture, non-point source pollution from agricultural activities has emerged as a significant contributor to the degradation of water quality in basins [1,2,3,4]. Nitrogen and phosphorus from these activities enter water bodies through surface runoff and subsurface flow, leading to ecological issues, such as algal blooms, hypoxia, and a decline in biodiversity [5,6,7]. Recent studies have highlighted the severity of these issues, particularly in regions with high agricultural activity and limited water resources. This is particularly pronounced in karst regions, where unique geological and hydrological features complicate the transport and transformation of these nutrients [8,9,10], exacerbating their impact on water quality [11,12,13]. Consequently, simulating the transport of agricultural nitrogen and phosphorus in karst watersheds, identifying key hydrological processes influencing their transport, and developing strategies for their effective control and reduction are crucial for safeguarding water resources and preserving the ecological balance.

In the realm of watershed-scale modeling of hydrological cycles and pollutant transport processes, the Soil and Water Assessment Tool (SWAT) is highly regarded for its robust physical mechanisms and broad application prospects [14,15,16]. Comprising three sub-models—hydrologic processes, soil erosion, and pollution load—SWAT is a powerful tool for long-term, distributed watershed simulations [17,18]. It has been extensively applied to simulate and evaluate watershed water quantity and quality across various geographical and environmental conditions [19,20,21] in different regions [22,23,24], achieving satisfactory results. The SWAT model simulates the migration process of agricultural non-point source pollutants as they move along with the hydrological activities within the watershed. It can effectively and quantitatively describe the migration pathways and influencing factors of agricultural non-point source pollutants in the study area, and is widely recognized as an effective tool [25,26,27].

However, the SWAT model does not fully account for karst terrains’ unique underground hydrological characteristics and the intricacies of non-point source pollutant migration. Underground conduits and fracture systems in karst formations facilitate rapid water infiltration and pollutant transport, with special groundwater dynamics that are challenging to replicate accurately within the traditional SWAT framework [28]. Moreover, the heterogeneity and anisotropy inherent in karst landscapes result in high spatial variability in pollutant migration rates and pathways across different media, further complicating the simulation efforts [29,30]. Simulation results of non-point source pollution in karst watersheds based on the SWAT model are not always satisfactory, and more specialized calibration is often required to achieve more accurate parameter settings. Despite advancements in tailoring the SWAT model for karst-specific hydrogeological conditions, these modifications are not universally applicable across diverse karst geomorphic settings [31,32,33]. For instance, Amatya et al. [34] noted that the model struggled to estimate karst spring discharge and effectively simulate nitrogen pollutant transport in the Chapel Branch Creek Karst basin. This highlights the need for SWAT model enhancements that are sensitive to karst hydrology’s nuances. Similarly, Baffaut and Benson [35] achieved promising results in modeling nitrogen transport in karst systems characterized by conduits and sinkholes. However, the generalizability of these improvements to other karst geomorphologies warrants further investigation.

In regions where exposed and covered karst landforms are interspersed, the transformation of precipitation into surface and shallow groundwater flow is particularly complex. The anisotropic transport of non-point source pollutants with water further complicates the exchange and transport of these pollutants between surface and shallow groundwater systems, exceeding the complexity observed in non-karst areas [36].

In this context, the objective was to address the limitations of the SWAT model in simulating the transport of nitrogen and phosphorus pollutants in agricultural areas with mixed exposed and covered karst landforms. By refining the original algorithms that govern the migration of non-point source pollutants with water, the SWAT model’s capability to accurately depict the dynamics of non-point source pollutants in karst water systems within agricultural watershed in Southwest China was enhanced. Taking the Mudong watershed, located in the Huixian karst wetland within the Lijiang River basin, as a case study, a modified version of the SWAT model was constructed and validated. The results of this study will contribute to a better understanding of nutrient transport dynamics in karst watersheds and provide a basis for improving the SWAT model’s performance in similar environments. Through simulation of agricultural non-point source pollution loads in this catchment, it provides a scientific basis for preventing, controlling, and managing non-point source pollution in karst agricultural catchments. It also offers a reference for watershed management under similar hydrogeological conditions.

2. Materials and Methods

2.1. Study Area

The Mudong watershed (110°09′ to 110°14′ E, 25°04′–25°09′ N) is located northwest of Guilin City in Guangxi (Figure 1). It resides in the core area of the Huixian wetland at the watershed divide between the Guilin (Lijiang River) and Liuzhou (Liujiang River) Basins. Huixian Wetland is one of the largest natural karst wetlands in Southwest China, playing a crucial role in regulating the regional ecological environment. In the regional topographic pattern, both the surface water systems and the drainage directions of the main groundwater drainage direction converge towards the central structural karst depression (i.e., the Huixian area) from the north and south ends. Under the continuous and complex interaction of surface and groundwater, numerous lakes and marshes have developed, ultimately forming the Huixian Karst Wetland.

The northern part of the Mudong watershed features a high-lying peak cluster valley. The central part consists of a low-lying flat karst peak cluster depression and solitary peak plain, with an altitude of 125–485 m and a catchment area of approximately 30 km².

The Mudong watershed is characterized by a subtropical monsoon climate, with an annual average temperature of approximately 20 °C and a potential evaporation of 1568 mm. Precipitation is abundant but unevenly distributed, with an average annual rainfall of 1800 mm, of which 75–80% occurs during the rainy season, typically from April to September.

As shown in Figure 2a, the predominant soil type within the watershed is paddy soil, (mean pH value ranging from 5.6 to 6.0 and containing a significant amount of calcium, iron, and manganese concretions). The secondary soil types are calcareous soil (mean pH value greater than 7.0 and an organic matter content between 1.4% and 5.2%), and lateritic soil (mean pH value ranging from 5.7 to 6.1). The steep karst limestone mountains are primarily located in the northern part of the watershed (Figure 2b), corresponding to the typical exposed karst landscape (Figure 2c). This typical geomorphology has shaped the pattern of land use distribution. Flat areas are distributed with arable land, water, grasslands, orchards, and villages, while shrubbery is found on the karst stone mountains (Figure 2d). The watershed is known for its agricultural productivity, primarily rice cultivation, and lack of industrial activity. However, the rural infrastructure remains underdeveloped. The recent expansion of the livestock, poultry, and fishing industries has led to increasingly prominent water environment issues in Huixian soluble rock wetlands. The current water quality is suboptimal, and water environmental protection is facing increasingly severe challenges.

2.2. Improvement of the SWAT Model

2.2.1. Improvement of Hydrological Process

The generalizations of hydrological cycle processes in the interleaved exposed and covered karst areas in Southwest China are shown in Figure 3, with the following characteristics. In areas with covered karst topography, after precipitation forms surface runoff, most infiltrates the loose Quaternary soil layer. After being regulated by the soil’s water-holding medium, one portion returns to the atmosphere, while the other part slowly seeps into the underlying epikarst aquifer. The remaining small portion formed slope runoff, some of which infiltrated and replenished the epikarst zone aquifer in areas with relatively low terrain and high karst development, with a hydrological cycle like that of the exposed karst topography area (Figure 4a). In areas with exposed karst topography, after precipitation forms surface runoff, it mainly infiltrates quickly into the epikarst zone through karst structures with good hydraulic connectivity and high permeability, forming a saturated zone. After being regulated by the epikarst zone [37,38,39], a small portion of the water returns to the surface through evaporation and biological processes, while the rest directly enters the saturated zone through karst structures. Another portion of the water was discharged to the surface as springs in relatively lower areas along the epikarst zone (Figure 4b).

In Southwest China, where exposed and covered karst landscapes are interspersed, the seepage of water through the epikarst zone is prone to encounter impure carbonate rock layers in the central part of the watershed. This geological feature significantly impedes the flow, leading to the emergence of karst groundwater through springs, which serves as a primary source of surface runoff. Additionally, a minor portion of the underground flow directly joins the river channels in the form of underflow. In response to these hydrological characteristics, this study has developed and integrated a surface karst module within the SWAT model framework. An epikarst module was added to the SWAT model to improve its water flow production process. Figure 5 illustrates the computational process of the modified SWAT model. Compared to the SWAT model, the modified SWAT model incorporates a hydrological process structure of the epikarst zone between the soil hydrological process and the shallow groundwater hydrological process, allowing the program to first perform epikarst flow calculations after running soil water flow calculations, and then proceed to the aquifer flow calculation phase. The source code of the SWAT was modified according to the new model calculation structure to accurately reflect the groundwater production characteristics of karst watersheds in Southwest China.

Based on the hydrological response unit, the water balance equation of the newly added epikarst zone was defined as follows:

$$W_{ep,j}=W_{ep,j-1}+W_{seep, j}-W_{epseep,j}-Q_{spr,j}$$

(1)

where $W_{ep,j}$ and $W_{ep,j-1}$ are the water contents in the epikarst zone on days j and (j − 1), respectively. $W_{seep, j}$ is the amount of water entering the epikarst zone through the bottom of the soil profile on day j, which is equal to the daily precipitation $R_{day,j}$ when the vertical distribution of the hydrological response unit consists of only the epikarst zone and bedrock. $W_{epseep,j}$ is the amount of water seeping from the epikarst zone into the shallow aquifer on day j, and $Q_{spr,j}$ is the spring flow discharged from the epikarst zone on day j. The unit for all variables was mm. For a visual representation of these parameters, see Figure 6.

Because the karst spring outflow in the Mudong watershed is caused by differences in geological structure and stratigraphic lithology and has a time-lag effect, the linear seepage law is not applicable. Therefore, the linear reservoir model was used to calculate the amount of spring water $Q_{spr,j}$ discharged in the epikarst zone:

$$Q_{spr,j}=Q_{spr,j-1}{·e}^{-{\alpha }_{spr}∆t}+W_{seep,j}(1-e^{-{\alpha }_{spr}∆t})$$

(2)

where $Q_{spr,j}$ and $Q_{spr,j-1}$ are the spring discharge from the epikarst zone on day j and day (j − 1); $W_{seep,j}$ refers to Formula (1) for details; ${\alpha }_{spr}$ is the attenuation coefficient of the spring flow (days); $∆t$ is the time step in days. The amount of water seeping from the epikarst zone into the shallow aquifer $W_{epseep}$ is calculated as follows:

$$W_{epseep}=W_{cra}+W_{slow}$$

(3)

where $W_{cra}$ and $W_{slow}$ are the fast- and slow-flow water seepage into shallow aquifers through large cracks and small solution gaps of the epikarst bandwidth (mm), respectively.

The $W_{cra}$ of the fast-flow water seepage into a shallow aquifer with large karst bandwidths is defined as follows:

$$W_{cra}={\theta }_{ep}·Z_{ep}$$

(4)

where ${\theta }_{ep}$ is the wide fracture rate and $Z_{ep}$ is the thickness of the epikarst zone (mm).

The slow flow process of the seepage of small solution gaps in the epikarst zone into the shallow aquifer is assumed to be a water column process with a variable head:

$$W_{slow}={\mathrm{W}}_{ep}·\left[1-\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{-∆t}{{TT}_{ep}}}\right)\right]$$

(5)

where $W_{slow}$ is the slow flow water seepage into the shallow aquifer through the small solution gap in the epikarst zone (mm) and $W_{ep}$ is the water content of the epikarst zone (mm). ∆t and ${TT}_{ep}$ are the time step and total water movement time (d), respectively.

The replenishment amount of shallow aquifer after increasing the epikarst zone was calculated by the following:

$$W_{rechrg, j}=W_{rechrg,j-1}{·e}^{-1/{\delta }_{gw}}+W_{epseep,j}(1-e^{-1/{\delta }_{gw}})$$

(6)

where $W_{rechrg, j}$ and $W_{rechrg,j-1}$ are the replenishment amounts of shallow aquifers on days j and (j − 1), respectively (mm); $W_{epseep,j}$ is the amount of water entering the shallow aquifer through the bottom of the epikarst zone on day J (mm); ${\alpha }_{gw}$ is the regression coefficient of the base flow; and ${\delta }_{gw}$ and ∆t are the recharge lag time and time step (d), respectively.

In summary, compared with the SWAT model (Figure 6a), the hydrological process of the land surface stage in the modified SWAT model after adding the epikarst zones (Figure 6b) showed an increase in precipitation infiltration and underground runoff.

2.2.2. Improvement of Non-Point Source Pollution Calculation

The influence of the addition of the epikarst zone on nitrogen and phosphorus nutrients is mainly reflected in the change in dissolved nitrogen load in the lateral flow and lower percolation. In the modified SWAT model, the dissolved nitrogen loads ${P_{ep,NO}}_{3,lat}$ migrating with lateral flow were calculated as follows:

$${P_{ep,NO}}_{3,lat}={{conc}_{NO}}_{3,mobile}·(Q_{lat}+Q_{spr})$$

(7)

Dissolved nitrogen load with downward seepage migration ${P_{ep,NO}}_{3,epseep}$ was defined as the following:

$${P_{ep,NO}}_{3,epseep}={{conc}_{NO}}_{3,mobile}·W_{epseep}$$

(8)

where ${P_{ep,NO}}_{3,lat}$ is the dissolved nitrogen loads of the lateral flow in the soil layer and epikarst zone (kg/hm²), ${{conc}_{NO}}_{3,mobile}$ is the concentration of dissolved nitrogen in free water in the soil layer and epikarst zone (kg/mm), $Q_{lat}$ is the amount of water flowing out of the soil layer through the soil flow (mm), $Q_{spr}$ is the spring discharge through the epikarst zone (mm), ${P_{ep,NO}}_{3,epseep}$ is the dissolved nitrogen load migrating with the middle-lower percolation in the epikarst zone (kg/hm²), and $W_{epseep}$ is the amount of water seepage from the epikarst zone into the shallow aquifer (mm).

Figure 7a shows the calculation of nitrogen and phosphorus moving along the water movement path [40]. In the modified SWAT model, the nitrogen and phosphorus loads varied with the amount of free water in the soil layer and epikarst zone (Figure 7b).

2.2.3. Improvement of the Hydrological Response Unit

In the SWAT model, the basic operation of hydrology and the nitrogen and phosphorus transport module take the Hydrological Response Unit (HRU) as the basic calculation unit [41,42,43]. HRUs were assigned according to the underlying land use/cover, soil physical and chemical properties, and slope. The soil layer of the bare karst geomorphic area in the Mudong watershed is shallow and discontinuous, and an epikarst zone forms on the surface in the bare karst geomorphic area or under the soil layer in the covered karst geomorphic area. To depict this feature, the area with a slope >50° in the exposed karst geomorphic area is generalized as a no soil layer on the surface, the bedrock is overlaid with bare rock of the developing epikarst zone, and the area with a slope ≤50° has a soil layer on the surface. Therefore, before spatial discretization, a new soil-type map (Figure 8) was created by overlapping the original and karst geomorphic distribution maps on the ArcGIS platform (https://www.arcgis.com, accessed on 5 December 2018). The detailed operations are as follows:

(1)

The slope of the entire basin was divided into ≤50° and >50° slopes.

(2)

The exposed and covered karst geomorphic distribution maps were overlapped with the slope distribution maps, and the exposed karst geomorphic areas with slopes > 50° and ≤50° were divided.

(3)

The slope map of the karst landform overlapped with the original soil-type distribution map, redefining the soil type with a slope >50° in the exposed karst landform area as bare rock and reclassifying it to assign a unique soil ID.

Therefore, the vertical distribution of vegetation cover, soil layer, and bedrock was distinct from that of the SWAT model (Figure 9a). The modified SWAT model has two types of hydrological response units: one type includes only the epikarst zone and bedrock, and the other is divided into vegetation cover, soil layer, epikarst zone, and bedrock (Figure 9b) in the vertical direction.

Finally, the SWAT source code is modified using the Fortran language in the Windows platform application development environment Visual Studio 2010 (https://learn.microsoft.com, accessed on 15 March 2019), according to the above structure. The modified subroutines include modparm.f, allocate_parms.f, percmain.f, zero (). f, gwmod.f, and sub-basing. f. After the modifications were successfully compiled, the original SWAT.exe executable file is replaced.

2.3. Model Validity Evaluation

The validity of the model was assessed using the Coefficient of Determination R² and the Nash–Sutcliffe Efficiency coefficient (NSE). The R² value represents the correlation between the simulated and observed values; the closer R² is to 1, the higher the consistency between the simulated and observed values. The NSE characterizes the overall efficiency of the model’s simulation. A NSE value of 1 indicates a perfect match between observed and simulated values, NSE = 0 suggests that the simulated results are equivalent to the mean of the observed values, and NSE < 0 indicates that the simulation is less accurate than the mean observed value, hence considered ineffective. A model is considered suitable for the study area and the simulation results are deemed acceptable when R² > 0.6 and NSE > 0.5. The formulas for calculating R² and NSE are as follows:

$${\mathrm{R}}^2={\left[{\displaystyle \frac{{\sum }_{m=1}^n\left(O_m-\overline{O}\right)\left(P_m-\overline{P}\right)}{{\sum }_{m=1}^n{\left(O_m-\overline{O}\right)}^2{\sum }_{m=1}^n{\left(P_m-\overline{P}\right)}^2}}\right]}^2$$

(9)

$$\mathrm{N}\mathrm{S}\mathrm{E}=1-{\displaystyle \frac{{\sum }_{m=1}^n{\left(O_m-P_m\right)}^2}{{\sum }_{m=1}^n{\left(O_m-\overline{O}\right)}^2}}$$

(10)

where $O_m$ represents the observed value for the m-th time step, $\overline{O}$ denotes the mean of all observed data over the simulation period, $P_m$ represents the simulated value for the m-th time step, and $\overline{P}$ denotes the mean of all simulated values over the simulation period.

3. SWAT Model Setup and Calibration

3.1. SWAT Model Setup

Based on the use of spatial data (30 m resolution global digital elevation model dataset of Geospatial Data Cloud of Computer Network Information Center of Chinese Academy of Sciences, https://www.gscloud.cn, (accessed on 15 October 2020)), soil attribute data (divided into five types of lime soil, paddy soil, red loam, water area, and bare rock), and land use-type data (30 m × 30 m resolution land use type data from the Environmental and Ecological Science Data Center of Western China, http://www.nieer.cas.cn/, (accessed on 3 December 2020)). Based on the spatial distribution of natural river channels in the Mudong watershed, burn-in was applied to import the stream network drawn by Google Earth. The water flow direction was then calculated, and the river network was extracted with a threshold of 60 hm² catchment area to generate a river network distribution that was more consistent with the actual situation of the basin. By manually adding river cross section monitoring points, the SWAT model divides the Mudong watershed into seven sub-basins with a catchment area ranging from 0.9 km² to 12.6 km². Multiple gradients (0–24.9°, 24.9–50°, and >50°) were used to divide the hydrological response units of the basin. The Mudong watershed is divided into 105 hydrological response units, including 44 exposed and 61 covered karst landforms (Figure 8).

Using hydrometeorological data (meteorological data from meteorological stations in Guilin from 2017 to 2019, http://data.cma.cn/data, (accessed on 3 December 2020)) as the driving data of the model, the modified SWAT model was constructed using meteorological data from 2019 to 2021 as validation data for the model, non-point source pollutants (survey statistics of rural domestic pollution, nitrogen and phosphorus pollution of livestock and poultry breeding, and aquaculture), and other attributed data.

3.2. SWAT Model Calibration and Validation

Based on the daily monitoring value of Levelogger, a three-parameter automatic water level meter from the Institute of Karst Geology, Chinese Academy of Geological Sciences [44], the measured data combined with precipitation data were used to modify the monthly average discharge data of the watershed outlet. In the modified SWAT model, we used the measured average monthly discharge data from January 2017 to December 2019, along with the model parameters related to the monthly load rates of total nitrogen and total phosphorus, for calibrating the model’s monthly flow and nutrient load simulations. Then, we selected the average monthly discharge from May 2020 to December 2021, as well as the monthly loads of total nitrogen and total phosphorus, to validate the simulation results of the model. To obtain reasonable initial conditions, the preheating period of the model was set to two years.

The Sequential Uncertainty Fitting Version 2 (SUFI-2) algorithm, which is part of the SWAT-CUP (SWAT Uncertainty Program) software for automatic calibration and uncertainty analysis, was utilized to conduct a sensitivity analysis of the model parameters [45]. Based on hydrological processes [46,47], solute transport characteristics [48,49,50] and modeling experience [51,52,53] of karst watersheds, 40 parameters related to water quantity simulation and nitrogen and phosphorus water quality simulations were first selected in the SWAT model for global sensitivity analysis. After the Global Sensitivity Analysis, One-at-a-time Sensitivity Analysis, and iterative simulation, 15 parameters were selected for calibration based on the characteristics of the karst watershed in the study area (

Table 1. Calibration results of the sensitivity parameter used in the modified SWAT model.

Notation	Parameter	Range	Final Value
CN2.mgt	Moisture condition Il curve number	35 to 98	75.5
ALPHA_BF.gw	αgw: Baseflow recession constant	0 to 1	0.75
GW_DELAY.gw	δgw: Delay time for aquifer recharge (days)	0 to 500	2.8
REVAPMN.gw	Threshold water level in the shallow aquifer for revap (mm)	0 to 600	578
RCHRG_DP.gw	βdeep: Aquifer percolation coefficient	0 to 1	0.015
CH_N2.rte	Manning’s “n” value for the main channel	0.01 to 0.3	0.13
SOL_AWC.sol	AWCly: available water capacity	0 to 1	0.08
GWQMN.gw	Threshold water level in shallow aquifer for base flow(mm)	0 to 5000	4010.9
SHALLST_N.gw	NO3sh: Amount of nitrate in the shallow aquifer (kg N/ha)	0.34 to 0.68	0.5
GWSOLP.gw	Soluble phosphorus concentration in groundwater flow (mg P/L)	0.1 to 0.5	0.26
ERORGN.hru	ƐN,sed: Nitrogens enrichment ratio	0 to 5	3.55
ERORGP.hru	ƐP,sed: Phosphorus enrichment ratio	0 to 5	1.45
NPERCO.bsn	βNO3: Nitrate percolation coefficient	0.01 to 1.0	0.64
PPERCO.bsn	kd,perc: Phosphorus percolation coefficient (m3/Mg)	10 to 17.5	12.97
BIOMIX.mgt	Biological mixing efficiency	0 to 1	0.027

). The coefficient of determination (R²) and Nash–Sutcliffe efficiency coefficient (NSE) were selected to evaluate the efficiency of the model [54].

3.2.1. Flow Rate

Figure 10a of the simulation results of the modified SWAT model shows that the coefficient R² of the determination of the mean monthly runoff in the rate period and verification period is greater than 0.86, and the Nash–Sutcliffe efficiency coefficient (NSE) is greater than 0.76, indicating that the simulated value is close to the measured value and the correlation is good.

In terms of monthly runoff volume prediction accuracy, the modified SWAT model demonstrated an average increase of approximately 6.84% in the coefficient of determination (R²) and about 14.79% in the NSE compared to the original SWAT model. According to Figure 10a, compared with the SWAT model, the average monthly runoff simulated by the modified SWAT model increases at its peak in the rainy season. It decreases at its base in the dry season.

3.2.2. Nutrient Load

The simulation results of the modified SWAT model in Figure 10b and Figure 10c show that the R² and NSE of TN monthly load in the calibration period are 0.77 and 0.70, respectively, while in the validation period, they are 0.73 and 0.68, respectively. The R² and NSE of TP monthly load in the calibration period are 0.68 and 0.67, respectively, while in the validation period, they are 0.67 and 0.66, respectively, indicating that the simulated values of TN and TP monthly load in the basin are well matched with the measured values and basically meet the accuracy requirements of the simulation.

As shown in Figure 10b,c, the basin’s monthly TN and TP load presents a similar change rule under the influence of karst underground aquatic flow. The comparison of simulation results from 2017 to 2021 shows (

Table 2. Annual performance of the modified SWAT model compared with the SWAT model inflow and nutrient load simulations.

	SWAT			Modified SWAT
Year	Flow (106 m3)	TN Load (t)	TP Load (t)	Flow		TN		TP
				Flow (106 m3)	Rate of Increase (%)	Load (t)	Rate of Increase (%)	Load (t)	Rate of Increase (%)
2017	22.8	40.9	6.6	23.7	3.7	44.1	7.3	6.8	3.5
2018	19.3	55.7	6.8	20.1	4.0	59.5	6.4	7.1	3.7
2019	26.9	62.2	7.5	27.9	3.6	67.6	8.0	7.8	3.4
2020	33.9	75.4	8.2	35.2	5.2	83.1	9.3	8.6	4.8
2021	18.4	57.1	5.7	19.1	3.5	59.8	4.5	5.9	3.2

) that the epikarst zone plays an incremental role in the discharge and N and P load of the drainage outlet, with the discharge increasing by 3.5~5.2%, and the total N and P load increasing by 4.5~9.3% and 3.2~4.8%, respectively. The increment in total nitrogen load was the largest, and the increment in total phosphorus load was between 0.81 and 2.33 kg/hm². The increment of total phosphorus load was the smallest. The increment of total nitrogen load was mainly due to the loss of dissolved nitric nitrogen, which is easy to migrate with karst groundwater. The increment of total phosphorus load was between 0.06 and 0.12 kg/hm², which may be mainly because the adsorption state of phosphorus is mainly migrated with surface runoff.

Table 3. Comparison of simulation results of nitrogen and phosphorus in exposed karst landform area between SWAT and modified SWAT model.

Model	Year	NSUPQ (Nitrate Migrating to Rivers from HRU Surface Runoff)	NLATQ (Nitrate Migrating to Rivers from Lateral Flow)	NO3GW (Nitrate Migrating to Rivers from Groundwater Flow)	P_GW (Dissolved Phase Phosphorus Migrating to the River from the HRU Groundwater Flow)	SOLP (Dissolved Phase Phosphorus Migrating to the River from Surface Runoff)	OrgN (Annual Amount of Organic Nitrogen)	OrgP (Annual Amount of Organic Phosphorus)
				(kg/ha)
SWAT	2017	0.011	0.451	0.023	128.219	0.514	182.235	26.181
	2018	0.032	0.509	0.159	458.593	0.481	182.410	23.604
	2019	0.011	0.345	0.106	620.035	1.743	233.084	41.187
	2020	0.023	0.468	0.185	719.372	2.387	223.227	35.330
	2021	0.045	0.416	0.141	494.558	0.560	191.885	22.140
Modified SWAT	2017	0.006	0.564	0.031	130.836	0.534	260.335	26.196
	2018	0.023	0.606	0.198	472.776	0.490	246.501	23.611
	2019	0.006	0.426	0.131	645.87	1.756	306.689	41.207
	2020	0.011	0.602	0.237	734.053	2.432	286.188	35.354
	2021	0.022	0.501	0.174	509.854	0.598	231.187	22.281

lists the amount of nitrate migrating to rivers from HRU surface runoff (NSURQ), lateral flow (NLATQ), and groundwater flow (NO₃GW), the amount of dissolved phase phosphorus migrating to the river from the HRU groundwater flow (P_GW) and surface runoff (SOLP), and the annual amount of organic nitrogen and phosphorus (OrgN and OrgP) migrated from the sub-catchment to the river in the exposed karst landform area between the SWAT and modified SWAT models.

4. Discussion

4.1. Flow Simulation

The flow simulation results show that the modified SWAT model more accurately reflects the hydrological process in which the drainage responds quickly to precipitation and holds less water after regulating and storing the epikarst zone. By adding the epikarst module, the modified SWAT model showed higher accuracy in capturing the hydrological characteristics of the cross-distributed areas of bare and covered karst landforms and improved the adaptability of the SWAT model to simulate hydrology in the southwest karst area. This method is consistent with the research conclusion of Wang et al. [55], which further verifies the key role of hydrological module adjustment based on specific geomorphic features on the adaptability of the SWAT model.

In comparison with other modifications proposed for SWAT in karst areas, our approach specifically targets the rapid response and storage characteristics of the epikarst zone. For instance, some studies focus on simulating ecohydrological processes under different vegetation coverages [56,57], while others concentrate on integrating groundwater flow dynamics using MODFLOW to better simulate the interaction between surface and underground flows in karst systems [27,58]. Although these methods enhance the model’s performance in simulating vegetation dynamics or its ability to represent vertical groundwater flow in rock layers, our modification emphasizes the horizontal flow dynamics within the epikarst layer, which is crucial for accurately capturing the rapid hydrological responses observed in Southwest China with alternating exposed and covered karst landforms.

4.2. Nutrient Load Simulation

The improvement of the SWAT model is mainly aimed at the horizontal and vertical migration processes and the path of surface sources of nitrogen and phosphorus with water in the exposed karst geomorphic area of the Mudong watershed. The bare karst geomorphic features have a certain positive effect on nitrate discharge in karst groundwater within the hydrological response unit but have no obvious effect on soluble phosphorus and organophosphorus discharge in karst groundwater. Zhang et al. [59] also found that the nitrate migration process in the southwest karst area was closely related to the water flow characteristics of different types of karst fracture pipes. The results show that the modified SWAT model based on karst geomorphologic features has advantages in capturing the dynamic processes of nitrate, a key agricultural non-point source pollutant.

In contrast to other modifications that focus on nutrient algorithms and transport pathways (such as incorporating new nitrification and denitrification algorithms [53], or soil–water process modules to affect nutrient fluxes [60]), our study specifically addresses the unique challenges posed by the karst topography in Southwest China. This approach supplements the improvement of the SWAT model by reflecting how the characteristics of the epikarst zone influence the transport of nutrients, particularly the output characteristics of nitrate under rapid hydrological responses in the southwest karst system.

5. Conclusions

In this study, we addressed the interspersed exposed and covered karst landform characteristics in agricultural watersheds of the southwest karst region and improved the SWAT model to more accurately simulate the transport of agricultural non-point source pollutants of nitrogen and phosphorus. By incorporating an epikarst module and adjusting the delineation of hydrological response units, the SWAT model can more reasonably reflect hydrological and non-point source pollution processes in karst watersheds. The main findings are summarized as follows:

(1)

Effectiveness of the Modified Model: The application of the improved SWAT model in the Mudong River watershed demonstrated good simulation performance. The prediction accuracy of monthly runoff, total nitrogen, and total phosphorus showed that the R² values of the improved model increased by approximately 6.8%, 10.3%, and 9.7%, respectively, and the NSE increased by 14.8%, 11.3%, and 9.9%, respectively. These improvements underscore the model’s enhanced applicability and precision in simulating hydrological and non-point source pollution dynamics within karst watersheds.

(2)

Hydrological Significance of Karst Features: The comparison between the original and modified SWAT models revealed that the latter captured the rapid hydrological response of karst watersheds to precipitation events more precisely. This underscores the significance of incorporating karst landform characteristics into hydrological modeling.

(3)

Specialty of Nitrate Migration under Karst Features: The modified SWAT model demonstrated a more pronounced positive response in simulating nitrate export under exposed karst landform conditions, in contrast to the simulation outcomes for soluble phosphorus and organic phosphorus. This highlights the necessity of tailoring model improvements to the specific landform attributes of karst regions to understand better and predict the dynamics of agricultural non-point source pollutants.

In conclusion, although this study has achieved notable progress in model refinement and application, there is still potential for further optimization and expansion. Future work will focus on two main areas: (1) We will further verify and refine the model parameters to enhance the model’s adaptability to various karst geomorphological and environmental settings; (2) We will examine the model’s performance in simulating different agricultural management strategies in dynamic environments. Our goal is to provide more targeted guidance for minimizing agricultural non-point source pollution and improving watershed management. These efforts are expected to contribute to more scientifically sound and effective decision-making for water resource management and environmental conservation in karst watersheds.