Improving the Applicability of the SWAT Model to Simulate Flow and Nitrate Dynamics in a Flat Data-Scarce Agricultural Region in the Mediterranean

Abstract

Understanding the soil and hydrologic processes in agricultural watersheds are vital for reliable assessments of water quantity and quality to support integrated river basin management. However, deriving hydrology-relevant information is complicated in flat data-scarce agricultural watersheds due to constraints in watershed delineation, flat topography, poor natural drainage, and irregular irrigation schedules by human intervention. The study aimed to improve the applicability of the Soil and Water Assessment Tool (SWAT) model to simulate daily flow and NO₃ concentrations in a flat data-scarce agricultural watershed in the Lower Seyhan Plain (LSP) in Turkey. Refined digitized stream networks, discharge data derived from fully equipped gauging stations, and satellite data (Landsat 7 ETM+, Aster GDEM, etc.) had to be integrated into the modeling process to improve the simulation quality. The model was calibrated using a 2-year (2011–2012) dataset of streamflow and NO₃ using the Sequential Uncertainty Fitting (SUFI-2) approach and validated from 2013 to 2018. Daily water yields were predicted with a reasonable simulation accuracy (E values ranging from 0.53 to 0.82 and percent bias (PBIAS) from 0 to +4.1). The results proved that integrating redefined stream networks to SWAT within a Geographic Information System (GIS) environment increases the simulation capability of flow and nitrate dynamics efficiently. Automated delineation of these networks and sub-basins at low topographic transitions limits the SWAT accuracy.

1. Introduction

Mesoscale hydrological models can be useful tools to answer land use and water management related questions in agricultural regions worldwide, mostly with mountainous and hilly topography [1]. However, the application of hydrological models in flat agricultural watersheds is still challenging due to anthropogenic disturbance of drainage systems (e.g., alteration of irrigation water, hydroelectricity) and plain topography. Drainage network extraction from the digital elevation model (DEM) is the critical input for delineating the number of sub-units, which affects the spatial variability of the datasets and the simulated spatial distribution of processes in the entire watershed [2]. Using the stream networks derived from the DEM in automated delineation procedures in hydraulic-hydrological models can lead to inaccuracies in low-relief regions [3,4,5]. These inaccuracies are mainly caused by utilizing the DEMs with a low resolution, which are usually free. DEMs with higher resolution (Light Detection and Ranging (LIDAR), IKONOS stereo images, etc.) could provide more accurate and precise elevation information; however, they are mostly expensive. Commonly used low resolution DEMs (i.e., Shuttle Radar Topography Mission (SRTM)) might ignore detailed topographical information and affect the water flows with minimal elevation gradients in these regions [6]. This is especially true for the Soil and Water Assessment Tool (SWAT); [7], which has not been frequently applied to coastal agricultural plains. Schmalz et al. [8] tested the capability of using SWAT in the flat lowland catchments of northern Germany with low hydraulic gradients. They found that topography was one of the most critical inputs for parameterization directly affecting the model behavior.

Habeck et al. [9] simulated the nitrogen flows well in the Nuthe low land catchment in northeastern Germany. They used the Soil and Water Integrated Model (SWIM) [10], which is based on SWAT and MATSALU [11] and uses a similar spatial disaggregation scheme like SWAT (basin—subbasins—hydrotopes). Water flows, vegetation growth, and nutrient cycling processes are calculated for each hydrotope (or Hydrological Response Unit (HRU) in SWAT) and then aggregated at the subbasin level, from where lateral flows of water and nutrients are routed to the basin outlet [12].

In conceptual semi-distributed hydrologic models such as SWAT, a sound structuring of watersheds into smaller units representing heterogeneity and topography of the region is essential for achieving accuracy [13]. Khafaji and Sweiti [14] assessed the complemental interactive effects of topography using DEM and land cover on the estimated runoff using SWAT in three watersheds in Iraq. They concluded that the watershed delineation and stream networks are profoundly affected by the DEM’s quality, particularly in flat watersheds with a limited slope and low gradients variance. Digitized stream networks derived by a stream burning method to artificially lower and fit the relief data along the existing channels in the catchment can account for this problem recently [15]. This technique can increase the accuracy of representing the input channels in the model—in case the stream networks may not be consistent with the DEM due to the flat topography of the study area [3].

Data scarcity is another problem in hydrological modeling [16]. The shortage of monitoring data generally limits the calibration of the model parameters, which are typically adjusted against observed time series values [17,18,19,20]. This is particularly challenging in Turkey, where gauging stations are sparse in agricultural watersheds, discharge measurements are limited, and the quality of measured data is generally low [20,21]. Although setting up well-equipped gauging stations is expensive, it is essential to modelling to gain hourly or daily discharge data in data-scarce watersheds [22,23].

The study investigates the suitability of SWAT to simulate daily flow and NO₃ concentrations in a flat, data-scarce agricultural watershed in the Lower Seyhan Plain (LSP) located in the Eastern Mediterranean part of Turkey. By redefining digitized stream networks, using discharge data derived from fully equipped gauging stations, and integrating Landsat 7 ETM+, Aster GDEM satellite images into the modeling procedure, we show how to overcome data scarcity and the problems caused by the flat topography of the study region. We tested how far the use of the redefined digitized stream networks and the additional data would improve the model’s capability to simulate the water flows, which would be essential to recommend best management practices for sustainable water use. Moreover, the changes in the NO₃ cycle were simulated to check if the model can represent both the spatial distribution of nutrient concentrations and the number of days with nitrate and water stress that affect plant growth. The study shows ways to improve the applicability of SWAT and similar models in flat data-scarce agricultural regions in the Mediterranean to predict flows and nutrient loads. It contributes to integrated river basin management in such areas.

2. Materials and Methods

2.1. Study Area

