Hydrological simulation and evaluation of drought conditions in the ungauged watershed Parishan lake Iran, using the SWAT model

Global climate change is primarily characterized by a notable rise in temperature and unequal distribution of precipitation, both of which pose challenges to sustainable development (Zhao et al. 2015; Dong et al. 2018). Alterations in precipitation patterns will have a profound impact on environmental processes and the utilization of natural resources, particularly water resources (Jiang et al. 2016; Javan & Movaghari 2022).

The physical, chemical, and biological characteristics of lakes are directly influenced by changes in air temperature and precipitation. Additionally, these changes indirectly impact lakes through modifications in the surrounding watershed. This includes shifts in hydrological flow pathways, landscape weathering, catchment erosion, soil properties, and vegetation. Limnologists, who study lakes and rivers, are particularly interested in the interaction between these variables. They also focus on the feedback effects that either amplify or mitigate environmental changes, as well as the threshold effects that can cause lakes to suddenly transition from one environmental state to another (Vincent 2009).

Parishan is a permanent freshwater lake located in the southwest of Iran and unfortunately has disappeared in recent years. Based on multiple studies, the western regions of Iran have experienced the most notable increase in temperature and decrease in precipitation. This pattern is expected to persist in the future. (e.g., Najafi & Moazami 2016; Ahmadi et al. 2018; Alizadeh-Choobari & Najafi 2018; Moghim 2018; Nazeri Tahroudi et al. 2019; Abbasian et al. 2021; Doulabian et al. 2021; Radmanesh et al. 2022).

The existence of wetlands is essential due to their ability to fulfill human requirements and provide ecosystem services. Parishan lake is the only freshwater lake in the Iranian plateau and is under continued changes. Although it does not directly serve as a water source, it plays a significant role in supporting the survival of plant and animal communities. Additionally, it has controller effects on the groundwater level in aquifers and climatic balance (Zamani & Mousavi Nasab 2012). Climate change, stemming from human activities, is not the only influence on water resources. Population growth, economic development, urban expansion, changes in vegetation and land use, construction of dams, and diverse water usage are among the direct human factors that impact water resources (Hosseinpour et al. 2020).

According to previous research, the rise in long-term temperature and variations in precipitation across different regions of the world, along with their effects on the environment, is becoming increasingly apparent. Droughts are the result of various hydrometeorological processes that restrict the availability of surface water or groundwater, leading to conditions that are considerably drier than usual (WMO & GWP 2016; Yimer & Yimam 2021). Hydrological drought characterized by a shortage of streamflow water and a decrease in reservoir and lake water levels (Botai et al. 2021; Kolachian & Saghafian 2021), in contrast to meteorological and agricultural droughts which are mainly responsible for crop failures. Hydrological drought entails a lack of water supply, a reduction in reservoir and groundwater levels, as well as diminished irrigation and hydropower production (Van Loon 2015).

While there are several drought indices available, the reconnaissance drought index (RDI) (Tsakiris et al. 2007) is a widely used indicator for monitoring and measuring the intensity of meteorological drought events, which are caused by the cumulative deficits between precipitation and the atmospheric evaporative demand (Belayneh et al. 2014). So, for meteorological drought, RDI was considered as an effective indicator due to its inherent severity and strong correlation with the actual duration of drought, its inclusion of potential evaporation (PET) as an additional meteorological parameter, in addition to rainfall (Tigkas et al. 2013; Orke & Li 2022) and is especially beneficial in semi-arid and arid areas dealing with climate unpredictability (Tsakiris et al. 2007).

Hydrological drought has significant impacts on both ecosystems and society. These impacts specifically relate to droughts in rivers, lakes, and groundwater. Therefore, it is important to comprehend, develop, and enhance our understanding of hydrological droughts (Van Loon 2015). In Iran, there is high scarcity of data for analyzing streamflow. The majority of reservoirs in the country do not have gauge staff, and groundwater levels are often not measured accurately. In such circumstances, the streamflow drought index (SDI) serves as the most suitable alternative for analyzing hydrological droughts due to its minimal input requirements and straightforward analysis and interpretation (Nalbantis 2008). This study utilizes streamflow data obtained from a donor catchment, as will be defined in the subsequent sections.

Drought assessment enables the monitoring of drought conditions in a cost-effective manner, thereby enhancing the adaptability of water retention. Utilizing satellite-based indices is a prevalent approach for drought monitoring, as it has been demonstrated to be both effective and user-friendly in previous research (Dalezios et al. 2014; Amalo et al. 2018). The normalized difference water index (NDWI) is a remote sensing technique utilized for evaluating water resources (McFeeters 1996) and NDWI serves as a sensitive indicator for detecting changes in water content (Gao 1996).

The evaluation of surface runoff is important for watershed management. Surface runoff data collected during the occurrence of events is important information needed for the planning and design of any watershed project. However, obtaining this data through intensive recording can be costly due to the extensive instrumentation required and the need for long-term monitoring (Sahoo et al. 2020; Panda et al. 2021).

Parishan wetland is located in a semi-arid region characterized by short winters and long, hot, and dry summers. In arid or semi-arid climates, the amount of rainfall is generally low, but its intensity is typically high. Consequently, heavy rainfall in these regions results in rapid hydrological response, leading to higher peak discharges. So, studying the various influential factors that contribute to changes in the water balance of the watershed and its lake has been a significant challenge for researchers. This is particularly true for the Parishan catchment, which is ungauged and remains largely unexplored in terms of its hydrological dynamics.

This problem is termed the ‘prediction in ungauged basin’ (PUB) in hydrology. In order to address the PUB issues, different regionalization methods are commonly employed to replicate the flow of water in ungauged catchments by transferring the model parameters from comparable catchments to those that are ungauged (Samuel et al. 2011; Golian et al. 2021; Wu et al. 2023). The three commonly used techniques for regionalizing parameters include regression-based, physical similarity-based, and spatial proximity-based approaches. Among these, the regression analysis method stands out as the most popular and extensively researched approach (Knight et al. 2012; Yang et al. 2019; Wu et al. 2023).

