Modeling of Hydrogeological Processes in the Zeta Plain, Catchment Area of Lake Skadar in Montenegro

Goran Sekulić1, Milan Radulović1


1 Faculty of Civil Engineering, University of Montenegro, 81000 Podgorica, Montenegro; E-mail: This e-mail address is being protected from spambots. You need JavaScript enabled to view it



The paper deals with the modeling of hydrogeological processes in the Zeta plain, one of the most important parts of the catchment area of Lake Skadar in Montenegro. This lake and its catchment area represents a very significant water resource in Montenegro. The conceptual hydrogeological model for the Zeta plain has been developed on the basis of the interpretation of hydrogeological characteristics of this area, which have been a subject of study for many years.  Results of the mathematical model have been used in assessing the groundwater discharge from the Zeta plain intergranular aquifer into the Lake Skadar. A calibration of the model showed that the obtained results are sufficiently acceptable and that the model can be a good basis for future hydrological forecasts.

Keywords: Lake Skadar, hydrogeological processes, groundwater mathematical model, karst.



Lake Skadar is a transboundary lake in the southeastern part of the Dinarides, which is divided between Montenegro (around 65%) and Albania (around 35%). The first task of all management plan activities for Lake Skadar would be to examine the possibility of maintaining controlled or planned water levels for the lake. The lake's water levels depend solely on natural processes. Water levels are the cause of frequent flooding and insufficient utilization of the lake's capacity for various purposes: tourism, use of water resources and potentials, biodiversity and ecology, and many other positive aspects. These aspects do not only pertain to the lake basin but also to the whole River Bojana area, partly due to the same reasons: flooding and other flood-related difficulties for inhabitants on both banks.

Controlled lake level can be achieved only if the transient water balance of the lake is known and if the Bojana River riverbed is regulated, so that an appropriate quantity of water can be evacuated in any given time, which would enable the lake level to be maintained at the planned elevation. Lake level monitoring is performed at several hydrometric stations set up along the shoreline, on both sides of the border.


Figure 1: Lake Skadar catchment area in the Montenegrin part (red line with river network (blue lines)).


An inflow of groundwater from the surrounding aquifers into Lake Skadar is an extremely complex and insufficiently studied issue (Sekulic 2015). A specific part of this issue is the intergranular aquifer of the Zeta plain. This aquifier, which is one of the richest reservoirs of groundwater in Europe, has been formed within alluvial and glaciofluvial sediments (the area of the aquifer covers over 200 km2, and the aquifer thickness is in the range of 30-100 m). The sediments lay over the limestone paleorelief which characterize a very permeable karst aquifer. The intergranular aquifer of the Zeta Plain represents a precious water resource, especially in hydrological minimum.

The Zeta-Skadar Valley is one of the largest tectonic depressions in the Dinarides. The lowest parts of the depression lie beneath the water of Lake Skadar. The bed of the lake is below sea level, i.e. it represents a crypto-depression (e.g. the sublacustrine Raduš spring discharges 66 m below sea level).

The Zeta-Skadar depression was formed due to the deviation of general strike of Dinaric structures. The Dinaric structures turned from the known regional strike of the Dinarides (northwest-southeast) to the west-east and southwest-northeast directions (Sekulic 2015).

Glaciofluvial sediments of Zeta Plain are characterized by intergranular porosity. The sediments are very permeable, composed of gravel, sand, sandy gravel and weakly cemented conglomerates. Alluvial sediments of the Morača River and Cijevna River are represented by sand on the surface and sandy gravel in the lower zones (below 3 meters).

The intergranular aquifer of the Zeta Plain is characterized by high permeability, with mean value of hydraulic conductivity of around 5.0 x 10-3 m/s, and with transmissivity of 1.79 x 10-2 m2/s. The hydraulic parameters of the aquifer are obtained on the basis of numerous pumping tests performed in wells distributed across the plain (Sekulic 2015). The vertical hydraulic conductivity (Kz) is assessed on the basis of horizontal hydraulic conductivity (Kx,y). For the vertical hydraulic conductivity a ten times lower value is assigned (Kz=Kx,y/10).