The LSP is located at the Eastern Mediterranean coast of Turkey. It covers 2132 km², with elevations ranging from 0 to 23 m. Around 81% of the entire basin area is irrigated. Open irrigation systems and drainage infrastructures are managed by the State Hydraulic Works (DSI) of Turkey. Mediterranean climate dominates the region with an annual total of 730 mm precipitation and a mean annual temperature of 30 ℃ for summer and 7 ℃ for winter. Due to the privileged climate and fertile soils, LSP is one of Turkey’s most productive agricultural areas, providing approximately 10% of the total agricultural production. Large-scale irrigated agriculture extends throughout the lower delta, where maize, citrus fruit, cotton, wheat, and vegetables are cultivated and exported to the European Union, Russia, and other neighboring countries. These crops depend on water provision supplied by reservoirs that store the runoff of winter, rain, and snow in the upper mountainous areas.

The Seyhan Dam supplies irrigation water to the LSP through a two-conveyance canal and the Seyhan River. The flow downstream of the diversion dam is used as one of the main sources for irrigation and is deficient. This water deficiency will be increased according to the rising demand to maintain the river ecosystem in the future [24]. Besides, the impacts of fertilizer load on salinity, water pollution, and other environmental issues are still unclear, especially on the Ramsar wetlands, including Akyatan and Agyatan Lagoons of Cukurova Delta along the coastal belt of LSP. A Ramsar site is a wetland site designated to be of international importance under the Ramsar Convention. Thus, the LSP represents an ideal testbed for examining the capability of the SWAT model to simulate hydrology and water quality in such a complex agricultural plain. Figure 1 shows both the location of the LSP and the gauging stations used in the study.

There are four main discharges of the drainage networks of the LSP (three to the Mediterranean Ocean, one to Akyatan Lagoon). We set up four real-time stream gauging stations (funded by the Turkish Scientific and Research Council (TUBITAK), Project No: 115Y063) at these locations to accurately measure the accumulated drainage flow from the entire study region. Each station included a Limnigraph and a multi-parameter device that measured Electrical Conductivity (EC), temperature, Total Dissolved Soils (TDS), Turbidity, Dissolved Oxygen, pH, and Nitrate on an hourly basis. A data logger on the station transferred the measured data to the server in our Remote Sensing and Geographic Information System laboratory at Cukurova University via data modem and General Packet Radio Service (GPRS). A separate concrete and steel construction platform was built for the Yemisli station (at Akyatan Lagoon) to cover all of the water inlets to the lagoon, while other stations were built on fixed structures such as bridges. The water level is measured by the radar system instead of a limnigraph at Bebeli station since the drainage channel in this area consists of mud after 2–3 m depth; this increased the station cost.

2.2. Model Description

SWAT is a process-based hydrological model that runs on a daily time step to simulate the impact of agricultural management practices on water, sediment, plant growth, nutrient cycling, and agricultural chemical yields at the catchment scale [7]. It has been applied to different sizes and scales of watersheds worldwide [25,26,27].

The model subdivides a watershed into sub-basins based on topography, which are connected to stream networks. Sub-basins are further partitioned to HRUs, which are heterogeneous land units with land cover, soils, slope, and aspect. The simulation outputs, including soil water content, surface runoff, nutrient cycling, and management practices, are estimated on an HRU basis. Weighted averages are applied to aggregate the results for the sub-basins.

The hydrological cycle is simulated in two major sections as land and routing phases. Runoff, erosion evapotranspiration, nutrient, and soil-water processes are simulated in the land phase. The routing phase can define the water movement through the channel network to the outlet [28]. Surface runoff from daily precipitation is estimated by the Soil Conservation Service Curve Number method (SCS-CN) [29]. Evapotranspiration was calculated using the Penman–Monteith method. Outflow from a channel is adjusted for transmission losses and return flow [7,28]. Stream networks, as model inputs, have critical importance in the SWAT modeling process. A command structure is used for routing runoff through streams and inputting measuring data on point sources. In this case, the model can simulate the discharge in the study basin by subdividing it into sub-watersheds based on these stream networks [30,31,32,33].

2.2.1. Model Inputs

Figure 2 and

Table 1. Spatial and temporal datasets used in the study.

Data	Data Type	Source	Scale/Resolution	Data Description/Properties
Spatial Data	DEM	ASTER GDEM	15 m	Stream network, sub-basin delineation, slope derivation
	Satellite images	LANDSAT 7 ETM+	30 m	LULC, cultivated parcels
	LULC	Classified from LANDSAT 7 ETM+	30 m	15 land cover/use classes. 95% (Kappa statistic)
	Soil	Derived from TUBITAK Project (No: 115Y063)	1/250.000	17 soil series Soil physical properties including depth, saturated hydraulic conductivity, texture
	Stream networks	Digitized using stream burning method	30 m	Including main drainage channels
	Hydrological Response Unit (HRU)	Derived using QGIS	30 m	Comprising spatial information on heterogeneous units
Hydro-meteorological Time-series	Temperature	Turkish State Meteorological Service	13 stations	Daily climate data (2011–2018)
	Relative humidity
	Wind speed
	Precipitation
	Solar radiation
	Observed flow	Real-time stream gauging stations by TUBITAK Project (No: 115Y063)	Four gauging stations	Flow rates were calculated using this 15 min flow rates (ft3 s−1) to Derive mean daily outflow rates (m3 s−1) and daily total (mm) for each sub-basin (Sari, 2018).
	NO3 records	Real-time stream gauging stations by TUBITAK Project (No: 115Y063)	Four gauging stations
Land Management Information	Management/Point sources	Water use associations, expert opinion, and field campaigns		Planting, harvest, irrigation applications, fertilizer/pesticide/tillage applications for each crop

provide an overview of the full range of spatial and temporal datasets that have been used as model input. Hydrological Response Units (HRUs) were derived independently by overlapping slope, soil, and land use/land cover (LULC) in quantum GIS and integrated into the modeling (see Figure 2). Spatial inputs of the model were derived using the free and open source Geographic Information System quantum GIS (QGIS).

