Simulation of ecohydrological processes influencing water supplies in the Tuul River watershed of Mongolia

Maintaining reliable water supplies is a significant policy issue facing many countries. The uncertainty in the global climate aggravates these problems. Growing populations can accelerate the water demand, while cost-effective, supply-side options remain limited (Dharmaratna & Harris 2012). It is essential to have sufficient water for multiple uses in the watershed to have sustainable water resources. Water resources are complex systems influenced by many interacting factors like terrain, precipitation, humidity, air temperature, soil and vegetation type, land use, and land cover. Understanding the current ecohydrologic processes for a watershed is essential to assess future climate and its potential impacts on hydrologic regimes (Hongfu et al. 2012).

Many researchers have used different methods to assess water resources, from simple water balance equations to complex hydrological models incorporating various water resource system components. The development of accurate and informative integrated watershed models that help water managers better understand the issues within their basin currently and in the future is essential. Successful watershed models can lead to the development of timely projects for securing future water resources. Watershed modeling is a standard tool to understand a catchment's behavior under dynamic processes (Fiseha et al. 2013; Tsvetkova & Randhir 2019). The best model fit aims to simulate reality with minimal parameters and reduced model complexity. Watershed simulation using various parameters and equations helps predict impacts of management practices, climate, and land use (Marshall & Randhir 2008; Devia et al. 2015; Talib & Randhir 2017). Wang et al. (2021a, 2021b) used an ensemble data-driven hydrologic model to study climate boundary conditions in 16 watersheds in China. Lindenschmidt & Rokaya (2019) used a stochastic modeling approach to study ice jam flooding. Wang et al. (2021a, 2021b) used a clustered polynomial chaos expansion model to explore uncertainty in hydrologic prediction in the Ruihe Watershed in China.

This study simulates ecohydrological processes in a watershed with an extreme continental climate. The Soil Water Assessment Tool (SWAT) model is used to simulate the hydrologic processes of the Tuul River (TR) watershed and how those hydrological processes influence the water resources. The model is widely used for the planning and management of watersheds throughout the world. The model components include hydrology, land use, weather, crop growth, erosion/sedimentation, nutrients, pesticides, and agricultural management (Arnold et al. 1998). The SWAT has several advantages: (i) can simulate small to large complex watersheds (Arnold et al. 2012); (ii) has potential for use in ungauged or watersheds with limited monitoring data (Gassman et al. 2007); (iii) performs well in studies worldwide; and (iv) capable of predicting changes to water quality and quantity (Douglas-Mankin et al. 2010). This study is unique in applying the SWAT model to assess spatial and temporal hydrologic processes under semiarid, cold conditions in mixed-use, large river basins. The TR Basin is also less studied in literature at a watershed scale.

This study aims to simulate processes in the TR watershed and subbasins and analyze the influences of those processes on the water resources. Specifically, to: (i) model ecohydrological processes in the study watershed; (ii) analyze temporal and spatial variations in hydrologic flows; and (iii) evaluate factors influencing water supply at subbasin and watershed scales.

Alternate hypotheses are: (i) hydrological models can simulate watershed processes to reasonable accuracy; (ii) temporal and spatial variations of water flows are significant in the watershed; and (iii) specific subbasin processes impact local water supplies.

The TR watershed is in Mongolia's central and northern parts (Figure 1). The watershed covers 49,994 km² and 704 km long. It flows from east to west, covering central parts of Mongolia, and is the freshwater source for its capital, Ulaanbaatar. The TR watershed has 69% rainwater, 25% groundwater, and 6% snowmelt (Dolgorsuren et al. 2012; Dalai et al. 2019).

The TR flows down to Ulaanbaatar from the mountain taiga and forest-steppe region. Steppes cover about 80% of the watershed. The upper river basin consists of steep rocks and forests with a valley width of approximately 1–3 km in width. Downstream of the city of Ulaanbaatar, the valley becomes broad, and the valley width expands to 8–10 km. The geographical coordinates of the point considered as the TR's origin are 108°13′20″ E, 48°30′39″ N, and the coordinates of the river confluence point are 104°47′52″, 48°56′55″ (Dolgorsuren et al. 2012).