The intergranular aquifer of the Zeta Plain is recharged by the atmospheric waters, by a sinking of the Morača River and Cijevna River, and through subterranean inflow from the bottom karst aquifer. The intergranular aquifer is discharged mostly at the south edge of the plain where a large number of small streams flow into the lake.

The precipitation data from the climate station "Podgorica" have been collected, and potential evapotranspiration has been assessed by applying the Penman method (Table 1). During the summer months, when potential evapotranspiration exceeds precipitation, atmospheric water is not infiltrated into the ground at all.


Table 1: The mean monthly and annual precipitation (mm), evapotranspiration (mm), effective infiltration (mm and %)


Groundwater Mathematical Model of the Zeta Plain

The regional groundwater mathematical model for the wider area of the Zeta Plain has been primarily created to build a data base for solving numerous water management issues. Results of the mathematical model have been used in assessing the discharge from the Zeta Plain intergranular aquifer into the lake.

The regional groundwater model in the area of the Zeta Plain has been developed through the following steps:

  • Collecting and studying the available data;
  • Creating the Conceptual Hydrogeological Model (CHM);
  • Selecting an appropriate computer code (software) for groundwater modelling;
  • Additional on-field data collection;
  • Interpreting all data and preparing them to fit appropriate formats to be input into the model;
  • Developing the groundwater mathematical model,
  • Sensitivity analysis and model calibration;
  • Model verification.


Collecting and Studying the Available Data

A significant number of literature and unpublished materials, as well as the results of numerous field and laboratory measurements have been used for the purpose of developing the regional model.


The Conceptual Hydrogeological Model of the Zeta Plain region

The conceptual hydrogeological model has been developed on the basis of the interpretation of hydrogeological characteristics of this area, which have been a subject of study for many years (Mijovic et al. 2014, Frankl et al. 2016, Radulović and Sekulić 2015, Sekulić 2015, Radulović et al. 2016).

The geometry of the hydrogeological system has been schematized into two layers. The first layer is represented by clastic sediments (glaciofluvial, alluvial, terraced and limnoglacial sediments), and by impermeable sediments (Pl) which occur below the intergranular aquifer in the south part of the plain (Fig. 2). The thickness of the first layer varies, as shown in the Paleorelief Map (Fig. 3); however, the average thickness is around 40 m. The second layer is "cut" at 200 m below sea level (Fig. 2).

The first layer is schematized as an unconfined aquifer, with variable hydraulic conductivity. The values obtained by pumping tests have been used as initial values of hydraulic conductivities. The hydraulic conductivity of 10-3 m/s have been adopted for the area characterized by distribution of glaciofluvial, alluvial and terrace sediments, while the hydraulic conductivity values rate of 10-6 m/s has been adopted for the south part of the plain in the Lake Skadar area, where limnoglacial sediments are distributed (Frankl et al. 201). The data on hydraulic conductivity of the second layer are scarce. Conclusions can be drawn on the basis of the existing knowledge which indicates that the carbonate rocks are very permeable (distributed in the central, northern, eastern and western part), whilst clastic sediments (claystone, marlstone and sandstone; Pl) distributed in the costal part of the lake and below lake floor are almost completely impermeable.

The groundwater flow regime is conditioned by the precipitation regime, the water levels of the Morača River and the Cijevna River, the water level of Lake Skadar, and by the regime of groundwater abstraction.

The boundary conditions of the conceptual hydrogeological model are schematized by external (divides, lake, rivers at the entering sections) and the internal boundaries (rivers across the area of the model, rivers ("drains") of the southern edge, abstraction wells and natural springs). One of the aims of creating the regional model was to avoid "the artificial", i.e. hydraulic boundaries when creating detailed models of the future investigation sites (construction sites, protection zones of the abstraction wells, etc.) (Moore and Doherty 2006).

