MODELLING OF HYDROLOGICAL PROCESSES IN THE CATCHMENT AREA OF LAKE SKADAR

Original scientific paper Paper deals with the modelling of hydrological processes in the catchment area of Lake Skadar in Montenegro. This lake and its catchment area is one of the most important natural resources of Montenegro. We used an existing model (Mike SHE (DHI)) to simulate the hydrological processes. Its implementation has shown some limitations that are primarily related to insufficient funds of available input data. A particular problem is the lack of data on groundwater regime as one of the key model and water balance factors. However, calibration of the model showed that the results obtained were sufficiently acceptable and that the model can be a good basis for future hydrological forecasts of the basin size. Also, the model can still be upgraded with the new data and results, as soon as the conditions call for it.


Introduction
Hydrological models are important for a wide range of applications, including water resources planning, development and management, flood prediction and design, and coupled systems modelling including, for example, water quality, hydro-ecology and climate. However, due to resource constraints and the limited range of available measurement techniques, there are limitations to the availability of spatial-temporal data; hence a need exists to extrapolate information from the available measurements in space and time; in addition there is a need to assess the likely hydrological impact of future system response, for example to climate and land management change. In this paper hydrological modelling is applied to verify the hydrological process understanding developed during the study period.
Several problems arise when working with hydrological models, such as data availability (especially when working with physically based models), spatial diversity between observed parameters and model parameters, differences between hydrological process scales and modelling scales.
The choice of a model is determined by the purpose of the model and the availability of data [1]. There are numerous watershed-scale hydrologic models. The choice of models should be based on the objectives of use. Literature reviews suggest that the MIKE SHE⁄MIKE 11 modelling package (DHI, Hørsholm, Denmark) [1] has several advantages over other hydrologic models for estimating watershed runoff: (1) it is a distributed model and most of the algorithms in describing the water movements are based on physical processes, (2) it simulates the overland flow processes commonly found in dry regions, and (3) it has been commercialized and a GIS user interface was built in the system that can directly use spatial GIS databases for model inputs. Also, the model has a strong visualization utility that makes interpretation of modelling outputs much easier.
MIKE SHE covers the major processes in the hydrologic cycle and includes process models for evapotranspiration, overland flow, unsaturated flow, groundwater flow, and channel flow and their interactions. [2,3] Each of these processes can be represented at different levels of spatial distribution and complexity, according to the goals of the modelling study, the availability of field data and the modeller's choices. The MIKE SHE user interface allows the user to intuitively build the model description based on the user's conceptual model of the watershed. The model data is specified in a variety of formats independent of the model domain and grid, including native GIS formats. At run time, the spatial data is mapped onto the numerical grid, which makes it easy to change the spatial discretization.
MIKE SHE, in its original formulation, could be characterized as a deterministic, physics-based, distributed model code [2]. It was developed as a fully integrated alternative to the more traditional lumped, conceptual rainfall-runoff models. A physics-based code is one that solves the partial differential equations describing mass flow and momentum transfer. The parameters in these equations can be obtained from measurements and used in the model. There are, however, important limitations to the applicability of such physics based models. For example: • it is widely recognized that such models require a significant amount of data and the cost of data acquisition may be high; • the relative complexity of the physics-based solution requires substantial execution time; • the relative complexity may lead to overparameterized descriptions for simple applications; and • a physics-based model attempts to represent flow processes at the grid scale with mathematical descriptions that, at best, are valid for small-scale experimental conditions. Catchments typically have a high level of spatial heterogeneity which can be prohibitively expensive to observe or comprehensively represent in the model. This is most obvious in the representation of subsurface processes because of the difficulty of observation and the high degree of soil/aquifer heterogeneity which often exists. In principle the parameters of physics-based models are measurable, but in practice this cannot be achieved at the scale of modelling application, because such measurements are essentially made at a point. Therefore, these models use averaged variables and parameters at grid or element scales which are greater than the scale of variation of the processes. Even under a "full" physics representation, parameterization (including spatial variability of parameters) of material properties does not represent catchment heterogeneity. [4] 2 Catchment properties and hydrometeorological data of the Skadar lake basin