The TR watershed is of high elevation surrounded by mountains with elevation ranging between 2,792 and 773 m (Figure 2). The climate is influenced by differences in day and night temperatures, long winters, and short summers. Most of the precipitation is in the summer months, dominated by warm, dry air and frequent thunderstorms. A sudden cold is observed in late August or early September, decreasing precipitation in the autumn months (Dolgorsuren et al. 2012).

The TR watershed has a continental climate with a wide range of temperatures, low humidity; varying precipitation; cold and long winter; and warm summer months. The rainy season (June–August) is when 90% of precipitation is received. The average annual flow in the TR at the Ulaanbaatar gaging station is 25.6 m³/s (Sukhbaatar et al. 2017). On the other hand, the river starts to freeze in early October until mid-April (Dolgorsuren et al. 2012), with a maximum ice depth of 1.16 m in mid-February (Sukhbaatar et al. 2017). Sukhbaatar et al. (2017) used the SWAT to model hydrological processes and climate change impact in the Upper TR watershed. They found that there have been cyclic variabilities such as several wet and dry periods. The latest dry period was observed during 1996–2015, where the average air temperature increased and the precipitation pattern reduced.

The average air temperature in the study area is −0.38 °C. The average temperature is from a minimum of 39.0 °C in January to a maximum of +37.20 °C in July. Warming trends are increasing by 2.0 °C from the norm in the TR watershed as observed in meteorological stations. With the increase in the number of hot days and highest temperatures since 1940 concentrating in the last decade, climate change impacts are increasing in the watershed (Sukhbaatar et al. 2017).

The TR watershed has an annual precipitation range from 253 to 275 mm (Sukhbaatar et al. 2017). Roughly 85–90% of the total precipitation is during the growing season (Dolgorsuren et al. 2012). In the TR watershed, the total rainfall does not change much, but rain duration is decreased during the summer. The aridity in the watershed is attributed to climate change, increasing evaporation (E), and causing aridity in the basin. With deficient water supplies to support plant cover, most areas have a low vegetative cover. The difference between evaporation and precipitation (E–P) is trending to be substantial in the last three decades.

Mountain soils evenly dominate more than 56% of the TR watershed area (Figure 3). The mountain soils include mountain dernotaiga, soddy taiga, forest dark, mountain meadow, and mountain dark chestnut soil classes. Along the river, meadow, swamp cryomorphic, meadow cryomorphic, and meadow solonchak soils are distributed. The soil types include mountain soil (56.3%), soils in steppe valley (26.2%), low hills soil (8.3%), humid region soils (6.6%), and the remaining area covered by other soils. Agriculture occupies 91.2% of the watershed, forest covers by 6.8%, and the remaining regions are settlements (Figure 3). According to the unified land classification of the Mongolian Law of Land, individual needs include agricultural land, urban and local settlement area, roads, forests, and water (Dolgorsuren et al. 2012).

Issues: The TR, once known for its clean and pristine water, is environmentally degraded in the last several decades. The river is the primary source (94%) of recharge to the alluvial aquifer (Tsujimura et al. 2013). It is the primary source of drinking water to the city of Ulaanbaatar, the capital of Mongolia. Therefore, water sources and the river flow are the most important to the city's estimate 2021 population of 1.6 million residents (UN 2019). Over the past two decades, the river is facing low flows and declining groundwater levels. The TR watershed faces some critical issues regionally. The Upper TR watershed is forested and has issues regarding deforestation due to the expansion of tourism and camping. The Middle TR watershed downstream of the Ulaanbaatar city is highly polluted with poorly treated wastewater discharged from Waste Water Treatment Plant. In addition, there have been issues regarding the overgrazing of livestock and pollution from the Zaamar gold mine in the lower TR watershed (Dolgorsuren et al. 2012). Several studies have assessed the water resources in the headwaters (Sukhbaatar et al. 2017) and groundwater abstraction (JICA 2010; Tsujimura et al. 2013). However, subbasin and river basin-level analyses were not conducted for watershed management.