For each sub-basin, a monthly observed flow point source dataset was also used. The Water Use Associations provided management operation data for the various crops in the study area for each crop (e.g., corn, citrus, wheat, soybean, peanut, vegetables) and added them into the database. The management operations included planting, harvest, irrigation, fertilizer and pesticide applications, and tillage operations.

2.2.2. Data Adjustments, Model Parametrization, and Setup

Ancillary Data

Ancillary data used for model parameterization in the study comprised LULC data and soil information. SWAT already includes an extensive default database that we adjusted using four data sources: Landsat 7 ETM+ images, topographic maps, State Hydraulic Works land cover records, and ground data from field surveys. LULC maps were derived seasonally using six Landsat 7 ETM+ images acquired in June, July, and August 2016 (summer) and December, January, and February 2016 (winter). These images were selected according to the harvest, sowing, and tillering periods of the crop plants and geometrically corrected to the Universal Transverse Mercator (UTM) coordinate system and WGS 84 datum. An object-oriented classification was used to derive seasonal LULC maps for winter and summer. The images were segmented at the first step of LULC classification and processed using the nearest neighbor technique as corn, cotton, grassland, forest, and settlements and used in the model. Corn covered 35% of the LSP, where cotton was planted on 30% of the area.

Besides the LULC information, an extensive and detailed soil database was used in the study; a 1:250.000 scaled soil map was provided by our project (no: 115Y063) supported by the Turkish Scientific and Research Council (TUBITAK). The shapefile was converted to a raster format (30 m resolution) in a GIS environment. Thirteen soil series were defined, and model parameters were included as depth, pH, water holding capacity, total porosity, sand, clay and silt percentages, hydraulic permeability coefficient, organic matter amount, and electrical conductivity. The definition of the soil parameters was supported and proven by field campaigns within our TUBITAK project (no: 115Y063).

Drainage Networks

During the model setup, numerous problems occurred due to the flat topography of the LSP. The elevation of the terrain ranges from 0 to 40 m, maximum slopes are only 5%. The pre-processing procedure in SWAT utilizes the DEM data to delineate the sub-basins, which is problematic in such flat regions as the LSP. For instance, the delineated sub-basins and stream networks did not represent the real conditions since the flat terrain also caused weak and inconsistent (delineated) stream networks. The SWAT pre-processing procedure offers to burn in streams in its standard settings. However, our attempts to derive these networks using SWAT were unsuccessful. The drainage networks mostly cause this failure to settle in contrary to slope direction and their disconnected branches. Thus, SWAT was unable to derive these networks satisfactorily due to their inconsistent flow direction (see also next subsection on sub-basin delineation). Hence, these networks were externally derived (see Figure 3) and included in the modelling process.

Moreover, the amount of water supplied daily to the channels was included in the model system as point sources. The State Hydraulic Works provided these point source data. The plant management database was also updated with current data from the field to define the crop-specific agricultural management practices in SWAT.

Sub-Basin Delineation

In the first step of the SWAT pre-processing procedure, the LSP was divided into four sub-basins (Baharli, Karagocer, Yemisli, Bebeli). We used the gauging stations on the stream network as main outlets of divided sub-basins. Four gauging stations set up for the study were used as outlets; sub-basins were delineated in a GIS environment and introduced into the model along with sub-basins. Furthermore, stream networks were delineated using default stream burning in SWAT. However, the sub-basins were not successfully derived since the stream networks were mostly disconnected. As previously mentioned, externally derived drainage networks were used, and the subbasins were then properly delineated using stream burning (Figure 3).

The coverage of the sub-basins delineated in the study ranged between 12,945 ha (Yemisli) and 61,439 ha (Baharli). The main channel of these sub-basins discharges directly into the Mediterranean Sea. They are fed by the return irrigation flow from agricultural lands. Gauging stations as the main outlets of each basin were installed on the bridges over the drainage networks.

Hydrological Response Units

HRUs were delineated for each sub-basin by overlapping spatial inputs as LULC, soil, geology, and topography. Each class layer was reclassified to derive different classes of hydrological significance. The total number of HRUs was reduced by eliminating units with redundant information to increase the model performance’s efficiency. HRU definition included minimum area thresholds (50 ha) for the LULC, soil, and slope to eliminate negligible areas and helped to reduce the total number of HRUs to increase the computational efficiency of the model simulations. The minimum threshold for HRU delineation was defined as 50 ha. Overall, 125 HRU for Baharli, 102 for Karagocer, 83 for Yemisli, and 107 for Bebeli sub-basins were delineated.

Time Series and Management Data

Meteorological data used in the study comprised precipitation (mm), minimum, mean, maximum temperature (℃), wind speed (m/s), net radiation (MJ/m²/day), and humidity (%). These data were derived from 5 climate stations in the LSP, including Adana (no: 17351), Karatas (no: 17981), Karaisalı (no: 17936), Ceyhan (no: 17960), Yumurtalik (no: 17979), and Seyhan (no: 17352), and was imported into SWAT.

Management data included plant sowing time, irrigation frequency and type, fertilization amount and type, soil tillage frequency considering the crop rotations [22,34]. The required information for each crop was obtained from the Turkish Ministry of Food and Agriculture database. Data on crop management included irrigation, pesticide, harvest, tillage, annual grazing operations in plant growing, and harvesting seasons. Special attention was given to citrus production that covers over 50% of the LSP. The fertilizer amount for each Citrus tree was defined as 10 kg/per tree in October and November (Sari 2018).

After setting up the spatial and time-series inputs, the model was set up and run in three main steps; (i) model simulations with default parameters, (ii) model calibration, and (iii) model simulations with improved parameters for the entire basin. The SWAT model was run monthly for six years from 2011 to 2018 with a two-year-warm-up period to minimize possible simulation and initialization problems. Streamflow and nutrient calibration was carried out from 2011 to 2012. The model was then validated for 2013 to 2018 using the data obtained from the gauging stations for each sub-basin.