The boundaries of the second layer, i.e. of the karst aquifer, are mostly divides (watersheds) of the immediate catchment area of the Zeta Plain determined on the basis of geomorphological, geological and hydrogeological criteria. However, the boundaries of the model do not represent the complete basin of the Zeta Plain, which would include the complete basins of the Morača River, Cijevna River and Sitnica River, but only up to the points from which these rivers begin sinking in such a way that they have an influence on the flow pattern of the the Zeta Plain intergranular aquifer. Therefore, the north boundary for the Morača River is placed in the zone 40 km upstream of the lake, at the point where the Morača River starts to sink. Upstream from this point, according to the available data (Sekulic 2015), the Morača cannot influence the flow pattern of the Zeta Plain intergranular aquifer, since it has an impermeable bed. Following this principle the entering points are determined for the Cijevna River and Sitnica River. The karst aquifer is bordered by Lake Skadar to the south in Hot and Hum bays, as well as by Malo Blato bay, where the aquifer discharge takes place, while part of the karst aquifer water has the possibility to continue flowing below Lake Skadar. Besides the above mentioned external borders of the second layer, there are also internal boundaries represented mainly by the Morača River and Cijevna River. Moreover, the internal boundaries are represented by karst springs such as Ribnica Springs, Milješ Spring, Krvenica Spring and Straganica Spring. Therefore, the mentioned karst springs are also the points through which karst aquifer discharge may take place during the rainy period of the year (Frankl et al. 2016).

The external boundary of the first layer, i.e. the intergranular aquifer, is mainly a geological boundary with the karst aquifer from which the intergranular aquifer is also recharged. Lake Skadar represents the southern boundary of the intergranular aquifer. The internal boundaries are defined by the Morača River, Cijevna River, Sitnica River, Ribnica River and Urelja River, out of which only the Morača River does not dry up during the summer months. Besides the aforementioned rivers, the boundary of the first layer follows smaller streams along the southern edge, which are schematized as the "drains" of the intergranular aquifer. The artificial discharge of the aquifer takes place through many abstraction wells which one should also bear in mind when defining the interior boundary conditions.

As for the water balance constituents, the atmospheric water infiltration value has been relatively well known (recharge rate), which is assessed by subtracting the calculated evapotranspiration (ET) from the quantity of the measured precipitation (P). The data from the climatologic station "Podgorica" are used for defining the recharge rate for the first layer, and for the second layer the data from the "Orahovo" precipitation station are used. The recharge rate for the first layer, i.e. for the intergranular aquifer, is assessed by calculating evapotranspiration (Table 1). For the second layer the recharge rate of 70% of precipitation is adopted based on the results obtained in previous research in the neighboring karst terrains (Radulović et al. 2010, Radulović et al. 2015). The artificial discharge of the intergranular aquifer by abstraction of groundwater from drilled wells (Qexp) is also a parameter which can be approximately determined on the basis of the data from the records of the biggest users (Podgorica Water supply System, the Podgorica Aluminum Plant, vineyards, etc.). The other elements of the water balance are mostly obtained in the process of development and calibration of the mathematical model, and in this case they include: inflow from surface waters (Qdpv), outflow to surface waters (Qopv), discharge into the Lake (Qoj), subterranean recharge from the karst aquifer into the intergranular aquifer (Qok), subterranean discharge from the intergranular aquifer into the karst aquifer (Qdk), and artificial discharge by wells (Qexp). Thus, the water budget equation for the intergranular aquifer is:

P + Qdpv + Qdk = Qoj + Qok + Qexp + ET


Figure 2: Conceptual hydrogeological model of the Zeta Plain region (Radulović at al. 2010);
Legend: 1. Intergranular aquifer, 2. Karst aquifer, 3. Impermeable sediments, 4. Stream, 5. Spring, 6. Groundwater flow direction, 7. Upper layer, 8. Lower layer, 9. Precipitation, 10. Evapotranspiration, 11. Inflow from surface waters, 12. Subterranean Inflow from the karst aquifer to the intergranular aquifer, 13. Subterranean outflow from intergranular aquifer to the karst aquifer, 14. Outflow from intergranular aquifer to the lake, 15. Hydraulic conductivity for glaciofluvial sediments, 16. Hydraulic conductivity for limno-glacial sediments.