A methodological framework (Figure 4) shows hydrological components (runoff, infiltration, subsurface, groundwater, and evapotranspiration) used in this study. The framework includes data collection, pre-processing, calibration, validation, model performance, and simulation of the hydrological components.

Hydrological models are commonly used to understand a watershed system and its hydrological processes better. The SWAT is a comprehensive model to provide a detailed analysis of water resources and easy to modify any changes to the system (Arnold et al. 1998). The SWAT model has been increasingly and successfully used worldwide for various purposes, e.g., assessing water resources, planning, and management, studying hydrologic impacts of climate or land-use changes (Ficklin et al. 2013) from global to local scales. The SWAT is an ecohydrologic model (Devia et al. 2015) useful in studying the effect of management and climate on water resources. The model segments the simulation area into subwatersheds and hydrological response units (HRUs) that incorporate unique land-use, soil attributes (Abouabdillah et al. 2014). It uses daily meteorological inputs to replicate the watershed mechanisms; the output could be summarized monthly or yearly for an extended simulation period. The SWAT model is also commonly used to simulate ungauged watersheds as there are many countries where hydrological and meteorological monitoring is limited.

The SWAT model has vast applications and is widely used worldwide for the assessment of water resources (Dessu et al. 2014; Jujnovsky et al. 2017), water quality (Abbaspour et al. 2015; Fan & Shibata 2015), land-use changes (Wang et al. 2014; Krysanova & Srinivasan 2015; Neupane & Kumar 2015; Sajikumar & Remya 2015), and climate change (Johannsen et al. 2016; Faramarzi et al. 2017).

Although the SWAT can require many inputs to simulate the hydrological process, it has been successfully applied in data-limited studies and is a feasible and valuable approach for strengthening water resources management (Niu et al. 2014; Fukunaga et al. 2015).

The following hydrological components were used to simulate the hydrology of the study watershed. The water balance equation is represented by ⁠, in which and are soil water at the beginning and end of the day i, t is the time interval in days, and R, Q, ET, P, and QR represent precipitation, runoff, ET, percolation, and return flow, respectively. The SCS curve number equation, ⁠, is used for predicting surface runoff (USDA-SCS 1972). Where Q, R, and s are runoff (mm), rainfall (mm), and a retention parameter that varies (1) with watershed because of land use, soil, management, and slope all range and (2) over time variation of soil water content. The retention parameter is associated with the curve number (CN) in the SCS formula (Arnold et al. 1998). The percolation is modeled as a storage routing technique, integrated with a crack-flow model, to simulate flows through each soil layer. Percolating water enters groundwater as a return flow downstream. The storage routing is based on in which SW₀ and SW are the soil water content at the beginning and end of a day, Δt is the time-step (24 h), and TT is time (h) through layer I (Arnold et al. 1998). The percolation rate (mm d⁻¹) is calculated as ⁠. In the soil profile from 0.0 to 2 m, lateral subsurface flow is estimated along with percolation. A kinematic storage model calculates the lateral flow in each soil layer. Return flow from the shallow aquifer to the stream is calculated as (Arnold et al. 2000). Where q_lat is the lateral flow (mm d^−l), S is the volume of soil water (m h⁻¹), α is the slope (mm⁻¹), θ_d is the drainage porosity (mm⁻¹), and L is the flow length (m) (Arnold et al. 1998).

This study uses Penman–Monteith method (Monteith 1965) for estimating potential evapotranspiration (PET). Exponential functions of soil depth and water content are used for predicting soil evaporation, and a linear function of PET and leaf area index is used for simulating plant water evaporation (Arnold et al. 1998). When the soil second-layer temperature is simulated more than 0.0°C, snow may be melted, which is simulated as where SML is the rate of snowmelt (mm d⁻¹), SPT is the snowpack temperature (°C), and T is the average daily temperature (°C). In many semiarid regions, river flow and streamflow are hydraulically linked with alluvial floodplain groundwater. Transmission is simulated using Lane's method (USDA 1983) to simulate losses as flood wave travels downstream (Arnold et al. 1998).