Catchment characteristics
The Morača river basin is located in the central part of Montenegro. The river originates in northern Montenegro, under Rzac mountain at an altitude of 975 masl. In the upper and middle part of its flow, the Morača river is a highly mountainous river that has cut a canyon with steep slopes. Its length is 113.4 km, and the area of the river basin to the hydrological station Podgorica is 2628 km 2 [5].
The Zeta is a right tributary of the river Morača. It appears as a stream Surduk, which in Nikšićko polje, after having received two small tributaries, becomes the Sušica river. Downstream of the river Rastovac confluence, the river changes its name to Zeta. During the low water period the river delves and continues to flow underground in a straight line, for about 5 km, and reemerges from the cave on the southern slope of the Planinica, at an altitude of 100 masl. At the village Tvorila the river creates the Slap waterfall. Through Bjelopavlići plain the Zeta meanders quietly. At the debris of the ancient city of Duklja, the river flows into the Morača river as its main tributary. Its length is 85 km, and the river basin area to the most downstream hydrological station Danilovgrad is 1216 km 2 . After the merge with its largest tributary, the river Zeta, the Morača enters the Zeta plain and flows to its confluence to the Lake Skadar. The Morača river is the largest tributary of the Skadar Lake (Fig. 1).
The steepest hillside slopes are located in the upstream and middle part of the Morača river flow. The inclinations are up to 69 degrees. The slopes of the Zeta river basin are significantly milder and they vary between 0 and 30 degrees. Some smaller areas of the basin have also high terrain slopes exceeding 45 degrees.

Figure 1 Flow Accumulation derived from DEM with Arc Hydro Tools
The spatial distribution of the height on the Morača catchment is given by the Digital. The watershed elevation ranges from 22 to 2217 masl. The mean catchment elevation is 981.53 masl.

Climate data
The precipitation data on the Morača river basin are  Precipitation stations are evenly distributed with altitude on the catchment area, but the temperature data are available only for the stations under 1000 masl. Since the catchment has a complex relief, it is necessary to have more precipitation and temperature data at all altitudes for the adequate simulation of the hydrological processes.
The precipitation distribution by months and seasons indicates that here a long rainy period with a humid climate and a shorter summer period with an extremely arid climate can be separated (Fig. 3). In the period from October till December the mean monthly precipitation is 873.08 mm, which is 40.16 % of the annual total and from January to April the average precipitation is 796.43 mm or 36.56 % of the total annual precipitation. The largest amounts of precipitation are measured at Andrijevo station, followed by Danilovgrad and Kolašin stations.

Hydrological data
The Moraca river catchment area for the purpose of this study is determined by hydrological profile of Podgorica. Monthly mean discharge values available for this study are those recorded from January 1948 till December 2012. Mean multiannual discharges for Podgorica profile are 158.05 m 3 /s and 189.66 mm respectively.
The seasonal cycle of monthly mean discharge at Podgorica station exhibits a high discharge of 184 m 3 /sto 284 m 3 /s from November till May and a relatively low discharge of 26 m 3 /s to 56 m 3 /s from June to October with maximum discharge usually in December and in August. The interannual variation of monthly mean discharge is generally small in the dry season and large in the wet season ( Fig. 4). MIKE SHE model integrates the entire land phase of the hydrologic cycle and can model interception, actual evapotranspiration, overland flow, channel flow, flow in the unsaturated zone, flow in the saturated zone and exchange between aquifers and rivers. MIKE SHE applied at a catchment scale implies the assumption that smaller scale equations are valid also at the larger scale; thus, it performs an upscaling operation using effective parameter values.
MIKE SHE consists of the following model components: • Precipitation (rain or snow), • Evapotranspiration, including canopy interception, which is calculated according to the principles of Kristensen and Jensen, • Overland flow, which is calculated with a 2D finite difference diffusive wave approximation of the Saint-Venant equations, using the same 2D mesh as the groundwater component. Overland flow interacts with rivers, the unsaturated zone, and the saturated (groundwater) zone, • Channel flow, which is described through the river modelling component, MIKE 11, which is a modelling system for river hydraulics. MIKE 11 is a dynamic, 1D modelling tool for the design, management and operation of river and channel systems. MIKE 11 supports any level of complexity and offers simulation tools that cover the entire range from simple Muskingum routing to high-order dynamic wave formulations of the Saint-Venant equations [6], • Unsaturated water flow, which in MIKE SHE is described as a vertical soil profile model that interacts with both the overland flow (through ponding) and the groundwater model (the groundwater table is the lower boundary condition for the unsaturated zone). MIKE SHE offers three different modelling approaches, including a simple 2-layer root-zone mass balance approach, a gravity flow model, and a full Richards's equation model, • Saturated (groundwater) flow, which allows for 3D flow in a heterogeneous aquifer, with conditions shifting between unconfined and confined. The spatial and temporal variations of the dependent variable (the hydraulic head) are described mathematically by the 3D Darcy equation and solved numerically by an iterative implicit finite difference technique.
For integrated process-oriented catchment modelling MIKE SHE has different numerical methods from simple conceptual approaches to advanced, physics-based methods. That possibility enables modelling in situations when some model parameters are not available, therefore some processes could be conceptualized. All modules for hydrologic processes can be linked or unlinked with each other easily. For the best simulation of the hydrological processes of the Morača river basin all available spatial and temporal data are included in the MIKE SHE Flow model.
The model for the Morača river basin was built with modules for Overland Flow, Unsaturated Flow and Evapotranspiration. Modelling of the Saturated Flow needed information and data that were not at the disposal at the moment. Simulated data were then used as input to the Water Balance Module, where the components of the water balance were calculated.

