Received 21 March 2016; accepted 9 July 2016; published 12 July 2016
The Republic of Uzbekistan is a doubly land locked country in central Asia. Therefore, most water resources are supplied from neighboring countries. In the last decades, water supply into to the country has significantly decreased due to the mismanagement of transboundary water resources by Central Asian countries (  -  ). Consequently, most river basins in the region suffer from water shortages and increased aridity. Global climate change threatens to worsen these life-threatening problems  . The development of plans for the sustainable usage of water resources, and mitigation and adaptation strategies is a key concern for all stakeholders across Uzbekistan  . The Chirchik River Basin (CRB) is considered one of the country’s largest and most important basins. It is located in the northeastern province of Tashkent (Figure 1), and its economy relies heavily on agricultural production; employing most of the region’s population and consuming the lion share of available water through irrigation practices  . Unsurprising, the CRB is facing several potentially severe water related problems, including increased aridity and, land salinization, and has been declining in agricultural production. Continuously, rising groundwater level and the inefficient use of river waters are assumed to further cause these problems  .
Poor management of water resources due to limited knowledge of hydrological processes has been identified as a root cause of water resource problems in the basin (   ). Efficient management requires accurate estimation and modeling of water balance and main hydrological parameters. In order to obtain useful results, irrigation processes, which play an important role in the hydrological cycle, should be included in the model to account for the contribution of irrigation water to other hydrological parameters. MIKE SHE coupled fully distributed hydrological models have been used widely and effectively by many researchers in producing detailed water balance estimations, examining hydrological responses to land use and cover change, and in groundwater and irrigation management (  -  ). Compared with other conventional methods, the model used in this study has the advantage of fully integrating the surface, subsurface, channel flow and their interactions in hydrological process simulation.
Moreover, this model also has an integrated system to calculate irrigation processes  . Further evaluation and comparison of the MIKE SHE and other hydrological models are provided by (  -  ). Integrated models have rarely been used in CRB studies, and there is a paucity of literature estimating water balance and hydrological parameters. In this concern, a main aim of this work is to study hydrological processes in the basin territory through detailed water balance estimation using an integrated hydrological model. This research aims to contribute to the knowledge of groundwater and surface water interaction and spatial variability of hydrological parameters.
2. Method and Materials
2.1. Study Area
Figure 1. Tashkent province and the Chirchik river basin, Uzbekistan.
cotton and wheat. The largest expanse of irrigated land is located in the middle and downstream areas of the basin  . Chirchik River water is distributed to agricultural lands and urban areas via the Bozsuv, the Qorasuv and Parkent’s main canals. The remaining volume is mixed with returned water from agricultural lands and urban areas, and discharged into the Syrdarya River.
2.2. Model Construction
The MIKE SHE model is a deterministic, fully distributed and physically based hydrological modeling tool. The model calculates the following major hydrologic process in the hydrological cycle by fully integrated basis: evapotranspiration, overland flow, unsaturated flow, saturated flow, and stream flow. MIKE SHE uses the MIKE 11 model to simulate stream flow, lakes, and channels in 1D. The model uses Kristen and Jensen method to calculate actual evapotranspiration (AET). This method requires leaf area index (LAI), and root depth (RD) parameters for each vegetation type. LAI is dimensionless quantity and defined the area of leaves per unit ground surface area (LAI = leaf area/ground area). Typical value of the LAI varies between 0 and 7. The root depth is defined depth of roots in the roots zone of each vegetation type. This is used to estimate the amount of water extracted by the roots between the ground surface and the lover level of saturated zone. In MIKE SHE, the temporal variation of the LAI and RD must be specified for each vegetation types.
The Richard’s 1D equation is used to calculate unsaturated flow process. The saturated flow component is calculated using the Darcy’s 3D method. Also, snow melting and freezing processes are simulated using a degree-day empirical approach, which requires a degree-day coefficient of the study area. The MIKE 11 model uses a 1D Saint-Venant equation to compute stream flow computations. Finite element methods are used to solve these partial differential equations. MIKE SHE uses a network of regular grids to discretize the horizontal plane of a watershed, and represent the spatial variability of the calculated hydrological process. The model also simulates irrigation processes, pumping wells and various water control structures. Simulation of the irrigation processes requires irrigation command areas, water sources (rivers, wells or lakes), irrigation time and demand. The model simulates water releases from stream flow with these releases then converted to irrigation depth to over-irrigation area. A detailed explanation of MIKE SHE and MIKE11 hydrological models can be found in  and  .
2.3. Model Set Up
The MIKE SHE hydrological model requires a spatially distributed data set, which includes topographical, soil, geological, land use, climate and initial potentiometric head data. These geospatial data were prepared with the ArcGIS 10.2 platform and converted to 500 m by 500 m grid format. All these geospatial data were projected to CRB’s coordinate system (WGS84 UTM 42N). Surface topographical information was obtained from ASTER (Advanced Space borne Thermal Emission and Reflection Radiometer) GDEM (Global Digital Elevation Model). All data output time steps were set to 24 hours.
2.4. Unsaturated Zone
Unsaturated flow computation in the model requires a spatially distributed soil map, and its hydraulic parameters. The distributed soil map was obtained by digitizing 1:800,000 scale soil map of the Tashkent province. Soil types vary depending on geomorphological zoning in the basin. As total 12 soil types were identified in the basin region (Figure 2). This is a generalized soil map as it groups closely similar soil types as one. In this study Richards’s equation was selected for calculating vertical flow in unsaturated zone, which are the functions of soil moisture retention curve and effective conductivity.
2.5. Saturated Zone
A 3D finite difference method was selected to compute saturated zone flow. This method requires the spatially distributed thickness of the computational geological layers. Only the surficial soil layer was considered in this research due to insufficient information of geological layers within the basin. Soil layer thickness was constructed using lithological information from a geological borehole dataset. The depth of the soil layer continuously increases from the mountainous to the downstream areas of the basin, ranging between 5 - 35 meters. The assumption of soil-layer thickness is consistent within the geological literature  . Initial groundwater table data was generated through a geostatistical interpolation method using data from 90 groundwater wells. The specified groundwater table is used as the lower boundary condition for the unsaturated model.
2.6. Land Use of the Basin
A land use map of the CRB was obtained from AVNIR-2 (Advanced Visible Near Infrared Radiometer) satellite image through a supervised classification method using the Image Analysis toolset in ArcGIS 10.2. The primary land use data were calibrated with minor modifications performed using paper-based cadastral land use maps of districts in the Tashkent province. The land use of the basin has been classified as urban, water body, farmland, forest, grassland and arable lands (Figure 3). The values of seasonal changes of LAI and RD data were obtained from the Institute of Water Problems of Uzbekistan, and assigned to each land use class.
Figure 2. Soil map of the Chirchik river basin.
Figure 3. Land use map of the Chirchik river basin.
2.7. Meteorological Input
As shown in Figure 1, the CRB contains only eight weather stations. Climate data from each station was provided by the Uzhydrometeorological Authority of Uzbekistan. Potential evapotranspiration (PET) was calculated using the Penman-Monteih FAO 56 model  . The Thiessen polygon method was used to spatially distribute daily-accumulated precipitation, PET, the degree-day melting coefficient and short wave solar radiation data. Proper lapse rates were set to precipitation and temperature data for correction of the Thiessen polygon method in mountainous areas.
2.8. MIKE 11 Model
MIKE SHE uses the MIKE 11 model to simulate channel flow components. The MIKE 11 requires river network, cross section and discharge data of simulated streams. River and channel network data were extracted via digitization of the topographical map of the Tashkent province. Cross-section data were collected by field surveys in each of the simulated branches. Daily average discharge data of the Ugam River and the Charvak water reservoir were assigned to the upstream boundaries of the network. The created Chirchik stream network consists of four branches and 28 cross-section data. Storing frequency of simulation time step was set to 6 sec 6 min. The constant Manning’s roughness coefficient was applied for each branch in the stream network. Time series data of water withdrawal rate of simulated rivers were set after coupling MIKE 11 flow simulation to the MIKE SHE model. After coupling MIKE 11, flow simulation to the MIKE SHE model irrigation area and demand were set to simulate irrigation processes. The model was set to obtain 100% of irrigation water from the Chirchik River.
3. Results and Discussion
Simulation results were calibrated and validated for 2009-2011 and 2012-2013 against stream flow and groundwater level data. Model performance was evaluated using mean error (ME), mean absolute error (MEA), root mean square error (RMSE), correlation coefficient (R2) and the Nash and Sutcliffe coefficient (EF)  . Initially, climate parameters for modelling snow melting and freezing, including degree-day coefficient and maximum wet snow fraction, were adjusted to climatic conditions of the basin. The next step of calibration process focused on adjusting horizontal and vertical conductivity, and specific storage of saturated zone parameters.
3.1. Stream Flow Hydrograph
Table 1 presents performance statistics of streamflow simulation. Figure 4 and Figure 5 presents the comparison of simulated and observed discharges for the calibration and validation period in the Chinoz gauging station. The model was generally good, simulating daily streamflow discharge with an average of 0.91 (2009-2013) and 0.77 (2009-2013) of R2 and EF. The hydrographs show over and under estimation of streamflow discharge. Simulated river flow shows numerical imbalance when river water starts to be withdrawn for irrigation. This is clearly shown in March 2013 (validation period). Relatively high overestimation was detected from July to
Table 1. The statistics of streamflow simulation at Chinoz gauging station.
Figure 4. Observed and simulated streamflow in Chinoz gauging station (2009-2011).
Figure 5. Observed and simulated streamflow in Chinoz gauging station (2012-2013).
November immediately after peak discharge rate. The main reason for this discrepancy is the absence of the operation information of small irrigation canals, reservoirs and hydro-engineering structures. Not accounting for the exchange between river water and saturated zone due to insufficient geologic data may also explain this.
3.2. Groundwater Dynamics
Simulations of groundwater level dynamics were calibrated and validated against average monthly observed well data for 2009-2011 and 2012-2013 respectively. In this study, data from a total of eight ground water wells was used.
These groundwater wells are located at up (W5, W6), middle (W4, W3) and downstream (W1, W2, W7 and W8) sites in the basin (Figure 1). The location of groundwater wells assists appraisal of model performance in different zones. The performance of groundwater simulations are provided in Table 2 and Table 3. Figure 6 and Figure 7 show comparison of observed and simulated groundwater depth. Simulated groundwater dynamics show a satisfactory match with observed data at each location. Hydrographs show the groundwater table reaching its maximum level from April to May, and dropping to minimum levels at the end of September and October. In general, the model underestimated the groundwater level during peak recharge period. All hydrographs show consistent dependence on precipitation in terms of timing and quantity. Simulated groundwater elevations across the middle and downstream sites varied up to an average of 0.66 m during the study period in response to precipitation events. This value was equal to an average of 1.7 m in upstream sites. This shows that the fluctuation amplitude of the groundwater table in upstream sites is relatively higher. Potentiometric surface maps show the
Table 2. Calibration results of groundwater simulation.
Figure 6. Observed and simulated groundwater depths in calibration period.
general direction of groundwater flows from Northeast to Southwest towards the Syrdarya River. This clearly underlines how groundwater flows contribute hugely to the Syrdarya River through base flows. This is an estimated average of 55 mm/year.
Table 3. Validation results of groundwater simulation.
Figure 7. Observed and simulated groundwater depths in validation period.
3.3. Water Balance
The hydrological balance of the CRB is strongly dependent on agricultural activity and the climatic conditions of the basin. The model simulated the main hydrological processes, including evapotranspiration, overland unsaturated and saturated zone storage changes for calibration (2009-2011), and validation (2012-2013) periods. Table 4 shows the contribution of each water balance component in the CRB (mm/year). In simulation, precipitation and irrigation water are the main hydrological inputs to the water cycle. The overland water from snow and glacier melting is estimated at around 4% of the total inflow. The results indicate the main proportion of this water was lost through AET during the vegetation period (March-September) due to applied irrigation water.
The total AET was estimated at an average of 821 mm/year. This corresponds to 77% of the total water budget. The estimated average AET is exceeds around 8% from the average annual precipitation. The AET from snow surface was an estimated average of 89 mm/year or 10% of the total average AET. AET is highly variable across the watershed, and it predominantly increases in downstream sites in the basin. This occurs because PET increases with decreasing elevation  . The spatial range of AET is a varied average of 1211 - 692 mm/year. High AET (1211 mm/year) was estimated in downstream irrigation areas. This is almost two times greater than precipitation, and it occurs during non-rainy periods. The high evaporative power of the atmosphere is to be compensated by irrigation water.
It should be stated that decreased irrigation water causes of increased aridity in the basin. Simulation results show that the main recharge of groundwater is mainly provided by precipitation and irrigation water. During the study period, the total recharges varied from 180 to 221 mm/year (the last row in Table 4). This ranges from 17 to 20% of the total inflow. Groundwater recharge largely occurs from March to May (average 90 - 60 mm/ month) when the amount of rainfall and snow melting increases in middle and downstream sites. In contrast, mountainous area recharge lasts until the end of June, and then starts to decrease, depending on variations in the amount of rainfall. A maximum recharge of 221 mm/year was estimated in 2010, when high amounts of precipitation occurred. The results indicate that in mountainous areas around 51% of total precipitation falls as a snow during the cold season. This proportion equaled an average of 29% in downstream sites. Snow storage mainly occurs from December to the end of February. Additionally, the results indicate overland flow peaks in mountainous areas of the basin in spring and fall with seasonal heavy rains and the beginning of snow melt on low infiltration ground surfaces. Snow pack melt-water contributed an average of 34% of the total overland flows. In this way, overland flow substantially feeds the source of the Chirchik River. This amount is an estimated average of 77 mm/year, which is 3% of the average total of precipitation. A small amount of water outflows from the basin boundary, particularly from downstream sites. This is an estimated average of 89 mm/year, or 4% of total precipitation. The overall water balance error was averaged at 2% and 3% during the calibration and validation periods. This shows incoming and outgoing water balance parameters were well satisfied.
The MIKE SHE integrated hydrological model was used to study the hydrological balance of the CRB. The
Table 4. Contribution of hydrological parameters to water balance.
model was constructed with several assumptions when preparing the geospatial data set of unsaturated and saturated zones, and climate and land use parameters. This hydrological study of the CRB allows for the following main conclusions:
・ MIKE SHE as a physical-based mathematical model uses a large number of hydrological parameters. All need to be adjusted to the climatic conditions of the study area to reach satisfactory results that correspond to reality.
・ Evapotranspiration was found to be the main water loss factor among water balance components, with an average of 821 mm/year (77% of the total water budget). As an arid land, AET is strongly dependent on irrigation water quantity irrespective of rainfall in downstream sites.
・ AET is highly variable across the basin, with increases toward to downstream sites. The estimated AET in irrigation areas at downstream sites was almost two times higher than upstream sites. Therefore, decreased irrigation water causes increased aridity, particularly in downstream site.
・ Estimated groundwater recharge varied between 180 - 221 mm/year, making up 17% - 20% of the total water budget. The highest groundwater recharge occurs from March to May with an average of 90 - 60 mm/ month.
・ The general direction of groundwater flow is toward the Syrdarya River. Base flow from basin boundary into the Syrdarya River was estimated at an average of 55 mm/year.
・ The Chirchik River is gaining upper stream sites by overland flow on average 77 mm/year.
・ The main water balance error was obtained from simulated overland flow. Adjusting the climate parameters of the model to the basin environment optimizes this. The overall water balance error is estimated 2.5% on average, which demonstrates that interacting hydrological components and parameters are well matched.
・ Accurate irrigation schedule and operational information of hydro-engineering structures are critical model inputs in order to improve daily streamflow simulation and minimize numerical errors.
・ This study confirms that in arid and semi-arid regions with intensive agricultural, integrated modelling is valuable for understanding water cycles in large basins.
The novelty of this study is a development of a local hydrological model for a relatively large river basin covering from upstream to downstream site of the basin and considering agricultural water use and administrative boundaries of the districts. This modeling approach is better describing formation, utilization and discharging of water resources under the impact of different land use process. This novel tool is urgent for integrated water resources management, particularly assessment of quantitative status of surface and groundwater resources. The comprehensive part of the research was the processing all required geospatial data (soil type, land use, topography, geology, river network geometry and agricultural water use) into integrated numerical model. This is challenging because all these data have never been used together before in hydrological studies in Chirchik River Basin. Using the model and model results, further steps of study should focus on predicting effects of climate and land use change to the hydrological regime of the CHRB. A part from quantitative analyses, in future research should also consider surface and groundwater quality issues.