Selection of the appropriate computer code for creating the mathematical model

The groundwater mathematical (numerical) model is a system with n differential equations of groundwater flow and with as many unknowns. For resolving such a system of equations it is necessary to use software by which equations can be solved in a so-called numerical way (we used GV Vistas software). In this case, the Finite difference model was chosen, so that the differential equations are replaced by finite differences and solved using numerical methods. For the purpose of developing a groundwater mathematical model of the Zeta Plain area, the best known finite difference model (the MODFLOW) developed by the USGS (United States Geological Survey) (McDonald and Harbaugh 1988) has been selected. The MODFLOW belongs to a group of so-called "block-centered" models, the reason being that the equation of groundwater flow refers to the center of each block, i.e. cell. (Scanlon et al. 2002,McDonald and Harbaugh 1988), Pollock 1994).


Data processing and preparation of adequate file formats

All the data obtained through studying the documentation and through field work were processed and prepared in adequate formats suitable for direct entry into the model. Therefore, special digital files for each spatially variable parameter were created, as follows:

  • topography, i.e. the top of the first layer,
  • the bottom of the first layer, i.e. the top of the second layer (Paleorelief map),
  • recharge rates,
  • initial groundwater levels, hydraulic conductivity, effective porosity, specific yield, etc.

The data from the documentation have been prepared in digital formats, first within files containing the coordinates (X, Y) in the first two columns, and the examined parameter (elevation, initial GWL, etc.) in the third column. Files prepared in this way have been used for interpolation by using the kriging method, with appropriate subjective corrections. The final files cover the whole area of the layers for which the data have been entered. Figure 3 shows a visual display of the raster elevation file for the bottom of the first layer.


Figure 3: Digital elevation model of the bottom of the Zeta Plain intergranular aquifer (CORIN 2000, Scanlon 2002).


Creation of the groundwater mathematical model

The total modeled area covers the surface of about 593 km2. The Zeta Plain covers an area of around 297 km2. The area of the model has been discretized in a plan view by a grid with cells 200 x 200 m to 400 x 400 m. The total area of both layers has been divided in 10,600 active cells.

Creation of the mathematical model primarily comes down to converting a Conceptual Hydrogeological Model (Fig. 2) into a numerical model (Bredehoeft 2005).

After discretization of the area of the model, all the cells that were left outside the selected boundaries of the karst aquifer have been marked as inactive cells (Fig. 4).

As mentioned before, the model is divided into two layers. The first (upper) layer is represented by the intergranular aquifer so that it is defined as the unconfined layer. The second (lower) layer is represented by the karst aquifer, which is mostly the unconfined aquifer (Unconfined–T Varies). In the part where impermeable sediments (Pl) are present below the intergranular aquifer, inactive cells are defined for the second layer.

Absolute bottom and top elevations of the first and second layer have been defined for each cell, by using files prepared in advance. As the bottom elevation of the model, i.e. of the second layer, the elevation of -200 m has been selected.

Regarding the time discretization, the model has been marked as transient, and divided into 12 periods (Stress Periods), with a duration for each period of about 30 days (Period Length), and each period has been divided into 3 time steps (No. Time Steps).

Horizontal and vertical hydraulic conductivities have been entered for each cell of the model. Although the hydraulic conductivity varies across the area, in the regional scale the glaciofluvial and aluvial sediments could be analyzed as homogenous intergranular media. During the calibration of the model, the best match between measured and calculated groundwater levels is obtained with hydraulic conductivity of 9.3 x 10-3 m/s. The effective porosity value of 0.25 has also been specified. This regional model with homogenous hydraulic characteristics could be used only for regional analyses. From this model,detailed "telescopic" models could also be extracted which must be additionally calibrated.

The model boundary conditions have been specified to simulate, as realistically as possible, the boundary conditions of the Conceptual Hydrogeological Model. However, given the complexity of the natural conditions in general, there are no software packages which could map the real image from the field.