2.2.3. Model Calibration and Sensitivity Analysis

The program SWAT-CUP (SWAT Calibration and Uncertainty Procedures; [35]) was used for calibration, validation, and uncertainty analysis of the model. From the four approaches SWAT-CUP offers, we used the acknowledged and robust SUFI-2 routine [28,35,36] for model calibration and uncertainty analysis of each ensemble member. In SUFI-2, a multivariate uniform distribution in a parameter hypercube was used to describe parameter uncertainty. The cumulative distribution of the outputs variables was used to derive model uncertainty of model outputs (Abbaspour et al. 2007; Strauch et al. 2012). Global and temporal sensitivity analyses were applied to define the most and least sensitive parameters within the broad parameter set. Parameters higher than a p-value of t-stat were specified. In the first step of the calibration and sensitivity procedure, an objective function is defined, and a summation form of the squared error was selected in Equation (1):

$$A=\frac{1}{{\sigma }_{low}}{\displaystyle \sum }_{t=1}^T{\left(y_{tlow}-f_{tlow}\right)}^2+\frac{1}{{\sigma }_{low}}{\displaystyle \sum }_{t=1}^T{\left(y_{thigh}-f_{thigh}\right)}^2$$

(1)

where y_t and f_t are the daily observed and simulated flow. These variables are divided into two subsets by the threshold of 2 m³ s⁻¹. This threshold represents the mean streamflow during the calibration. If the y_t is lower than it, y_t and f_t are belonging to subset y_tlow, f_tlow, otherwise high. The standard deviation of the lower (q_low) and higher (q_high) was used as weights to prevent missing base flow in optimization [28]. In SUFI-2, the spread of uncertainties in the parameter ranges expresses the model output variables with 95% (95 PPU) probability distribution. It utilized Latin hypercube sampling with several iterations. For each iteration, SUFI-2 aims to achieve the reasonable values of P (the percentage of observed data surrounded by 95 PPU) and R (the thickness of the 95 PPU envelope) factors for the most accurate results.

Flow and NO₃ parameter sets were selected in the SWAT-CUP calibration procedure, and the model was initially run 500 times. In the global sensitivity analysis, the number of parameters was decreased to increase the calibration efficiency. While determining the parameters sensitive to these changes, the difference between the t-stat value (the coefficient divided by the standard error of a parameter) and the P-value (testing whether each parameter coefficient is equal to zero) was considered. The outputs of the global sensitivity analysis applied in the study for flow and NO₃ are given in Figure 4 and Figure 5.

The global analysis identified the channel parameters (channel width, length, and depth) as most sensitive and the HRU parameters as least sensitive for the flow-related parameter. As a result of these analyses, 35 sensitive parameters were set for flow data, while 11 parameters were set for NO₃. Moreover, SWAT-CUP was run separately for each parameter, and the sensitive parameters were determined by temporal precision (Figure 6). The calibration procedure was carried out five times for each parameter, where the most sensitive ones were determined. The main reason for the sensitivity of SWAT to the channel parameters was the region’s flat topography.

Model performance was evaluated by using standard criteria suggested by Moriasi et al. [37] and Cakir et al. [38], i.e., r² (Equation (2)), E (Equation (3)), percent bias (PBIAS) (Equation (4)) were used as the objective functions in model performance evaluation.

r²: The coefficient of determination [39].

$${\mathrm{R}}^2=\frac{\left[{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}({\mathrm{Q}}_{\mathrm{m}}^{\mathrm{t}}-{\mathrm{Q}}_{\mathrm{m}}^{-}\right){\left({\mathrm{Q}}_{\mathrm{o}}^{\mathrm{t}}-{\mathrm{Q}}_{\mathrm{o}}^{-}\right)}^2]}{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\mathrm{Q}}_{\mathrm{m}}^{\mathrm{t}}-{\mathrm{Q}}_{\mathrm{m}}^{-}\right)}^2{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\mathrm{Q}}_{\mathrm{o}}^{\mathrm{t}}-{\mathrm{Q}}_{\mathrm{o}}^{-}\right)}^2}$$

(2)

E: Nash-Sutcliffe coefficient (Nash Sutcliffe, 1970).

$$E=1-\frac{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\mathrm{Q}}_{\mathrm{o}}^{\mathrm{t}}-{\mathrm{Q}}_{\mathrm{m}}^{\mathrm{t}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{T}}{\left({\mathrm{Q}}_{\mathrm{o}}^{\mathrm{t}}-{\mathrm{Q}}_{\mathrm{m}}^{-}\right)}^2}$$

(3)

PBIAS [40].

$$\mathrm{PBIAS}=100\times \frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Q}}_{\mathrm{o}}-{\mathrm{Q}}_{\mathrm{m}}\right)}_{\mathrm{i}}}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{Q}}_{\mathrm{o}}}$$

(4)

where Q_o Q_m = mean observed discharge, Q_m = modeled discharge, Q^t = discharge at time t. Since only one technique might cause misleading results, all those functions were used for the performance evaluation of the model [41].

3. Results

This study intended to improve the SWAT model’s performance in a flat Mediterranean agricultural basin by introducing the digitized stream network input, Landsat 7 ETM+ satellite images, and involving gauging data from the stations set up for this research. The model was first run in daily time steps using default parameter conditions for each sub-basin and nutrients, runoff, actual evapotranspiration (ET), and water stress. The model was then calibrated in a proxy sub-basin in the LSP and validated within the entire basin.

3.1. Model Runs Using Default Parameter Conditions