Model setup
In most watershed problems one or two hydrologic processes dominate the watershed behaviour. Thus, a complete, physics-based flow description for all processes in one model is rarely necessary. A sensible way forward is to use physics-based flow descriptions for only the processes that are important, and simpler, faster, less data demanding methods for the less important processes. The downside is that the parameters in the simpler methods are usually no longer physically meaningful, but must be calibrated -based on experience.
Water flow on the ground surface is calculated using a finite-difference, diffusive wave approximation of the Saint -Venant equations. Due to spring snow melt discharge it was necessary to include the Snow Melt module and have adequate annual discharge distribution. The Snow Melt module uses a Degree-day melting algorithm. The melting of snow and ice is assumed to be related to air temperature as long as air temperature is above a critical threshold. The amount of snow or ice melted at a certain place, during a certain period, is assumed proportional to the sum of positive temperatures (on the Celsius scale) at the same place and in the same period. The amount of melt is linked to this positive degree-day sum by the degree-day factor.
The main purpose of the Two-Layer Water Balance ET method is to provide an estimate of the Actual evapotranspiration and the amount of water that recharges the saturated zone. The model divides the unsaturated zone into a root zone, from which evapotranspiration can be extracted, and a zone below the root zone, where evapotranspiration does not occur. Evapotranspiration is extracted first from intercepted water (based on the leaf area index), then ponded water and finally via transpiration from the root zone, based on an average water content in the root zone.
However, evapotranspiration reduces the water content in the root zone, creating unsaturated zone storage. The minimum water content in the root zone is the wilting point water content, but this can only occur when the water table is below the root zone. The Two-Layer Water Balance evapotranspiration method requires time series for the root depth and the leaf area index, as well as the Reference evapotranspiration.
Some components of hydrological circle were not possible to be modelled. The river network (MIKE 11) with the chainage locations and cross-section points is necessary for the flow exchange with Groundwater Flow Module and for simulation of the lateral flow. The input data to the MIKE SHE model include data on topography, land use, geology, hydrogeology and meteorology. Preparation of the input data in suitable format was the first step needed to be done for the model setup. With MIKE Zero and MIKE SHE Tools standard file formats were converted to a format readable by the model.

Model domain and grid
Global Digital Elevation Model (GDEM) released by NASA and the Ministry of Economy, Trade, and Industry (METI) of Japan were used for geomorphological analyses. The GDEM data are posted on a 1 arc -second (approximately 30 m at the equator) grid and referenced to the 1984 World Geodetic System (WGS84)/1996 Earth Gravitational Model (EGM96) geoid. For the Morača river basin, tiles ASTGM2_N42E18, ASTGM2_N42E19, ASTGM2_N43E18 and ASTGM2_N43E19 were downloaded from http://earthexplorer.usgs.gov/.
Generation of catchment area was carried out using the basin command inside of the hydrology toolset in Arc GIS. Polygon with the catchment area information defines the model domain. All other spatial data, such as topography, are interpolated during pre-processing to that domain.