For simulating the behavior of the Morača, Cijevna, Sitnica, Ribnica and Urelja rivers, the so called River Package has been used with reference to the perennial mean monthly water levels from the hydrometric stations "Zlatica" and "Podgorica" on the Morača River and "Trgaj" on the Cijevna River (The Hydrometeorological Institute of the Republic of Montenegro).

For simulation of the Lake Skadar behavior, the so-called General Head Boundary Package, (GHB) has been used, which determined the perennial mean monthly water levels of the lake measured at the hydrometric station "Plavnica" (The Hydrometeorological Institute of the Republic of Montenegro).

The behavior of small rivers of the southern edge (the Plavnica, the Zetica, the Gostiljska River, the Mala Mrka, the Velika Mrka, the Raičevića Žalica, etc.), through which the Zeta Plain intergranular aquifer discharges, has been simulated using the Drain Package. (Mehl et al. 2006, Neville and Tonkin 2004). The same package was used for simulation of the karst springs in the investigated area.

The groundwater abstraction from the intergranular aquifer of the Zeta Plain has been simulated using the Well Package, and only those groups of wells which can significantly influence the flow pattern in the regional scale have been included.

The aquifer recharge from precipitation (effective infiltration) has been simulated by using the Recharge Package and the perennial mean monthly precipitation from the climatologic station "Podgorica" have been used for the area covered by the intergranular aquifer, and for the area covered by karst aquifer on the terrain surface the same data from the rainfall gauging station "Orahovo" have been used.


The Results of the Model

The groundwater levels of the Zeta Plain have been computed using the calibrated regional mathematical model. This model can be used for the regional water balance assessment of the intergranular aquifer. On the basis of this regional model the detailed "telescopic" models could be extracted, but these new models have to be additionally calibrated. These "telescopic" models could be used for different hydraulic analyses of the groundwater flow such as: hydraulic analyses for the needs of opening a new and extending the existing sources for water supply of settlements; hydraulic analyses of impact of the planned facilities on the levels of groundwater of the Zeta Plain; delineating the protection zone of the water sources in the area of the Zeta Plain; etc.


Figure 4: The mathematical model grid with the defined boundary conditions for the Layer 1 (upper layer) (Sekulic 2015)
Legend: 1. Cells defined by the River Package, 2. Cells defined by the General Head Boundary Package (GHB), 3. Cells defined by the Drain Package, 4. Cells defined by the Well Package, 5. Inactive cells


Given all the advantages of hydraulic analyses on the basis of mathematical modelling, one should still bear in mind that this is a model which does not offer the possibility of complete mapping of exceptionally complex natural conditions. Therefore, all future analyses based on this model should be carefully conducted.

Contours of groundwater level obtained by mathematical modelling for the intergranular aquifer of the Zeta Plain have been presented in Figure 5.


Figure 5: Calculated contours of groundwater levels obtained by mathematical modeling of the Zeta Plain intergranular aquifer (for December).


Inflow Into the Lake from the Zeta Plain Intergranular Aquifer

The Zeta Plain intergranular aquifer mostly discharges alongside its southern edge where a great number of shorter streams are formed, some of which are: the Plavnica, the Zetica, the Gostiljska River, the Svinješ, the Pjavnik, the Velika and Mala Mrka. Discharge measurements of all these streams are very difficult to conduct since run off is diffused and takes place on a wide section.

For the needs of assessing the water discharge to Lake Skadar from the northern coast, the results of the previously created mathematical model of groundwater flow in the area of the Zeta Plain (Radulović et al. 2010, Radulović et al. 2015) have been used. On the basis of the mathematical model, the obtained value of the mean annual water discharge to Lake Skadar from the intergranular aquifer is 11.62 m3/s. In the total water balance of Lake Skadar this amount of water cannot be ignored, especially since it is part of the baseflow and is present throughout the whole year.



The presented mathematical model is only one part of a complex system of groundwater flow in the basin of Lake Skadar. The procedures performed during the formation of the model have revealed all the complexities of the problem, and it has been concluded that research must be continued in this direction. Bearing in mind the complexity of defining the groundwater balance, the obtained results should be used with caution, since in cases like this it is always realistic to expect a certain degree of error. Model results can be combined to further stages of research especially in planning of future states of the water level in the lake and their management.