Using the default parameters has led to the prediction of excessive or insufficient fertilization according to the changes in plant nitrogen and phosphorus stress days. In this regard, a particular effort was put into the model setup to indicate mean annual stress days, including water, temperature, nitrogen, and phosphorus, for four sub-basins in the study area (

Table 2. Mean annual stress days in the four sub-basins of the Lower Seyhan Plain (LSP).

Stress Types (Number of Days/Year)	Sub-Basins
	Baharli	Karagocer	Yemisli	Bebeli
Water Stress	73.85	84.62	59.87	56.38
Temperature Stress	88.86	89.09	99.81	84.44
Nitrogen Stress	73.72	40.66	67.15	50.89
Phosphorus Stress	14.33	8.45	6.27	8.87

). The average number of nutrient and water stress days per year is crucial since it directly affects plants’ growth.

Model simulations indicated that nitrogen and water stress days were relatively high in four sub-basins. The reason for high nitrogen and water stress days can be referred to as the dates of sowing and the variation of the harvesting dates. Additionally, the predominant clay soils in the study region limited the infiltration of irrigation water provided by rain and drainage channels.

Table 3. Total annual amount of the hydrological variables of the sub-basins in the LSP.

Hydrological Variables (mm)	Sub-Basins
	Baharli	Karagocer	Yemisli	Bebeli
Precipitation	802.8	766.3	766.3	628.2
Surface Runoff	374.84	377.81	339.6	244.31
Shallow Aquifer	57.68	64.01	163.68	NAN
Deep Aquifer Recharge	4.23	4.57	9.92	2.12
Total Aquifer Recharge	84.55	90.52	198.44	42.47
Total Water Yield	436.93	446.46	513.74	268.26
Percolation Out of the Soil	83.81	89.52	196.96	42.22
Actual Evapotranspiration (ET)	356.6	323.2	386.9	352.1
Potential Evapotranspiration (PET)	1334.4	1376.5	1380.1	1399.9

shows the total annual amount of the hydrological variables of the sub-basins in the LSP.

The simulated amount of water passing from the soil to the shallow and deep aquifers was limited. Discharge from groundwater to the stream was also low due to the impermeability and porosity structure of the clay soil. Annual total surface runoff was estimated between 244 to 374 mm, and surface runoff and precipitation ratio ranged from 30% to 45%, which can be considered reasonable for a Mediterranean region [22]. The total water yield was estimated as 513 mm for Yemisli sub-basin, whereas the highest runoff was estimated in the Karagocer sub-basin. Monthly water distributions of the sub-basins are given in Figure 7.

Evaporation and humidity reached their highest levels in the summer and the lowest in the winter. The surface flow depended on rain and irrigation water, which were at the lowest in summer and highest in winter. The surface runoff in summer was mostly due to irrigation contribution to the hydrological system. The impermeability of the clay soil led to an increase in the amount of evaporation and surface runoff. In addition to these hydrological dynamics, NO_3, nitrogen, and phosphorus of the sub-basins were also simulated (

Table 4. Simulated mean annual amounts of nutrients for the sub-basins of the LSP.

Average Annual Basin Values, Nutrients (kg/ha)	Sub-Basins
	Baharli	Karagocer	Yemisli	Bebeli
Organic N	9.881	1.604	1.145	1.404
Organic P	1.205	0.221	0.156	0.188
NO3 yield (SQ) (to Stream in Surface Runoff)	3.399	2.679	2.872	2.173
NO3 yield (LAT) (to Stream in Lateral Flow)	0.009	0.012	0.025	0.009
Sol-P yield	0.033	0	0.131	0.049
NO3 leached	76.846	70.083	99.665	66.547
P uptake	NaN	12.379	14.411	17.754
N uptake	47.895	38.251	48.681	57.843
P-Fertilizer applied	7.402	4.286	4.773	5.764
N-Fertilizer applied	259.672	226.447	262.881	226.179
Denitrification	27.688	14.863	15.926	18.676
NO3 in rainfall	0	0	0	0
Initial NO3 in soil	52.219	53.21	55.548	52.412
Final NO3 in soil	225.04	202.345	243.52	214.992
Initial org N in soil	8721.166	82.616	85.152	82.174
Final org N in soil	8130.764	135.558	193.541	149.786
NO3 in fertilizers	16.823	9.742	10.848	13.101
P uptake in yield	2.88	1.58	2.82	3.972
N removed in yield	19.638	12.615	22.263	31.812
Ammonia volatilization	102.663	111.7	117.634	99.893
Ammonia nitrification	140.351	105.306	134.719	113.005

).

Seasonal crop patterns were derived using Sentinel-2 data and the records of the Water Use Associations (WUAs) of the region. During the calibration step, fertilization and irrigation parameters, which significantly affect streamflow and nitrate, were derived from farmer’s questionaries. Accordingly, it was defined that the crops used approximately half of the nitrogen for the sub-basins. Therefore, the simulated nitrogen and NO₃ days were high, while phosphorus stress days were optimum. About half of the used nitrogenous fertilizers were lost due to denitrification. Most of them were accumulated in the soil or drainage system. Therefore, the solubility of nitrogen in the surface stream was relatively high. Simulations showed that the amount of NO₃, nitrogen, and phosphorus in the soil was relatively high compared to the initial amounts. These issues indicated an over-fertilization in the sub-basins.

3.1.1. Model Calibration

For the model calibration, a proxy basin was selected by referring to the LSP’s physical conditions. Then, calibration was performed in Bebeli sub-basin as a proxy area bordered by the Berdan and Seyhan Rivers and located at the western part of the LSP. Bebeli sub-basin was selected as a test area since it represented the common crops and water regime of the entire basin. The validation process was applied for the four sub-basins using the parameter set obtained from the calibration study in Bebeli. Calibration and validation processes were applied using gauging station records derived from the stations constructed for each sub-basin. The calibration procedure (period 2011 to 2012) has been carried out for the left bank irrigation, water, and NO₃ budget of the LSP.