Topography
In MIKE SHE, the topography defines the upper boundary of the model and the top elevation of both the Unsaturated and the Saturated Zone model. It also defines the drainage surface for overland flow. Since the MIKE SHE has limitations regarding the number of computational cells, the finest topography grid resolution that was possible to use is 100×100 m. The Bilinear Interpolation method was applied for interpolation of the gridded DEM data (Fig. 5) [7].

Climate data
The daily accumulated precipitation amounts were used for modelling of the climate data. Since the Snow Accumulation and Melting Model are included daily mean temperatures were also imported. Most synoptic and climatic stations have complete data sets for the period from January 2000 to December 2010. For that period lack of precipitation data is quite present in time series obtained from precipitation stations. In order to fill these gaps, correlation analysis was performed Time series with serious data deficiency or poor correlation with neighbouring stations were excluded from further work. As opposite to precipitations, temperatures show strong correlation; therefore missing data were replaced and used for further modelling.
In mountainous areas, precipitation and temperature can vary significantly with altitude. However, there are rarely enough gauging stations to measure their spatial variability. Since the station-based rainfall data were used, choosing to correct the precipitation and temperature for elevation allows us to define a spatially variable correction factor.
When snow melt is included, various snow melt parameters should be specified to present that process: threshold melting temperature, degree-day coefficient that defines the amount of the snowmelt runoff, melting coefficient for solar radiation, minimum snow thickness that covers the entire cell with snow, maximum wet snow fraction in snow storage, initial total snow storage and initial wet snow fraction.
Snow accumulates until mean daily temperature is lower than threshold melting temperature defined in the model, otherwise the snow melts. That process involves conversion from dry snow to wet snow, and then to overland flow when the parameter that represents maximum wet snow fraction in snow storage is exceeded. Spatial distribution of the precipitation and temperature is done using the Thiessen polygons method. Two different dispositions were constructed, one for the precipitation and the other for temperature stations (Fig. 6). The Reference evapotranspiration is multiplied by the Crop Coefficient to obtain the Crop Reference Evapotranspiration. An assumption that evapotranspiration is uniformly distributed over the catchment was made because the spatial distribution of evapotranspiration was not available.
The distribution of the vegetation over the catchment area is necessary for the simulation of canopy interception and the process of evapotranspiration (Fig. 7). It can be defined over the model with the Leaf Area Index (LAI) and Root Depth (RD). The spatial distribution of the forest, shrub/trees and crops/grass is obtained from Global Land Cover 2000 Project (GLC 2000). The time series of the root depth and leaf area index for defined vegetation types, for either one year or for the growing season, are taken from Vegetation Database included in the model [8].

Model calibration and validation
Since MIKE SHE is not coupled with any river flow model it was not possible to simulate river discharge and compare it with discharge registered at the Podgorica hydrological station. Thus, validation of the model was possible by comparison between the simulated annual evapotranspiration and the one obtained from similar studies. The expected value of total annual evapotranspiration, expressed as precipitation ratio should be around 30 %.
It is also expected that inter annual distribution corresponds to the one in the nature. In case of reaching accumulation of the snow from October till the end of April, that will also be an indicator that the water balance is correct. The model was initially run for a ten year period. The initial run was done to achieve the approximate expected actual evapotranspiration percentage, the proper annual snow distribution and the minimal water balance error. Ten years of data were split into two halves for calibration and validation. The period from October 2000 till October 2005 was used for calibration, and the period from October 2005 till October 2010 for validation period. The model parameters represent the physical or hydrologic characteristics of a watershed. To adequately represent the system being modelled, during calibration, an iterative process, model parameters are set within an appropriate range. During that process, each model parameter is varied following a trial and error procedure, with all other parameters being constant (Tab. 1).
Using parameters fine-tuned in the calibration process, the model was validated with the data from the first five years. The final parameter list is presented in Tab. 2. Technical Gazette 24, Suppl. 2(2017), 427-434