The primary actions involve creating regression equations that link model parameters to catchment descriptors in the gauged catchment, and determining model parameters in ungauged catchments using the established regression relationship (Boughton & Chiew 2007; Beck et al. 2016; Wu et al. 2023). Nevertheless, some research studies have indicated that the connections between model parameters and catchment descriptors can be intricate, and the process of estimating these parameters in ungauged catchments often results in significant errors (Jafarzadegan et al. 2020; Wu et al. 2023). The physical similarity approach considers that catchments sharing similar physical characteristics (such as climate, vegetation, and topography) exhibit comparable mechanisms of runoff generation and confluence processes (Oudin et al. 2010; Guo et al. 2021; Wu et al. 2023). The donor catchment is chosen based on the spatial distance between the neighboring observed and ungauged catchment, using the spatial proximity method. The parameters of the donor catchment are then transferred to the target catchment (Guiamel & Lee 2020; Wu et al. 2023).

These methods have gained significant popularity in recent years (Reichl et al. 2009; Sellami et al. 2014; Wu et al. 2023). However, the spatial proximity approach proves inadequate when dealing with the substantial spatial differences between neighboring basins (Ly et al. 2013), while the physical similarity approach is constrained by the rationality of selecting catchment characteristics (Heng & Suetsugi 2014). Many catchments face limitations due to geographical or economic factors, resulting in insufficient observed data to effectively calibrate model parameters. Generally, the optimization and calibration of model parameters rely on observed streamflow data at the basin's outlet.

Despite the ongoing dedication and financial investments made to gather hydrometeorological data of high resolution throughout the past century, numerous catchment remain unmonitored without any gauging station. However, the emergence of advanced satellite technologies has addressed this data acquisition gap, enabling the development of hydrological models. These models are designed to effectively simulate both stochastic and deterministic events in time, into the spatial and temporal scales required, ranging from catchment to regional levels (Choudhari et al. 2014; Panda et al. 2021). In this study, due to obtaining sufficient runoff results, both methods, donor catchment and Copernicus satellite data, were used, which the preference of them will discussed in subsequent sections.

Meteorological data, including precipitation, temperature, and evapotranspiration, are crucial factors for hydrological simulation. However, it is equally important to analyze the effects of land use land cover (LULC) changes in hydrological modeling. Therefore, conducting time series LULC mapping and taking into account the impacts of these changes is essential in the watershed simulation process. LULC changes are firstly due to human activities (Foley et al. 2005; Liu & Li 2008) and impact the distribution of water across various hydrological routes including interception, evapotranspiration, infiltration, and runoff (Sterling et al. 2012; Yin et al. 2017). For example, changes in LULC can have an impact on the capacity of surfaces to hold water, alter the speed at which water flows through the environment, and ultimately influence the water balance and hydrological processes in a given region (Zhang et al. 2020; Ware et al. 2024). The alteration in LULC has a significant impact on the runoff behavior of a drainage basin (Getahun & Haj 2015; Sajikumar & Remya 2015; Manderso 2019). Surface runoff is directly influenced by alterations in land use, with the most impactful factors being modifications in forests, agriculture, and settlements.

The emergence of new software and models in hydrology and water resources engineering enables the analysis of stream flows and hydrological simulations at low cost without the need for field studies. It is crucial to conduct hydrological modeling in areas prone to extreme events and natural hazards such as flooding and droughts. By understanding the availability of water resources, stakeholders and policymakers can effectively plan and develop an area. There are several hydrological models available to estimate water resource availability, including lumped models, physically distributed and semi-distributed models, empirical models, and statistical models.

In this study, the Soil and Water Assessment Tool (SWAT) (Arnold et al. 1998) was utilized for rainfall-runoff modeling. SWAT was specifically developed to quantify the runoff and concentration load resulting from distributed precipitation and other meteorological data, taking into account watershed topography, soil, and land use conditions (Rafiei Emam et al. 2017). The SWAT model is extensively utilized in environmental studies due to its accessibility, computational efficiency, and user-friendly interface, which greatly facilitates its application in various hydrologic and environmental issues. This model enables the simulation of water quantity, water quality, and climate change at different temporal scales (annual, monthly, daily, and sub-hourly) over extended periods within a specific basin. Moreover, it can be employed to simulate hydrological processes considering future changes, such as alterations in climate and LULC, and their impact on different variables of the water cycle (Wagner et al. 2017; Schurz et al. 2019; Sanchez-Gomez et al. 2022). Several research studies have been conducted on the SWAT model, including Aghakhani et al. (2019) in the Taleghan watershed, Tajbakhsh et al. (2018) in Zoshk, Nouri et al. (2019) in the Mehrgard watershed of Semirom, Rezaei Moghaddam et al. (2019) in the Lenbaran Chay watershed, and Gorgij et al. (2020) in the Sarbaz watershed. These studies have collectively found that the SWAT model performs satisfactorily in estimating runoff. Furthermore, the GLUE and sequential uncertainty fitting (SUFI-2) algorithms have successfully estimated the model parameters during the validation and calibration stages.

The study conducted by Jaberzadeh et al. (2022) revealed that when simulating runoff in the Dez watershed, Iran using the SWAT and IHACRES models, both models demonstrated good performance in modeling the rainfall-runoff process of the watershed. Nevertheless, the accuracy of the SWAT model was found to be superior to that of the IHACRES model when utilizing the NCEP CFSR climate database. These models have the potential to be utilized in regions with limited access to observational data. Other research aimed to simulate rainfall-runoff in arid and semi-arid regions to aid in decision-making and water resource management, as well as to prevent natural resource degradation. This study utilized the SWAT and IHACRES models, revealing that both models perform well even in areas without a hydrometric station at the output. Furthermore, this study found that the SWAT model outperformed IHACRES in terms of performance (Ahmadi et al. 2019). According to a recent study by Chathuranika et al. (2022), the comparison survey between the SWAT and HEC-HMS models revealed that the SWAT models demonstrated superior performance in seasonal scales. As a result, the SWAT model emerges as an appealing choice for simulating both wet and dry seasonal flows. Also, in a study by Anshuman et al. (2019), the Génie Rural à 4 paramètres Journalie (GR4J) and the Australian Water Balance Model (AWBM) were compared with the semi-distributed SWAT model. The findings revealed that all three models were capable of accurately predicting streamflow. Given their simplicity in structure, minimal data needs, and ease of calibration, it is recommended that conceptual models be favored over SWAT in areas with limited data availability.