The SWAT-CUP calibration module was run 12 times to achieve the best value range for 35 parameters, and 800 simulations were performed for each run after determining the temporal parameters for the calibration of the flow. The most accurate result was derived in run 8. Then, these 35 parameters were used in the modeling and validation for the whole basin. The results for the E value was 0.44, and the r² value was 0.47. The relatively low values are caused by the flat topography’s influence and the complex drainage network of the study region.

The drainage networks supply water from Seyhan and Catalan dams from April to November. These frequent water inlets were defined as point sources in the model and had a significant effect (30%) on the simulated flow data. Moreover, the previously described additional information on soils, routing, groundwater, HRUs, and management parameters were also used for the calibration studies. The additional information gained for the test area was used in the validation process for the entire basin.

After determining 46 parameters (11 for NO₃ data, 35 for flow data), the SWAT-CUP was run 15 times to derive the optimal ranges for these parameters. Nine hundred simulations were performed for each run. The most accurate result was achieved in run 15 with an E value of 0.45 and an r² value of 0.46.

A comparison of the observed, uncalibrated, and validated flow and NO₃ outputs are given in Figure 8 and Figure 9.

The calibration process improved the accuracy of the model simulations compared to the uncalibrated results, as shown with the nitrate simulations in the Bebeli sub-basin (improvement of the coefficient of determination r² from 0.48 to 0.85, where E was increased from −0.18 to 0.82). A similar performance was achieved for other sub-basins as well: In Yemisli, E was increased from −0.44 to 0.72. In Karagocer, r² was improved from 0.10 to 0.61, which is reasonable for a flat agricultural basin, according to Kapur et al. [42]. Therefore, the model responded to the complex conditions of the LSP well and provided a meaningful and reasonable trend of the flow and NO₃ distribution [42,43,44].

Water

Calibrations were carried out in 2017, and validation was performed for Baharli by using the parameter ranges obtained from the Karagocer station. Numerous simulations showed a high sensitivity of groundwater and streamflow routing parameters. Selected flow parameters used in the model are shown in

Table 5. Selected flow (Flow_OUT) parameters used in the model.

Parameters (Part-A)	Best Value	Minimum Range	Maximum Range	Description
Management (Mgt)
IRR_ASO	0.4234	0	1	Irrigation/flow ratio
IRR_EFF	0.2939	0	1	Irrigation efficiency
IRR_AMT	31.51	0	100	Irrigation application depth
IRR_EFM	0.201	0	1	Irrigation productivity
AUTOWSTRS	0.179	0	1	Water stress factor of the plants
Direction (rte)
CH_D	16.891	0	30	Average main channel depth
CH_N2	0.1134	0.01	0.3	Manning roughness coefficient
CH_L2	50.281	0	500	Average main channel length
CH_W2	60.798	0	1000	Average main channel width
Groundwater (gwt)
GWQMN	1863	0	5000	Water depth coefficient in a shallow aquifer
GW_SPYLD	0.177	0	0.4	Specific yield of a shallow aquifer
GW_DELAY	15.654	0	500	Groundwater delay
GWREVAP	0.12	0.02	0.2	Groundwater “Revap” coefficient
REVAPMN	224.02	0	500	Threshold in a shallow aquifer for “revap” level depth
SHALLST	12344	0	50000	Initial depth of the shallow aquifer
DEEPST	21234.5	0	50000	Initial depth of the deep aquifer
GWHT	12.62	0	25	Water table depth
ALPHA_BF	0.92	0	1	Groundwater alpha factor
RCHRG_DP	0.877	0	1	Deep aquifer percolation factor
Parameters (Part-B)	Best value	Minimum Range	Maximum Range	Description
Management (Mgt)
TRSRCH	0.544	0	1	Quantity function of water passing through the main channel to the deep aquifer
SLSUBBSN	89.116	10	150	Average slope length
SURLAG	13.725	0.05	24	Surface flow delay
MSK_X	0.2382	0	0.3	Weight factor Affecting of storage of water entering the river
MSK_CO1	2.105	0	10	Controlling coefficient of storage timing of the water in the river
MSK_CO2	4.39	0	10	Controlling coefficient of storage timing of the water in the river
EPCO	0.395	0	1	Plant uptake factor
EVLAI.	0.121	0	10	Leaf area Index with no evaporation
Soil
SOL_K	272.43	0	2000	Soil hydraulic conductivity coefficient
SOL_AWC	0.143	0	1	Existing soil water holding capacity
SOL_CBN	8.378	0.05	10	Amount of soil organic carbon
SOL_BD	1.432	0.9	2.5	Soil bulk density
HRU
ESCO	0.0337	0	1	Soil evaporation factor
OV_N	19.654	0.01	30	Manning “n” value for surface flow
LAT_TTIME	118.52	0	180	Lateral flow movement time

.

E and r² values were 0.64 and 0.67, respectively, in the 12th run with 800 simulations, and 35 parameters were used for calibration of flow data. These results can be considered reasonable concerning the simulation problems caused by the flat topography and the complex drainage networks of the LSP (

Table 6. Simulated mean annual hydrological variables of the sub-basins.

Hydrological Variables (mm)	Sub-Basins
	Baharli	Karagocer	Yemisli	Bebeli
Precipitation	802.8	766.3	766.3	628.2
Flow generation	374.84	377.81	339.6	244.31
Contribution of the Subsurface Current to Streams	0.05	0.08	0.21	0.06
Contribution of the groundwater Current to Streams	4.27	4.53	10.01	2.15
Amount of water from the shallow aquifer to plant/soil profile	25.26	25.52	26.96	20.81
Shallow aquifer discharge	57.68	64.01	163.68	NAN
Deep aquifer discharge	4.23	4.57	9.92	2.12
Amount of water entering both aquifers	84.55	90.52	198.44	42.47
Water yield from HRUs to Streams	436.93	446.46	513.74	268.26
Infiltrated water to the soil profile	83.81	89.52	196.96	42.22
Actual Evapotranspiration (ET)	356.6	323.2	386.9	352.1
Potential Evapotranspiration (PET)	1334.4	1376.5	1380.1	1399.9

). The LSP shows a high proportion of irrigated agricultural land with a wide range of crop rotations. Irrigation networks were usually supplied from the dams in specific periods, particularly in summer, to support agricultural irrigation. These periodic water inflows defined as a point source in the model significantly affected simulated flow data and considerably improved the simulation quality. The impact of point source data on the system was estimated at approximately 30% [22].