Data used in modeling include climate (rainfall, air temperature), hydrologic processes (flow rate), and spatial (topography, land use, soil). The type and sources of the data used in this study are provided in Table 1.

The ecohydrological model for the entire TR watershed and its subwatersheds was developed in this study based on the available data. The SWAT is used to simulate hydrological processes in the watershed.

As it has been widely used in scarce data watersheds to simulate various processes in watersheds, the SWAT model is selecting to simulate and understand the ecohydrological processes in the TR watershed. The SWAT model was previously used for two watersheds in Mongolia. The effect of subarctic conditions on the Kharaa river basin, located in the northern part of Mongolia, was simulated using the SWAT model. It was found that SWAT satisfactorily reflects streamflow for particular years but is not reliable for a more extended period (Hülsmann et al. 2015). Sukhbaatar et al. (2017) developed the SWAT model for the Upper TR watershed to analyze the impact of climate change on hydrological processes. They found that the model sufficiently predicted the hydrological processes in the watershed and further concluded that the river flow was primarily influenced by precipitation.

The TR watershed was analyzed based on the HRU that represents homogeneity in land and soil characteristics. Water resources for the TR watershed were simulated from 1996 to 2011 using observed precipitation, temperature, wind speed, and air humidity inputs. The streamflow data from two stations (Ulaanbaatar and Terelj) in the watershed are used to calibrate and validate SWAT model results.

The SWAT model was developed using the topography, land use, soil, and meteorological data from the TR watershed. SWAT parameters were derived from spatial analysis of the watershed using GIS and procedures suggested in the SWAT manual. Spatial inputs to the SWAT model (Figure 5) include (a) topography, (b) land use/cover, (c) soil class, (d) surface slope, and (e) subbasins. The global digital elevation model, generated from satellite imagery using an Advanced Spaceborne Thermal Emission and Reflection Radiometer (NASA JPL 2009), was used. The topography has a resolution of 30 × 30 m. Land-use/cover data of the European Space Agency (ESA) was used in modeling. The GlobCover initiative of ESA uses the Envisat MERIS data of Fine Resolution (300 m) (Arino et al. 2012; Sukhbaatar et al. 2017). The soil map is from the Harmonized World Soil Database from FAO (Food and Agriculture Organization). The TR watershed is divided into 27 subwatersheds based on the unique HRUs (Figure 5).

Calibration and validation of the SWAT model performed using the coefficient of determination (R²), RMSE of observations standard deviation ratio (RSR) (Moriasi et al. 2007), and Nash–Sutcliffe model efficiency (NSE) (Nash & Sutcliffe 1970). The R² is the degree of agreement between simulated and measured data, measured as ⁠, where Q is the dependent variable, and m and s stand for observed and simulated for data pair. Higher values of R² indicate higher explanatory power. NSE determines the relation between residual variance compared with measured variance. ⁠, where the bar represents an average value. NSE ranges between −∞ and 1, with an optimal value of NSE = 1.

The varies from 0 to large positive values, with lower values of RSR indicating a better model fit (Moriasi et al. 2007).

Santhi et al. (2001) suggest that the SWAT calibration results are acceptable of R² and NSE. In addition to that, Moriasi et al. (2007) also summarized the overall performance rating for recommended statistics for monthly time-step (Table 2).