Hydrological models are commonly employed to discern the effects of LULC as well as climate changes on runoff (Praskievicz & Chang 2009). These models offer significant frameworks for investigating the changes in different hydrological pathways resulting from both climate and human activities (Leavesley 1994; Jiang et al. 2007; Wang et al. 2010). Distributed hydrological models have been utilized to evaluate the effects of LULC on runoff in water resource management regions (Yang et al. 2008, 2014; Chen et al. 2016). For instance, the hydrological implications of LULC alterations were investigated by Martínez-Retureta et al. (2020) using the SWAT model (Ware et al. 2024). The current study utilized LULC changes over the SWAT model.

Several studies conducted on Parishan lake, for instance in research by Shafiee et al. (2010, 2014), Aghaei-Malekabadi & Khodagholi (2013), and Shafiei & Raeini-Sarjaz (2015), drought conditions have been assessed using PNI, SPI, CZI, MCZI, ZSI, PN, and DI. The present study not only evaluates the impacts of meteorological drought but also considers hydrological drought, as well as examining the relationship between them. Also, the lake area was monitored utilizing remote sensing (Goudarzi-Mehr et al. 2012a, 2012b; Ebadi & Golzar 2016; Khosravian et al. 2017). NDWI can detect the water body of Parishan lake more accurate in comparison with several indices (Khosravian et al. 2017). In the present study, we also used NDWI to monitor the lake but using 112 images to derive lake area changes seasonally throughout the 28 years.

There are some notable innovations in current research; the Parishan watershed is an ungauged watershed, meaning that there is no hydrometry station and only limited data on precipitation and temperature are available. According to the previous studies on Parishan lake and its challenging catchment, there is no research in order to investigate the hydrological process, so we simulated water balance components (WBCs) by using the SWAT model. The main objective of this study is to examine the factors which accelerated the shrinking of Parishan water body, so influences of LULC changes on the catchment's runoff are quantified. Therefore, we calibrated the hydrological model using a regionalization approach, also applying SWAT Calibration Uncertainty Procedures (SWAT-CUP) and SUFI-2 methods, to assess the model's performance and NSE were used for runoff prediction. Furthermore, we conducted an analysis of the model's uncertainty using PBIAS, RSR, p-factor, and r-factor statistics. Additionally, for a better understanding of the watershed condition, metrological and hydrological droughts which are calculated by RDI and SDI are applied to analyze drought conditions during the period 1989–2017. Subsequently, the lake water body was extracted seasonally using NDWI, followed by an analysis of the correlation between drought and changes in the water body. The findings of this study will be of great significance to water resource planners, managers, and extension officials involved in watershed development.

The model requires important details regarding the climate of the watershed, meteorological data, topography data, soil properties, as well as LULC information. The spatial data consists of the digital elevation model (DEM), soil, and LULC maps, which were projected onto UTM zone 39 in the North hemisphere of WGS 1984. The daily meteorological data, specifically rainfall and temperature, from the years 1989 to 2018 were obtained from the Parishan meteorologic station. In cases where data was missing, the weather generator of SWAT was utilized. Since Parishan is an ungauged catchment, we used runoff data from a nearby donor catchment that shares similar characteristics in terms of area, LULC, slope, and soil properties. This data was used for both model calibration and evaluation purposes. The following is a detailed examination of the data used in this research.

Digital elevation data is essential for understanding the spatial distribution of regional topography; it plays a key role in defining the basin's boundary, dividing sub-basin, generating river networks, extracting hydrological parameters and dividing the hydrological response unit (HRU) (Lu et al. 2017). In order to carry out this study, we provide DEM of the Parishan catchment with a resolution of 12.5 m from the Alaska website (https://asf.alaska.edu/). This DEM was used to analyze topographic features and terrain settings.

Since one of the research objectives is to examine the role of LULC changes during the simulation period, its impact on lake drying, regional temperature, and precipitation occurrence has been assessed. According to Ngo et al. (2015), annual surface runoff will increase when forest lands are converted to urban areas and will decrease when forest lands experience significant expansion. On the other hand, by increasing lands with sparse vegetation cover, raindrops directly fall on the soil leading to rapid evaporation and increased losses. Subsequently, a decrease in vegetation cover leads to temperature rise and consequently increased evaporation. Studies show that for every degree increase of temperature in the region, evaporation may increase by 1–3% (Yan et al. 2013; IPCC 2021). LULC for the years 1989, 1994, 1999, 2004, 2009, 2014, and 2017 were extracted and the map of LULC 1989 was introduced to the SWAT model for initial running and creating HRUs. Then, the percentage of changes for each land use relative to the corresponding time step is defined in the land use update section. By this approach, the model can identify the watershed physical condition better, through which we can determine the impact of this factor on calibration and have more accurate consideration in watershed hydrograph changes (Table 1).

Table 1

Percent of LULC changes according to the year 1989

View Large

Meteorological or weather information encompasses the daily highest and lowest air temperatures (Tmax and Tmin), daily rainfall, and solar radiation. Since Parishan is an ungauged catchment, we used runoff data from a nearby donor catchment that shares similar characteristics in terms of area, LULC, slope, and soil properties. In this research daily meteorological data, specifically rainfall and temperature, from the years 1989 to 2018 were obtained from the Parishan meteorologic station. Stream discharge data from Iran's Water Resources Department was used for Fars province and Chamchit hydrometric station. In this study, due to limitations in statistical flow, the data from the website https://cds.climate.copernicus.eu was employed to extract runoff from the output Parishan watershed. Following the research findings, it is recommended whether this approach could be suitable or not for other research considering limited data.

The SWAT model has been used in this study (SWAT version 2012). It is a physically based hydrological model, allowing for the direct modeling of processes such as water and sediment movement, plant growth, and nutrient cycles. These processes are rooted in the concept of water balance (Gassman et al. 2007). The SWAT model was created to assess the influence of management practices in extensive and intricate watersheds encompassing diverse soil types, land uses, and management scenarios over an extended period. This model is semi-distributed, providing the flexibility to analyze data on an annual, monthly, or daily basis (Douglas-Mankin et al. 2010; Janjić & Tadić 2023).