A total annual runoff of 348.3 mm was simulated for the Bebeli sub-basin, with a maximum runoff of 515.49 mm for bare grounds and a minimum of 200.35 mm for citrus cultivated lands. The average surface flow of the Karagocer sub-basin was 320.29 mm, with a maximum flow of 465.69 mm and 189.82 mm as the lowest value (for the citrus plantations). According to the model results, Yemisli and Bebeli sub-basin had 486.74 mm and 247.18 mm flow, respectively.

Spatially-distributed flow outputs of the four sub-basins showed that the lowest flow occurred in citrus plantations since they are located in protected areas and watered using drip irrigation, limiting the runoff considerably. In contrast, the highest flow values were simulated on bare grounds due to the absence of vegetation and the soil’s low hydraulic conductivity. In other words, soil porosity and permeability were low, and the water table was high. Thus, more than 50% of the rainwater occurred as surface flow and less than 22% groundwater.

Nutrients

Flow and nitrate (NO₃) calibration were performed for each of the four sub-basins of the LSP separately. A total of 46 parameters from the test area runs were used for calibration and validation of the model, which lead to satisfactory results. NO₃ concentration in drainage networks or channels increased during the spring months when rainfall was abundant and decreased gradually in summer months. While irrigation of agricultural areas was relatively high in summer, its effect on NO₃ concentration was limited. Since the LSP is almost flat, the discharge into the drainage channels was relatively low. However, due to excess fertilization and less permeability of the clayey soils, NO₃ remained in the soil and accumulated gradually. More than 25% of the applied fertilizers were denitrified in this way. The simulated NO₃ infiltration into the soil was less than 21% of the applied fertilizer. Less NO₃ infiltration was referred to a high amount of nitrogen solubility in the surface flow and the evaporation of this flow before penetrating the soil.

Due to the tides occurring at specific periods in Baharli and Karagocer drainage channels, the salinity increased considerably, and NO₃ records also increased during these periods. The reason for using flow parameters was to simulate the transport rate of NO₃ to water accurately.

The amount of NO₃ in the soil increased with fertilizers and is directly related to plant growth and, thus, for agricultural production and management. After simulating the flow of sub-basins, sensitive parameters for calibration of NO₃ simulations were also defined. The list of sensitive parameters for nitrate simulations is given in

Table 7. Selected nitrate parameters used in the model.

Parameters	Best Value	Minimum Range	Maximum Range	Description
Basin (Bsn)
CH_ONCO_BSN	81.1828	0	100	Organic nitrogen in the channel
NFIXMX.bsn	10.255	1	20	Daily maximum nitrogen fixing
FIXCO.bsn	0.263	0	1	Nitrogen fixation coefficient
N_UPDIS.bsn	28.238	0	100	Nitrogen uptake and distribution parameters
NPERCO.bsn	0.3938	0	1	Nitrogen infiltration coefficient
RCN.bsn	1.1891	0	15	Nitrogen in the rain
Groundwater (gwt)
SHALLST_N	740.31	0	1000	NO3 water in groundwater contribution to flow
HLIFE_NGW	105.744	0	200	Half-life of NO3 in a shallow aquifer
HRU
SOLN_CON	5.157	0	10	Soluble nitrogen concentration after applying urban Best Management Practices (BMP)
Chemical (chm)
SOL_NO3	28.93	0	100	The initial concentration of NO3 in soil
SOL_ORGN	11.24	0	100	The initial amount of NO3 in soil

.

The most accurate results were derived in Run 15 with an E value of 0.45 and an r² value of 0.46. The results were satisfactory for each sub-basin (Table 5). The model performance of nitrate simulations in the sub-basins is provided in

Table 8. Model performance for calibration and validation of the flow and NO₃ simulations in the study area.

Sub-Basins	Flow						NO3
	Calibration (2011–2012)			Validation (2013–2018)			Calibration (2011–2012)			Validation (2013–2018)
	E	r2	PBIAS	E	r2	PBIAS	E	r2	PBIAS	E	r2	PBIAS
Bebeli	0.44	0.47	−0.1	0.64	0.70	+4.0	0.45	0.48	−9.4	0.82	0.85	+3.4
Yemişli	0.15	0.22	−0.63	0.60	0.68	+3.3	−0.44	0.23	−15.9	0.72	0.75	+2.9
Karagöçer	0.01	0.08	−2.9	0.57	0.66	0	0.10	0.18	−4.8	0.53	0.61	+3.8
Baharlı	0.13	0.44	−0.1	0.63	0.82	+4.1	0.21	0.43	0	0.66	0.85	+7

and Figure 10.

Since we have integrated digitized stream networks into the model, the flow was reasonably estimated even in uncalibrated mode. However, the NO₃ showed an overestimation, especially without calibration. This can be explained as the effect of anthropogenic disturbance (over-irrigation) and, therefore, nitrate percolation in the region. The calibration of NO₃ seemed reasonable and helpful to improve its simulation results in this case.

The amount of applied fertilizers was high in the areas with extensive agriculture. Most of the applied fertilizers were accumulated in the soil, while some were uptaken by the plants. According to the model outputs, the mean annual NO₃ rate in the soil of the Bebeli sub-basin was 1.931 kg N/ha. In this sub-basin, the highest NO₃ amount was defined as 17.68 kg/ha as nitrate content in citrus plantations of the Incirlik subregion and whereas the lowest was 0.1 kg N/ha in Arikli. The mean annual total NO₃ in the soil of the Karagocer sub-basin was 1,895 kg N/ha, and the NO₃ rates ranged from 0.01 to 15.22 kg N/ha. The mean annual total NO₃ in Yemisli and Baharli sub-basins are 1.486 kg N/ha and 1.500 kg N/ha, respectively. The NO₃ maps of the sub-basins derived from the SWAT simulations are given in Figure 11.