Water balance analysis
The purpose of the water balance component of this assessment is to characterize the climatic, surface water, and groundwater components of the watershed. The water balance is intended for use as a screening tool to further evaluate water resources allocations within the watershed and to identify water balance components that may require further analysis during the next levels of watershed planning. The main analytical formulation of annual water balance for large areas is: Where: P -represents precipitation, ETevapotranspiration, RO -runoff.
ΔS change of water resources, which can be neglected if we analyze water balance for a long time period. But if we need to analyze the water budget on catchment scale, the equation becomes more complex: Where: I -plant interception, U -base flow, L -leakage to or from the catchment.
Since the river flow module is not included, the validation of the model was done comparing simulated annual evapotranspiration with the value obtained in earlier studies. In case of reaching accumulation of the snow from October till the end of April, that will be also an indicator that the water balance is correct.
The annual water balance summaries for the calibration and validation periods are presented in Tab. 3 and Tab. 4. The runoff (RO) is presented as the discharge observed at the Podgorica station. The year in the table represents a water year, the period of 12 months, the period between October 1 st of one year and September 30 of the next. The water year is designated by the calendar year in which it ends.
In the calibration period, the annual precipitation varies between 2480 and 4194 mm and from 1814 to 3519 mm in verification period. The discharge in the calibration period is in the interval from 1502 to 2988 mm, making the runoff coefficient within 60 to 86 % (Tabs. 5 and 6). The runoff coefficient in the verification period is around 74 %. Unfortunately, some of the very important rainfall stations on the catchment have significant data gaps and they are excluded from this study. The existing data point to conclusion that there is much more precipitation fell over the catchment than the one calculated by the model, so we can assume that the runoff coefficient could be lower. The total calculated error is less than 0.3 %. Infiltration varies both spatially and temporally but in general it depends on soil moisture condition, soil type and type of vegetation. Earlier studies indicate that infiltration over the Morača catchment is large and can be up to 80 % of falling rain [9]. The model indicates that the infiltration accounts for 71.12 % of the total precipitation amount for the period 2001-2005. Plant transpiration, evaporation from intercepted water, ponded water and evaporation from snow are included in the calculation of evapotranspiration. The loss of rain water through evapotranspiration is 28.88 %. Those figures are 71.98 % and 28.02 % for the verification period.
Although adequate snow accumulation needs temperature data on higher altitudes, the model indicates that the formation of the snow cover on the Moraca catchment mostly starts in November, when temperatures fall under the threshold value of 0.°C for snow accumulation. The model includes melting due to solar radiation and energy in rain besides melting due to air temperature. The biggest snow accumulation is in January and February. With increase in temperature during the spring the process of snowmelt starts and can last till the end of May.
Snow cover on the highest peaks melts till June and July.
The seasonal water balances for calibration and verification period are given in Tabs. 5 and 6. Table 3 Annual water balance summaries for calibration period Seasonal distributions of the water balance components for both periods are respectively presented in Figs. 8 and 9. Soil type and soil moisture condition define the percentage of precipitation that will infiltrate to the ground. As we noticed before 70 % of the precipitation goes to infiltration, thereby the smallest amounts infiltrate to the ground in summer. Significant amounts of water infiltrate during the rest of the year.

Conclusion
Uncertainties in the estimation of the water balance components depend on the accuracy of the measured input data into the model: the one with temporal variability (e.g. rainfall) and the one with spatial variability (e.g. rainfall, vegetation, groundwater levels, soil hydraulic properties etc.). Developing a water balance or water budget for the Skadar lake basin was difficult due to the absence of adequate data and information. MIKE SHE model for the Skadar lake catchment area shows satisfactory results and could be improved in the future with new input data, if they become available meanwhile.
It is no longer acceptable to manage groundwater and surface water independently of one another. Advances in data collection and availability, as well as computer resources, have now made distributed, physics-based watershed modelling feasible in a wide range of applications. MIKE SHE is one of the few commercially available codes that has been widely used for integrated hydrologic modelling. MIKE SHE's process based framework allows each hydrologic process to be represented according to the problem needs at different spatial and temporal scales. This flexibility has allowed MIKE SHE to be applied at spatial scales ranging from single soil profiles, to the field scale, and up to the watershed scale. Furthermore, each process can be represented at different levels of complexity. MIKE SHE has a modern, windows-based user interface that includes advanced tools for water quality, parameter estimation and water budget analysis.