The hydrologic cycle is influenced by the climate and contributes to the water balance through various factors including daily rainfall, maximum and minimum air temperatures, solar radiation, wind speed, and relative humidity. These inputs play an important role in the control of water balance (Wang et al. 2019; Janjić & Tadić 2023). The SWAT model directly extracts data from the files, and in case any information is missing, can generate it. Within the SWAT framework, a watershed is divided into numerous sub-basins, and each sub-basin is further subdivided into multiple HRUs. Each HRU consists of specific combinations of land use, soil properties, slopes, and hydrological components computed for both surface water and groundwater (Neitsch et al. 2011; Yuan & Forshay 2020).

In the given equation, the subscript hru represents the HRU (hydrologic response unit), shows the area of an HRU (L2), represents the fraction of sub-basin area that drains into the sub-basin-scale wetland (i.e., the wetland's catchment area), represents the water surface area of the wetland (L2) at the HRU scale, and other symbols have been previously defined. The SWAT model initially calculates surface runoff and interflow without considering any wetland within an HRU. Afterwards, the equivalent amount of water that would have been generated from the area occupied by the HRU scale wetland (i.e., the fraction of the total wetland area within an HRU) is subtracted from the flows generated across the entire HRU area (Rahman et al. 2016).

HRUs are areas of land that are expected to exhibit similar responses to weather conditions. These HRUs serve as the fundamental units for simulating the land-based aspects of the hydrologic cycle, encompassing processes such as evapotranspiration, infiltration, percolation, and surface runoff. The HRUs involve overlaying soil, land use, and slope classifications within each sub-basin, followed by determining the criteria for defining HRUs. While the implementation of HRU thresholds enhances the computational efficiency of SWAT modeling by reducing the quantity of HRUs, it may also lead to a potential accuracy of watershed characteristics (Gassman et al. 2007; Janjić & Tadić 2023). Then the model can run, initiating the calibration and evaluation processes. SWAT was originally developed and enhanced from a previous model during the early 1990s by Arnold, an individual affiliated with the United States Department of Agriculture (USDA) (Krysanova & Jeffrey 2008; Janjić & Tadić 2023). Different simulated hydrologic processes include the estimation of surface runoff using either the SCS curve number or the Green-Ampt infiltration equation. Percolation is modeled using a combination of layered storage routing technique and a crack flow model. Lateral subsurface flow, groundwater flow to streams from shallow aquifers, and potential evapotranspiration are also considered. The potential evapotranspiration is estimated using Hargreaves, Priestley-Taylor, and Penman-Monteith methods. Additionally, the simulation takes into account snowmelt, transmission losses from streams and water storage, as well as losses from ponds (Arnold et al. 1998; Yuan & Forshay 2020; Janjić & Tadić 2023).

The SWAT model requires process-based input parameters that should be within a realistic uncertainty range. The initial step in calibrating and validating the SWAT model involves identifying the most sensitive parameters for a specific watershed or sub-watershed. The user can determine which variables to adjust based on expert judgment or sensitivity analysis. Sensitivity analysis involves assessing how changes in model inputs (parameters) affect the output of the model. There are two types of sensitivity analysis: local, which involves changing one value at a time; and global, which allows for changes in all parameter values. These two types of analysis can yield different results because the sensitivity of one parameter often depends on the value of other related parameters. The drawback of the one-at-a-time analysis is that the correct values of other parameters are fixed and unknown. On the other hand, global sensitivity analysis requires a large number of simulations. Both types of analysis provide valuable insights into parameter sensitivity and are crucial steps in model calibration (Refsgaard 1997; Janjić & Tadić 2023).

SWAT-CUP is a software tool designed for calibrating SWAT models. It facilitates sensitivity analysis, calibration, validation, and uncertainty analysis of SWAT models (Abbaspour et al. 2007).

The SUFI-2 algorithm (Abbaspour et al. 2007) is chosen for flow calibration in SWAT-CUP software, integrating parameter calibration and uncertainty prediction. It uses uniform distributions to represent input parameter uncertainty and calculates model output uncertainty as 95% prediction uncertainty (95PPU) from the cumulative distribution of outputs via Latin hypercube sampling (Abbaspour et al. 2004). Global sensitivity analysis uses Latin hypercube regression analysis to rank parameters’ sensitivity to streamflow (Abbaspour 2015). The T-stat and P-value measure sensitivity, where a higher absolute T-stat and lower P-value (<0.05) indicate significant sensitivity (Thavhana et al. 2018).

In this study, by applying SWAT-CUP and SUFI-2 algorithms, the effectiveness and sensitivity of one parameter relative to another based on P-value and T-stat indices are investigated; then, the most influential ones on the stream flow were examined. Regarding P-value, if the value be closer to zero, the sensitivity of the parameter becomes higher. Similarly, for parameter T-stat when its absolute value being greater than one indicates a higher sensitivity of the parameter (Abbaspour et al. 2017).

We utilized 30 parameters, 17 of which were found to be the most effective as indicated in Table 3.

Calibration involves the process of enhancing the model's parameterization based on specific local conditions, thereby minimizing prediction uncertainty. This is achieved by selecting values for the model's input parameters within their respective uncertainty ranges and comparing the model's predictions under assumed conditions with observed data for those conditions (Refsgaard 1997; Janjić & Tadić 2023). To minimize uncertainty, the majority of calibrations were performed using the automated tools SWAT-CUP and SUFI-2. In SUFI-2, all uncertainties, such as those in input data, conceptual model, and parameters, are reflected in the model output.

The level of uncertainty in SUFI-2 is quantified using two efficiency criteria: the p-factor and the r-factor. The p-factor measures the percentage of observed data within the 95% prediction uncertainty (95PPU), ideally being 1, indicating 100% coverage. The r-factor reflects the thickness of the 95PPU band, calculated as the average distance between the upper and lower 95PPU divided by the standard deviation of observed data, ideally approaching zero for best accuracy (Yang et al. 2008). However, achieving a balance between these two factors is often necessary, across various conditions such as average, wet and dry, high and low temperature years. These calibrations utilized the same parameter ranges and iteration numbers. Manual calibration, on the other hand, was infrequently carried out.

SUFI-2 quantifies optimal parameter values through an iterative process aimed at optimizing a selected objective function. The key steps in the SUFI-2 calibration procedure include:

(1) Selection of Objective Functions: Choose 17 factors from different objective functions in SUFI-2. (2) Definition of Calibration Parameters: Identify the parameters to be calibrated and define their minimum and maximum ranges, assuming uniform distribution within the calibrated sub-basins. (3) Random Parameter Combination: Generate n combinations of random parameter values using Latin hypercube sampling (McKay et al. 1979), and run the model n times for each parameter set. (4) Calculation of Objective Functions: Compute the objective functions for each of the n simulations. (5) Post-Processing: Analyze the n simulations to determine the best-fit parameter sets for each sub-basin by comparing simulated and observed data. (6) Calculation of p and r factors: Determine the p-factor and r-factor. (7) Further Rounds of Sampling and Model Runs: Conduct additional rounds of sampling and model runs with updated parameter ranges to address large uncertainties in the initial parameters.

This process ensures accurate model calibration, accounting for uncertainties and improving the fit between simulated and observed data (McKay et al. 1979). In this study, we divided the streamflow into two sections and performed calibration using the SUFI-2 method with different objective functions. Following Abbaspour et al. (2007), we executed calibration and selection of parameter ranges over 100 iterations. The final calibrated parameter set, referred to as the near-optimal parameter set, was derived from 17 out of 30 parameters, which best reproduced the behavior and characteristics of the flow.

The performance of the SWAT model was assessed using three indices: the correlation coefficient (⁠⁠), Nash–Sutcliffe efficiency coefficient (NSE), and percent bias (PBIAS) (Abbaspour 2015; Lu et al. 2017). Objective functions, such as root mean squared error (RMSE) and NSE (Nash & Sutcliffe 1970), are commonly employed to minimize the difference between the observed and simulated flows. The calibration approach has been extensively used in various hydrologic models (Sorooshian 1991; Gan & Biftu 1996; Gupta et al. 1998; Gupta et al. 2009). Furthermore, evaluation is a customary procedure in hydrological modeling (Andréassian & Kaloustian 2009) employed to assess the model's performance using previously unprocessed data (Kim et al. 2018). Because there is no hydrometry station at the outlet of our study area and we cannot calibrate the relevant parameters directly, we used a regionalization approach to transfer the relevant parameters from the nearby catchment to our study area. To do that, the Chamchit catchment, 7 away from the Parishan watershed, with the observed discharge data was selected. Data was available at the daily scale from 1989 to 2017 in this catchment. The characteristics of the two basins are similar with a difference in annual precipitation of 2.6%. More physical characteristics are presented in Table 4.

For calibration, the data of donor catchment (Chamchit) were used from the year 1992 to 2013, and in evaluation, we utilized 2014–2017 data.

Sun et al. (2014) introduced a novel technique for calibrating hydrological parameters in ungauged basins. Rather than using river discharge data, they used river water level data obtained from satellite radar altimetry observations at the basin outlet. In a different study, Sun et al. (2010) utilized river flow width data from satellite observations as a proxy for river discharge records to calibrate a hydrological model (Rafiei Emam et al. 2017). The regionalization method allows for the prediction of hydrological model parameters in ungauged basins by leveraging the similarities in characteristics among watersheds. This approach, as discussed by various researchers (Vandewiele & Atlabachew 1995; Bárdossy 2007), involves transferring hydrological parameters from similar watersheds, a process known as regionalization. Factors such as soil composition, topography, land use, precipitation, and temperature play a crucial role in this transformation, as highlighted in the study by Samaniego & Bárdossy (2005). Different techniques, including regression methods (Bastola et al. 2008), spatial proximity (Merz & Blöschla 2004), and physical similarity (McIntyre et al. 2005), are employed for this purpose. Additionally, geostatistical methods such as Kriging have been utilized by Parajka et al. (2005) to transfer model parameters effectively.

The Copernicus Data Space Ecosystem provides access to vast quantities of open and free Earth observation data. This encompasses the Copernicus Sentinels 1,2,3,5P Missions, Copernicus Contributing Missions, Federated datasets, and Complementary data. The Copernicus Sentinels Missions specifically provide access to Earth observation data obtained from the Copernicus Sentinel satellites. The Copernicus Contributing Missions, on the other hand, offer data that complements the Copernicus Sentinel Missions and play a vital role in Earth observation. These missions are carried out by ESA, its Member States, and international third-party operators, and they provide high-resolution optical and radar data. In addition to the extensive online Sentinel data archive and access to the Copernicus Contributing Missions, more datasets will be incorporated in a federated manner over time. This integration allows for the utilization of other European back-ends within the Copernicus Data Space Ecosystem, enabling further advancements. Furthermore, the website also offers a wide range of supplementary data sources, including high-resolution satellite imagery from various providers and data offerings from different Copernicus services. This ensures a comprehensive and diverse dataset for your specific requirements (Copernicus Data Space Ecosystem 2014).

The observed data on the i-th day is denoted as ⁠, while the predicted/simulated value on the i-th day is represented as S_i. The mean of the measured/observed data is denoted as ⁠, and the mean of the predicted data is represented as ⁠. The total number of data is denoted as n. To assess the uncertainties associated with the SWAT model, the p-factor was used. The p-factor quantifies the percentage of measured data that falls within the 95% prediction boundary. It ranges from 0 to 1, with values closer to 1 indicating a higher model performance and efficiency. On the other hand, the r-factor measures the average width of the 95% PPU band divided by the standard deviation of the observed variable. It varies from 0 to ∞, indicating the variability of the model's performance. Both the p-factor and the r-factor are mathematically expressed in Equations (5) and (6), as presented by Abbaspour et al. (2007a, 2007b) and Yang et al. (2008, 2015).

The RDI, as proposed by Tsakiris et al. (2007), has been widely employed for the assessment of meteorological droughts (Fooladi et al. 2021). By utilizing precipitation and evapotranspiration data, it enables the quantification of drought severity and frequency (Barker et al. 2016; Abbas et al. 2021). The RDI can be calculated in three different forms: the alpha RDI (Equation (9)), normalized RDI (Equation (10)), and standard RDI (Equation (11)) (Tsakiris et al. 2008). The standard RDI and alpha RDI are specifically used for assessing drought severity and aridity classification, respectively. The computation of RDI for any given location starts with an aggregated form (i.e., alpha RDI) using a monthly time step and can be computed for each month of the hydrological year or for an entire year. The assessment of meteorological droughts was conducted using the RDI. In this research, a time scale of 12 months was employed and the average RDI values were utilized to identify the meteorological dry and wet periods in the study area.

In which is equal to the ln ⁠, while is its average and is its standard deviation, respectively.

A drought event occurs when the RDI values remain negative and below the threshold of −1.0 for more than two consecutive months. The event concludes once the indices turn positive (McKee et al. 1993; Schwalm et al. 2017). The severity of the drought is indicated by more negative values. Each index's results help in characterizing and categorizing the drought event. The classification varies based on the drought's impact on specific systems and factors related to the study objectives and location. The indices used in drought analysis are categorized as follows (Table 5): extremely wet conditions (≥2); severely wet (1.5–1.99); moderately wet (1.0–1.49); near-normal (0.0–0.99); mild drought (0.0 to −0.99); moderate drought (−1.0 to −1.49); severe drought (−1.5 to −1.99); and extreme drought (≤− 2.0) (Tigkas et al. 2013).

It is crucial to establish the reference flow levels and indicators of drought severity in order to effectively develop and manage water resources. These indicators help determine the severity of drought based on factors such as duration and flow deficit. Various indices are used to measure the deviation of streamflow from historical norms for a specific period. Water supply planners often rely on these indices to make informed decisions (Cancelliere et al. 2003; Wu et al. 2018). Streamflow, being the key variable in assessing the availability of surface water resources, is directly associated with hydrological drought occurrence under normal conditions (Van Lanen et al. 2013; Hasan et al. 2019).

and S_k represent the mean and standard deviation of cumulative streamflow volumes during the reference period k, calculated over an extended period. The SDI and RDI range from −2 to +2, indicating the spectrum of wetness and dryness. Values below −2 and above are considered extremely dry and wet, respectively. Nalbantis & Tsakiris (2009) provide a hydrological drought classification based on RDI and SDI as shown in Table 6.

The relationship between RDI (3-, 6-, 9-, and 12-month periods) and SDI (12-month period) is established using linear regression analysis of historical data. Forecasted values and their probable ranges are calculated based on confidence levels. Additionally, the correlation between hydrological and meteorological indices for 3- to 12-month time series was analyzed with and without delays of 1–5 months (Tigkas et al. 2012).

The index is specifically created to optimize the reflection of water using green wavelengths, represented by the green band, such as TM band 2. It aims to reduce the low reflection of water features in the near infrared band, represented by NIR, such as TM band 4. Additionally, it utilizes the high reflection of NIR by vegetation and soil features. Consequently, water features are enhanced and have positive values, while vegetation and soil typically have zero or negative values and are therefore suppressed (McFeeters 1996). NDWI ranges from −1 to +1. This product allows for water bodies to be observed and can be beneficial in evaluating the effects of flooding. The total number of pixels exceeding the threshold (zero is as threshold) are tallied to depict the water body, and subsequently, the surface area is computed. In this research, the area of lakes is determined by summing the positive digital number values (with specific thresholds for each map) and then multiplying them by the resolution number (based on Landsat image and corresponding band). The outcome of this calculation will indicate the area of lakes.

The ENVI (Environmental for Visualizing Images) product was utilized for processing and analyzing geospatial imagery obtained from satellite images extracted seasonally between 1988 and 2017 via the Global Visualization Viewer Glovis (Glovis 2017). The images were sourced from Landsat 4–5 Thematic Mapper I, Landsat 7 Enhanced Thematic Mapper plus (ETM+), Landsat 8 Operational Land Imager (OLI), and Thermal Infrared Sensor (TIRS). Before initiating image processing for monitoring lakes area, certain preprocessing procedures need to be carried out in accordance with the images. When sensors detect reflective waves, three types of errors may arise in the raw images. The first type is Sensor's error, the second type is Geometric errors, and the third type is linked to the influence of atmospheric agents on the sensor reflection gain. Preprocessing endeavors to reduce these effects to the required level for a particular application (Hadjimitsis et al. 2010).

The Chamchit hydrometric station was selected for the examination of the SWAT model results and analyzing sensitivity of flow parameters effectiveness. This is necessary because the model has not been calibrated precisely to describe these processes. The use of SUFI-2 algorithms in SWAT-CUP software and the identification of effective parameters are discussed in Section 2.3.1. As the study area includes changes in elevation and due to its being a lake in the Perishan watershed, it requires the selection of parameters that can accurately simulate the model output.

Therefore, considering that the topographic slope can induce significant changes in precipitation, snow, ice formation, and melting, in this study, following Abbaspour et al. (2017), snow parameters (SMTMP and SFTMP) were applied separately from others and parameters that play a significant role in the formation of surface runoff and base flow in the model were selected as illustrated in Table 3. Another sensitive parameter in simulating runoff is CN2, the curve number. The runoff curve number estimates runoff based on the relationship between precipitation, soil hydrological group, and land use. Other researchers have also found that the curve number SCS is the most sensitive parameter in hydrological modeling studies (Qiu et al. 2012; Abbaspour et al. 2015; Gyamfi et al. 2016). Other sensitive parameters include the value of Manning's N for the main channel (CH_N2), effective hydraulic conductivity in the main channel floodplain (CH_K2), coefficient of infiltration into deep groundwater (GW_REVAP), base flow alpha coefficient (ALPHA_BF), precipitation impact coefficient (RFINC), fraction of sub-basin area draining to ponds (PND_FR), and volume of water stored in ponds during emergency overflow (PND_EVOL). These parameters have had a greater impact on the calibration process. Although other parameters have been used in this study, those mentioned have been highlighted due to their significant influence on the calibration process. Additionally, other parameters were not considered sensitive to watershed stream flow because their P-values were more than 5% of the research (Asres & Awulachew 2011; Betrie et al. 2011; Gamvroudis et al. 2015). This is expected because conditions such as LULC, soil characteristics, and climatic factors vary from one watershed to another. Model output calibration based on criteria ⁠, NSE, PBIAS, r-factor, and p-factor calibration distance for flow at selected stations in the watershed is shown according to Table 7.

In the SUFI-2 program, parameter sensitivity is determined using the p-factor and r-factor (Abbaspour et al. 2007). Table 7 presents the criteria for the calibration period at selected watershed stations. For the calibration period from 1992 to 2013, the extracted range values for observed data and simulation estimate the R² and NSE ranges as 0.81 and 0.67, respectively. According to the scoring index in Table 8, this indicates a good calibration range. These results are also acceptable for the evaluation period.

In this study, parameter evaluation for basin discharge uses the p-factor and r-factor to assess uncertainty in the conceptual model. A p-factor close to 1 indicates high model performance, while an r-factor less than 1 suggests good calibration. After calibration and evaluation, p-factor values were 0.06 and 0.02, and r-factor values were 0.02 and 0, respectively. These results indicate effective runoff simulation and reduced data uncertainty, confirming the SUFI-2 algorithm's potential in SWAT-CUP for watershed runoff simulations and aligning with findings from other studies (Tang et al. 2012; Rostamian et al. 2013; Singh et al. 2013).

The PBIAS factor in the calibration extracted is the value of 24 which is desirable according to the previous studies. The main reason for the high value of PBIAS in this study can be attributed to utilizing statistics from the donor station. In terms of evaluation, the statistics from 2014 to 2017 also demonstrated that the model results were within a good range. Performance rating of model simulation (Congalton 1991; Moriasi et al. 2015; Leta et al. 2021; Houshmand Kouchi et al. 2017) is illustrated in Table 8.

As illustrated in the figure, the extracted information from the Copernicus site has shown adequate accuracy in estimating the minimum and normal flow values. However, it's noteworthy that the extraction of the maximum hydrograph has been the main concern. The primary reason for this lies in the extraction of runoff amounts based on images and moisture levels that have been more prevalent within that time frame. Therefore, it could not accurately predict the maximum runoff.

After conducting the accuracy assessment of runoff estimation by the site, it is now time to proceed with the calibration of the output SWAT in the software SWAT-CUP using the extracted parameters and the estimated runoff level by the Copernicus site. As observed in Table 9, within the calibration range, considering the values of the range NSE and ⁠, which are respectively 0.5 and 0.62, the model results have been acceptable. However, in the evaluation range, considering the values of the range NSE and ⁠, which are respectively 0.48 and 0.71, a good correlation value has been estimated, but the NSE coefficient indicates an acceptable range. Therefore, based on the evaluation index, which represents the assessment of the flow hydrograph, it can be concluded that using statistical information from processed satellite data of the Copernicus site is acceptable for some areas with limited or absence of runoff information, considering a percentage of possible error.

There are many hydrological drought indices, but except for SDI, all indices are more data-intensive (Nalbantis & Tsakiris 2009; Sutanto et al. 2020). Therefore, in this study, the historical hydrological drought trend in the Parishan watershed was analyzed using SDI. SDI is a point or site drought indicator that gives information about the stream temporal and special variation. The analysis was computed using the DrinC model (Drought Indicator Calculator), which is developed to determine SPI, RDI, and SDI using monthly input data. Application of DrinC is used which is widely used in several places, especially in arid and semi-arid regions, show that it is being achieved as a useful research and operational tool for drought analysis (Tigkas et al. 2013, 2015).

The RDI displays multiple fluctuations, with an extended period of wet and normal conditions from the initial stages up to 2006, except for the moderate drought in 2001. Subsequently, from 2006 until the end of the statistical study period, the RDI indicates a decreasing trend. Between 1992 and 1999, the watershed faced moderate wet and extremely wet conditions. The situation transitioned from normal to severe drought in 2006, lasting until 2011, and returning to normal conditions by 2012. The remaining years experienced fluctuations ranging from normal to moderate and severely drought conditions.

The SDI exhibits a decreasing trend and aligns with the RDI trend. However, there is a slight variation in 1992, where the RDI indicates a moderately wet condition while the SDI suggests that the watershed experiences non-drought conditions. The following year, the SDI shows a moderate drought, and from 1994 to 1998, the condition returns to non-drought. In 1997, there is a discrepancy between the RDI and SDI, the RDI indicating normal conditions while the SDI shows a significant peak, representing the best hydrological condition during the study period. Another difference occurs in hydrological drought, where the RDI indicates moderate wetness while the SDI suggests a mid-level drought. However, in the subsequent years, the condition returns to non-drought. Additionally, in 2002, the RDI shows moderate drought while the SDI indicates a non-drought condition. From 2002 to 2017, the SDI and RDI trends are similar, with only a slight difference in 2010, where the RDI shows extreme drought while the SDI suggests a mid-level drought.

The meteorological and hydrological droughts of the basin were assessed using the Hargreaves method, with correlations between RDI and SDI examined both with and without delays of 1–5 months. The results indicated a higher correlation between RDI and SDI in the no-delay state. Overall, the annual correlation between RDI and SDI was weak (R² = 0.3), and there was no significant correlation in the 3-, 6-, and 9-month time series. This suggests that reduced rainfall alone does not significantly impact streamflow reductions. The region's agricultural activities, including changes in vegetation cover and cropping patterns, also influence flow levels. This is further supported by the weak correlation between lake area and drought.

In the concluding part of this research, we utilized a calibrated SWAT model to investigate the water balance and analyzed the impacts of LULC and climate change on hydrological processes in the watershed. Therefore, through a holistic assessment, we examined several hydrological components such as SurQ (surface runoff generated in watershed), GwQ (groundwater contribution to stream in watershed), LatQ (lateral flow contribution to streamflow in watershed), WYLD (water yield to streamflow from HRUs in watershed), SW (amount of water stored in soil profile in watershed), and ET (evapotranspiration), as outlined in Table 10. Overall, these components provide a comprehensive understanding of the hydrological balance and potential changes in the watershed dynamics. It is worth mentioning that the basis of this comparison was the establishment of a baseline scenario, upon which potential changes resulting from different land uses and baseline climatic conditions were evaluated. This research approach allowed us to separate the actual climatic effects on watershed hydrology using the SWAT model and present a clearer picture of potential vulnerabilities and changes in hydrological processes, which can solely be attributed to climate variations (Berberoglu et al. 2020; Debnath et al. 2022).

As seen in Table 10, the average values of hydrological components were estimated relative to the total length of the statistical simulation period, which in the SWAT model was defined to reflect the impact of land use changes over a time series interval of 5 years, to illustrate the effect of these changes on drying of Lake Parishan over the entire statistical period. As observed, Lake Parishan had a significant impact on the water balance of the watershed, such that we witnessed high levels of moisture and evaporation from 1994 to 2010, followed by a reduction in the watershed's evaporation thereafter, which is primarily due to the absence of the lake. These analyses reveal similar patterns where temperature increase, reduced rainfall, changes in agricultural and urban land use in the watershed have led to increased water resource extraction, affecting on drying of the lake. The values of ET in this watershed indicate the positive impact of the lake on the watershed's climate, which has gradually decreased as the lake level has decreased. Furthermore, SW, SurQ, WYLD, GwQ, and LatQ values, with their decreasing trends, indicate a significant impact on vegetation cover, temperature, and precipitation changes on the hydrological components of the studied region.

The presence of wetlands and lakes is crucial due to their role in meeting human needs and providing essential ecosystem services. The phenomenon of lake drying can be considered a warning sign indicating the disruption of climatic balance. Therefore, precise understanding of hydrological variables and assessment of drought conditions are influential in identifying the causes of this issue.

Based on previous research, the current study takes into account certain differences. For instance, in previous studies, land use change and climate variations were often examined separately. However, in this research, an attempt has been made to bridge this gap by integrating various studies and utilizing the SWAT model. The goal is to offer new insights and a more comprehensive understanding of the combined impacts of land use changes and climate variability on soil and water resources. This approach is crucial as it reveals that land use changes and climatic variations are interconnected phenomena that can profoundly influence each other.

This study employs a two-step approach for simulating runoff in ungauged areas. First, a regionalization method is used, focusing on correlations between nearby stations based on spatial proximity and physical similarity. Second, Copernicus satellite data is integrated to enhance findings and reduce uncertainty by comparing with data from donor catchments. Additionally, incorporating LULC changes at various time steps improves accuracy, as variations in built-up areas, vegetation, and lake areas affect the hydrological processes in the watershed. Unlike studies with limited LULC change analysis, this research uses LULC maps at 5-year intervals to better plan for sustainable water resource development in the catchment.

This study aimed to investigate the impact of land use changes and climatic conditions in the Parishan watershed on the drying of the lake. The results indicated that the SWAT model can sufficiently and acceptably simulate the watershed runoff. Furthermore, in the absence of sufficient statistics and hydrometric stations, Copernicus satellite data can be examined and analyzed following the principles of station analysis and sampling. The conducted investigation showed that changes in temperature, precipitation, and vegetation cover had a significant impact on the lake level, especially on the hydrological balance components of the watershed.

Some studies have aimed to identify key factors contributing to the drying of Parishan lake. For instance, Ghazali (2012) highlighted the presence of underground connections between neighboring aquifers and the lake. It was observed that groundwater extraction doubled from 1992 to 2009, increasing from 9.5 to 43.6 ⁠, leading to a decline in the groundwater level. The impact of unregulated agricultural well withdrawals on the lake is significant. However, assuming all withdrawals be from the lake is not justified, as withdrawal patterns vary, with consistent withdrawals during planting seasons. Despite water withdrawals and low aquifer recharge in some years, the lake shows signs of restoration, indicating a trend toward increased water abundance.

Water balance changes, particularly during dry seasons, reveal that Parishan lake receives water not only from its own catchment but also from adjacent regions. Fluctuations in the lake's water levels are influenced by precipitation variations in neighboring watersheds, unregulated agricultural withdrawals, and aquifer depletion. Therefore, any human activities around the lake or in neighboring watersheds can impact the influx of springs that are crucial for the lake's water supply. Additionally, the Nargesi Dam in the neighboring watershed should be considered in water management, as changes in adjacent watersheds can affect the lake's water levels, making it important to secure water rights from the dam.

Given the severe drought conditions since 2007, the study finds an indirect relationship between drought and the reduction in the lake's area. However, changes in land use, lake surface area, and streamflow indicate that drought alone does not account for the lake's dryness. The study highlights the challenge of limited data in ungauged watersheds and advocates for using satellite data as a solution. Incorporating Copernicus satellite data improves the understanding of catchment characteristics and enhances result accuracy. The study also references Isaei & Isaei (2015), who proposed a lake restoration strategy involving a diversion channel. The current study's runoff volume calculations can aid in the efficient planning of this diversion channel for the Parishan region.

Based on the results obtained in this study, the following recommendations are proposed to assist restoration of the lake and efficient management:

• 1. Investigating hydrological relations of Lake Parishan with surrounding areas: Further research should be conducted to explore the hydrological connections between Lake Parishan and its adjacent regions. This involves studying groundwater flow patterns, surface water interactions, and the impact of land use changes on hydrological processes in the vicinity of the lake.

• 2. Controlling water withdrawals from Lake Parishan and its surrounding areas: Implementing measures to regulate water withdrawals from Lake Parishan and its adjacent areas is essential for maintaining the lake's water balance. This may include regulations on agricultural water usage, monitoring groundwater extraction, and enforcing water conservation practices.

• 3. Determining the lake's water rights: Establishing clear water rights for Lake Parishan is crucial for sustainable water management and preventing over-exploitation of its resources. Legal and policy measures are needed to allocate water resources fairly and safeguard the lake's ecosystem.

• 4. Extending management practices to adjacent areas: Given the significant influence of neighboring regions on Lake Parishan's water levels, it is important to extend water management practices beyond the lake's watershed. This could involve implementing similar water conservation measures and sustainable land use practices in surrounding areas to mitigate the impact of land use changes on the lake.

In conclusion, addressing such crises requires multifaceted management approaches based on practical decision-making and outcomes.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.