The SWAT model calibration was conducted using the SWAT-CUP (SWAT-Calibration Uncertainty Procedures) program. SWAT-CUP is designed for calibration and validation, sensitivity, and uncertainty analyses of SWAT results (Arnold et al. 2012; Abbaspour et al. 2015; Boithias et al. 2017). SWAT-CUP analysis and Sequential Uncertainty Fitting version 2 (SUFI2) algorithm (Abbaspour et al. 2007) were used for sensitivity analysis and calibration. The sensitive analysis was conducted to assess the influence of parameters on river flow prediction in the SUFI2 algorithm. Other methods for sensitivity analysis used in hydrologic studies include the Generalized Likelihood Uncertainty Estimation (GLUE) algorithm (Beven & Binley (1992) and statistical factorial methods (Wang et al. 2020). SUFI2 was found superior to GLUE in runoff simulations in forest watersheds (Nohegar et al. 2017). The obtained model is calibrated and validated using daily streamflow data from two hydrological gauges (Terelj and Ulaanbaatar stations) in the headwaters. The models' calibration and validation were conducted using the same period with two gauging stations available in the watershed. The model was calibrated with 1990–2011 data using the discharge data from the Ulaanbaatar station using statistical methods listed in Table 2. The model validation was conducted using the Terelj station data ranging from 1990 to 2011 as internal validation.

The SWAT model simulated the hydrological processes in each subwatershed for 25 years (1987–2011). The results are compiled for daily, monthly, and yearly values of flows. Daily values had low validation and calibration values. This was due to missing and limited observed values. Therefore, the initial three years are used for the warm-up. From model results, only annual and monthly average basin values were used to study seasonality and yearly trends in water balance components for every 27 subwatersheds.

The calibration results for the monthly time-step using the simulated and observed discharge data values at Terelj and Ulaanbaatar stream (Table 3). The statistical values for calibration and validation were good (R², NSE, and RSR) and satisfactory (PBIAS) (Moriasi et al. 2007).

The parameter global sensitivity analysis is conducted in the SUFI2 by regressing a Latin hypercube generated parameters against the objective function (Khalid et al. 2016). The most sensitive parameters are SOL_AWC(1).sol, GW_DELAY.gw, SOL_K(1).sol, ESCO.hru, CH_N2.rte, and CH_K2.rte. (Figure 6). The available watershed capacity (SOL_AWC) is a major factor in semiarid regions that influence the vegetative cover, water uptake, and runoff generation. The time of vadose zone flow past the lowest depth (GW_DELAY) is sensitive as this subsurface flow rate is influential in semiarid, cold climate watersheds. The saturated hydraulic conductivity (SOL_K) in layer 1 influences the rate of flow of water in the soil, which is critical in influencing runoff, subsurface flows, infiltration, and stream discharge. The soil evaporation compensation factor (ESCO) defines the soil evaporative demand and influences soil water through subsurface loss in water that influences the balance in subsurface flows. The channel parameters in the TR (CH_N2 and CH_K2) are sensitive given the variability in geomorphology and in-channel processes (Pietroń et al. 2015).

The observed and simulated monthly streamflow values at Terelj and Ulaanbaatar stations are shown in Figure 7. The average annual water balance components in the TR Basin are: precipitation = 277.5 mm, surface runoff Q = 0.44 mm, lateral soil flow = 58.14 mm, groundwater (shallow aquifer) Q = 4.50 mm, groundwater (deep aquifer) Q = 0.24 mm, revap (shallow aquifer ≥ soil/plants) = 0.66 mm, deep aquifer recharge = 0.24 mm, total aquifer recharge = 4.73 mm, total water yield = 63.31 mm, percolation out of soil = 4.73 mm, ET = 214.3 mm, PET = 696.0 mm, and total sediment loading = 0.16 t/ha. The model simulated precipitation = 277.5 mm and snowfall = 46.26 mm, which fall within the range reported by Dolgorsuren et al. (2012) and Sukhbaatar et al. (2017). Moreover, the TR watershed average has 45.15 annual basin water stress days and 133.79 temperature stress days.

The average monthly basin values are calculated and shown in Figure 8. The maximum amount of rain is simulated in July at 76.16 mm. The maximum amount of surface runoff was 0.12 mm in August when the lateral flow (16.43 mm) and water yield (17.74 mm) were also the highest. The lowest value of precipitation was observed in January. The average annual water yield was 63.31 mm. The seasonal flows show a typical water balance among the months.

Some variations are observed between the actual and potential ET measurements. This is due to input variables in Hamon PET that are based on day length and air temperature only. In contrast, actual ET is influenced by additional factors such as solar radiation, relative humidity, and wind speed. The changes in ET are directly correlated to the growing season and type of vegetation in the watershed. The ET gradually decreases from September to January, as deciduous trees lose their leaves, reduced demand for ET in September–October, and drop in air temperature. As the dormant season starts, the ET decreases and reaches a minimum in January. This is because of snow cover accumulating on evergreen trees and the ground. As days get warmer in spring and day length increases, ET also increases, starting in February. When the growing season begins in May, demand for ET also increases rapidly. The highest potential ET was observed in July, after which the ET again decreased as air temperatures subsided as days get shorter.

The time series of average monthly discharge and monthly mean air temperature, and average monthly precipitation in the TR watershed are shown in Figure 9. According to the simulation results of the average monthly discharge values in the TR watershed, 1993–1995 was a wet year. Starting from 1996 onwards, the simulated discharge values are lower than the previous years due to the warmer temperature and lower precipitation.

The time series of monthly average discharge values at each subbasin in the TR watershed is provided in Figure 10. There is a decreasing trend of monthly discharge over the years since 1996. The seasonal average discharge values at each subbasin in the TR watershed are provided in Figure 11.

The annual average discharge at each subbasin is spatially represented in Figure 12. The discharge values increase along the river and reach their highest value in subbasin one, the whole basin outlet.

Water scarcity is a critical issue in semiarid, cold climate regions with increasing land use. Simulation of the spatial and temporal processes can help understand the factors influencing scarcity and inform decisions at a watershed scale. This study contributes to the literature in simulating the dynamics to evaluate spatial and temporal patterns in the semiarid TR watershed less studied. Furthermore, the results will help manage water scarcity at a watershed scale for changing land-use conditions.

The SWAT model has performed well in simulating hydrologic processes in large watersheds worldwide. The model was used to simulate the watershed processes in the TR watershed with unique climatic and geomorphic conditions. Split-location, calibration, and validation are performed for monthly time-steps using observation at Terelj and Ulaanbaatar stations. Calibration and validation were conducted with SUFI2 within SWAT-CUP methods using both graphical and statistical techniques. The R², PBIAS, RSR, and NSE results demonstrate that the SWAT model effectively replicates annual, monthly streamflow values. In addition, we observe that the simulation model reasonably captures the watershed processes in the TR watershed.

Spatial and temporal variations in the streamflow discharge for the TR watershed were observed in the model simulation. The annual average discharge ranged from 0.3 to 96.87 m³/s in the basin, indicating a sizeable inter-annual variation in water flows. With water supplies facing high uncertainty, there is a need for watershed-wide management of water supplies and demand. The river flow was primarily influenced by climatic influences through precipitation and air temperature changes. This emphasizes the role of climate change in affecting the watershed hydrologic processes. With high actual ET in the watershed, decreasing river flows. We observe that spatial and temporal variations in hydrologic components are significant in the TR watershed.

The variability over geographic space and time is a fundamental issue in water management for scarcity. The methods and results of this study will be helpful for replication in similar semiarid watersheds elsewhere for addressing water scarcity at a watershed scale. In addition, information on variations over space can be used to target specific subwatersheds to address water scarcity issues.

The study has some limitations as it depends on the monitoring information for the calibration and validation of the simulation model. Data availability is often a major limitation in watersheds like TR, and future extensions can focus on developing and incorporating new datasets in the region. Other study extensions could quantify uncertainty in the area resulting from land use and climate change on watershed processes.

We thank the US Fulbright Program for providing the scholarship to pursue this research. Partial support is provided by NIFA-USDA, Massachusetts Agricultural Experiment Station (MAES) under projects MAS00036, MAS00035. We thank Dr David Bloniarz for his advice in the research.

All relevant data are available from publicly available datasets listed in methods (Table 1).