The amount of NO₃ in the soil of the four sub-basins was similar, and a large amount of NO₃ applied with agricultural fertilization was accumulated in the soil. This accumulation was mostly due to over-fertilization and the clay soil structure of NO₃ leakage to the soil.

The hydraulic conductivity coefficient and organic matter amount were relatively low in soils dominated by clay in the modeling process. In contrast, the water holding capacity of clay soil is high. During calibration, the hydraulic conductivity coefficient of the soil series and the amount of organic matter in the soil were increased by about 50% to increase the amount of lateral flow. In comparison, the water holding capacity of the soil was reduced by 20–30%. As a result, the amount of lateral flow was increased from about two mm to 15–20 mm approximately. The amount of water stored in the shallow aquifer was extended to allow water flow to reach the deeper aquifers.

4. Discussion and Conclusions

In this applied research, the problem of implementing the SWAT model in flat areas with composite irrigation networks, such as in the Mediterranean coastal areas, when high-resolution DEMs are lacking, is investigated. We showed that our suggested procedure of using accurate digitized stream networks and ancillary data could significantly improve the simulation quality of SWAT in a flat data-scarce basin. Our study has shown that a remarkable effort should be given to identify the most sensitive parameters to improve the calibration and validation of the model. Besides, longer-term time series data are of high value to achieve the model performance. The additional use of including field measurements on water quality and quantity, meteorological data, satellite data to define crop rotations, and local soil parameters are needed to enable accurate SWAT simulations on monthly and daily scales in our study area.

One of the most important steps was the calibration stage, carried out in a selected proxy sub-basin (Bebeli sub-basin) of the region. The simulated streamflow and nitrate loads improved in comparison to previous studies in this region with an overall E value of 0.69. Akgul [34] achieved SWAT simulation qualities with 0.56 E, whereas the simulations performed by Sari (2018) show an E of 0.62.

In the study of Sari (2018), the streamflow varied at Bebeli station between 200 to 300 mm annually, which is in the same range as our study (244 mm). Total annual streamflow generation changed between 339 to 375 mm. According to Tanaka et al. [45] this value is reasonable in this study area. Total annual precipitation inputs of our model application ranged from 628 to 804 mm, which are also in an acceptable range as defined in Hoshikawa et al. [46], with 744 mm per year.

During the calibration and validation procedure, many parameter changes had to be carried out to reduce the surface water flow amount and to increase the passage of the water to the drainage channels. The curve coefficient (CN) parameter affecting the surface water flow amount was decreased by 10%. Moreover, the evaporation compensation factor (ESCO) increased by 5% because the CN and ESCO parameters depend on the soil and LULC classes. These two parameters are fundamental since they determine the amount of surface flow. For instance, a high CN value indicated high surface runoff and well-developed soil structure, while the low CN value corresponds to high sub-soil percolation and low surface runoff.

Our results suggest using stream burning to derive digitized stream networks externally and include them in the SWAT modelling process to avoid inconsistencies in sub-basin delineation, which can also be used in similar regions. Basins such as the LSP with flat topographic structure, complex drainage networks, and diverse agricultural activities can cause severe hydrological modelling problems. Our literature review showed that the application of SWAT in those regions is still limited. Therefore, our study suggests methods to improve the applicability in such areas, in our case, also affected by data scarcity.

The LSP drainage network system is human made, unlike the common river network systems. Therefore, the connections of the drainage networks mostly consisted of artificial irrigation systems such as pipes or pump stations. Moreover, the subbasins delineated by SWAT stream burning did not represent the study sites since the model was not correctly able to produce the subbasins due to disconnected drainage networks. Hence, we included externally digitized stream networks to the model to enable precisely delineated sub-basins.

For deriving more precise spatial inputs, the Unmanned Air Vehicles and LIDAR systems might support the field surveys to map channel networks and sub-basin topography. Using these systems can be essential to produce more precise DEM data that can contribute to more accurate hydrological modeling in flat regions.

The sensitivity analyses showed that in our flat semi-arid basin, the most sensitive parameters were related to groundwater, surface runoff, and channel geometry. Since the model was not able to simulate drainage channels or network structure correctly, the most suitable intervals for channel parameters were set after determining the channel width, length, and depth parameters. The average depth of drainage channels is simulated between 10–50 cm.

The model results showed a high fertilization rate in the LSP, and plants used only 60% of the given fertilizer; the remaining part was accumulated in the soil. A high rate of fertilization might threaten the agricultural yield in the region in the next decades due to the degradation of the soils in the model.

It was also observed in the modeling that a small amount of nitrogen accumulated in the soil transferred to main drainage channels. The main reasons are assumed to be flat topography, clay soil structure, structural features of drainage channels, and water mobility. As a result, the calibration and validation graphs showed that the non-calibrated nitrate values were considerably lower than their observed data. The accumulation of NO₃ in the soil was ignored in model settings in calibration and validation. In this way, satisfactory results were derived in the validation of the sub-basins.

Soil and fertilizer type, quantity, and agricultural irrigation scheduling varied widely within the plant management database. This complexity may mislead to inaccuracies in model simulations. Irregular crop rotation might also cause an inconsistency in the parameters used during calibration and influence the results. Although the nitrate simulations were considered satisfactory, no point source data were used to model the nitrate loads, causing uncertainties. Understanding the nitrate and water balance of the flat agricultural plains could give us valuable insights into water management decisions and increase awareness for sustainable water resources management. It can be concluded that the SWAT can be a useful and adaptable tool to apply in flat plains of sub-tropical regions to simulate streamflow and NO₃ dynamics.