Bredehoeft, J., 2005. The conceptualization model problem—surprise. Hydrogeology journal, 13:37–46, ISSN: 1431-2174

CORIN (2000): Digital CORIN Land Cover maps, European Environment Agency EEA i Joint Research Centre JRC.

Frankl, A., Lenaerts T., Radusinović S., Spalevic V., Nyssen J., The regional geomorphology of Montenegro mapped using Land Surface Parameters, Zeitschrift für Geomorphologie, Volume 60, Number 1, March 2016, pp 21-34(14)

McDonald M. G., Harbaugh A. W., (1988): A modular three-dimensional finite-difference groundwater flow model: U.S. Geological survey Techniques of Water Resources Investigations, book 6, chap. A1, pp 586

Mehl, S. and M.C.Hill. (2010). Grid-size Dependence of Cauchy boundary conditions used to simulate Stream-Aqufer Interaction. Advances in Water Resources, ISSN: 0309-1708, 33, pp 430-442

Mehl, S., M. C. Hill, and S. A. Leake, 2006. Comparison of local grid refinement methods for MODFLOW, Groundwater 44(6): 792-796, ISSN: 1745-6584

Mijovic S., Sekulic G., Weaknesses in Skadar/Shkodra lake watershed assessment and management, International Journal of Ecosystems and Ecology Science ,(IJEES)Volume 4/1, pp 11-14, 2014, ISSN 2224-4980).

Moore, C., and Doherty, J. (2006). "The cost of uniqueness in groundwater model calibration." Advances in Water Resources, ISSN: 0309-1708, 29(4), 605-623.

Neville, C. J., and M. J. Tonkin, 2004. Modeling multiaquifer wells with MODFLOW. Groundwater 42(6):910-919, ISSN: 1745-6584

Pollock, D.W., 1994: User's Guide for MODPATH/MODPATH-PLOT, Version 3: A particle tracking post-processing package for MODFLOW, the U.S. Geological Survey finite-difference ground-water flow model: U.S. Geological Survey Open-File Report 94-464, pp 234

Radulović M. M, Stevanović Z., Radulović M. (2010): First outcomes from new approach in assessing recharge of highly karstified terrains – cases examples from Montenegro. In: Andreo B, Carrasco F, Durán J. J, LaMoreaux J.W (ed) Advances in Research in Karst Media, Springer, Berlin, pp 25-31

Radulović M., Sekulić G., Determining locations of sublacustrine springs by remote sensing: the Skadar Lake case example of Montenegro, Строительство уникальных зданий и сооружений. ISSN 2304-6295. 2 (29). 2015. pp 73-86,

Radulović M.M., Radulović M., Stevanović Z., Sekulić G., Radulović V., Burić M., Novaković D., Vako E., Blagojević M., Dević N., Radojević D. (2015) Hydrogeology of the Skadar Lake basin (Southeast Dinarides) with an assessment of considerable subterranean inflow. Environmental Earth Sciences, 74:71–82, DOI 10.1007/s12665-015-4090-7.

Radulović M.M., Sekulić G, Blagojević M., Krstajić J., Vako E. (2016) An assessment of territory participation in transboundary karst aquifer recharge: A case study from the Skadar Lake catchment area. In: Stevanović Z., Krešić N., Kukurić N. (eds) Karst without Boundaries, Internationan Association of Hydrogeologist, CRC Press, London, pp 87-100

Scanlon B.R., Healy R.W., Cook P.G., Choosing appropriate technuiqeus for quantifying groundwater recharge, Hydrogeology Journal 10, pp.18-39, 2002.

Sekulić G., (project leader), group of authors, Development of hydrological study of regulation of Skadar Lake and Bojana river water regime, Volume I and Volume II, Montenegrian Academy of Science and Arts, Special editions (Monographies and Studies), Knjiga 111, ISBN 978-86-7215-347-7 (Vol I), ISBN 978-86-7215-348-4 (Vol. II), Podgorica 2015