Original scientific paper Received: August 13, 2015 Accepted: September 24, 2015 Data preparation for groundwater modelling - Ljubljansko polje aquifer system Priprava podatkov za modeliranje podzemne vode -vodonosni sistem Ljubljansko polje Janja Vrzel1- 2 *, Nives Ogrinc2- 3, Goran Vižintin4 Ecological Engineering Institute, Ljubljanska 9, 2000 Maribor, Slovenia international Postgraduate School Jožef Stefan, Jamova 39, 1000 Ljubljana, Slovenia 3Jožef Stefan Institute, Jamova 39, 1000 Ljubljana, Slovenia 4University of Ljubljana, Faculty of Natural Sciences and Engineering, Aškerčeva 12, 1000 Ljubljana, Slovenia Corresponding author. E-mail: janja.vrzel@iei.si Abstract The paper presents data preparation for groundwater modelling of Ljubljansko polje aquifer system in central part of Slovenia. Data preparation is part of the conceptualisation, which is the first step for developing a regional steady state model. This model will provide a framework for the development of transient groundwater flow model for the whole Ljubljansko polje aquifer system. The FEFLOW software is used for the visualization of the groundwater flow model of the Ljubljansko polje aquifer system and its geometry is based on the geological composition and hydrological characteristics of the studied area. Already known fact that too precise model geometry leads to the instability of the model was proved in ours study case, therefore many deliberate and rational decisions should be made during the conceptualisation. The model will extend our knowledge about interactions between surface and groundwater. Key words: FEFLOW, GIS, modelling, Ljubljansko polje aquifer system Izvleček V prispevku je opisana priprava podatkov za modeliranje vodonosnega sistema Ljubljansko polje, ki se nahaja v osrednji Sloveniji. Priprava podatkov je del koncep-tualizacije obravnavanega območja - prvega koraka v pripravi regionalnega stacionarnega modela. Le-ta pa bo pozneje nadgrajen v nestacionarni hidrodinamski model vodonosnega sistema Ljubljansko polje. Model pripravljamo v računalniškem programu FEFLOW, definicija njegove geometrije pa temelji na geološki sestavi in hidroloških lastnostih vodonosnika. Potrdili smo že znano dejstvo, da preveč natančna geometrija modela povzroča njegovo nestabilnost pri umerjanju, zato je potrebno med konceptualizacije sprejeti dobro premišljene odločitve, s katerimi se model poenostavi do takšne mere, da še lahko pričakujemo zadovoljive rezultate. Nastajajoč model bo razširil naše znanje o interakciji med površinsko in podzemno vodo. Ključne besede: FEFLOW, GIS, modeliranje, vodonosni sistem Ljubljansko polje Introduction Hydrological cycle is one of the most important cycles for people and their environment, therefore a lot of attention, efforts and money were directed to control its processes through the history. Special attention should be paid to groundwater research, where it is important to know its sources, infiltration areas, its dynamics and vulnerability to possible contamination from industrial and agricultural activities, and to assure its sustainable development for future generations. Nowadays, new technologies involving different models enable us to better understand rather complex processes of the hydrologic cycle [1]. The model is defined as a simplified representation of a complex natural system with three different types of applications: (1) predictive, used as a framework for prediction of future conditions, (2) interpretative, used for studying system dynamics or organizing data, and (3) generic, used to analyse the hypothetical hydrogeologic systems [2]. A regional steady state model for the whole Ljubljansko polje aquifer system is under preparation, which is the main drinking water source for Slovenian capital city Ljubljana and its surroundings. The aquifer has also economic importance, since this aquifer is also an invaluable water resource for the Union brewery (Pivovarna Union, d. d.) located within a highly urbanized and industrialized area near the centre of Ljubljana. The groundwater model will be extended into a transient transport model describing the transport of stable isotopes of O and H and 3H activity. The transport of dissolved species in groundwater systems facilitate to assess potential hazards to the public, when contaminants are used as transport material. The vulnerability of aquifer to environmentally harmful substances was observed in previous investigations: higher concentrations of herbicide diklobenil were detected in drinking water in 2001 due to maintenance of cemeteries; and high concentrations of atrazine could still be detected in drinking water in 2003 due to intensive agriculture [3, 4]. Few models of groundwater flow in Ljubljansko polje aquifer system were already developed in the past [5-7]. The aim of our work is to upgrade the existing models of Ljubljansko polje system with finite elements numerical model able to represent Sava river with higher accuracy than any other numerical model build until now. Numerical modelling of a physical system requires very well understanding of the system behaviours. The present work is motivated by the collection and processing of huge required data and the necessity to generate a complete DTM (Digital Terrain Model) of the entire computational domain, including the wetted region of the Sava river channel. The reason for it is in the great influence of the Sava on the elevation of the water table in the Ljubljansko polje aquifer system with intergranular porosity promoting the exchange between groundwater and Sava through the stream bed. Because of the high resolution of the Sava's bed geometry treatment in the groundwater model more stable model is expect during the calibration. Background At the Ljubljansko polje aquifer system (Figure 1) extensive investigation has been conducted in the last fifty years [8-10], therefore many geological and hydrological data are available for the formulation of a good conceptual model. Mathematical model simulates a natural system indirectly using a set of governing differential equations that describe physical processes occurring in the system, and numerical methods yield approximate solutions to the governing equation through the discretization of space and time [2, 11]. Many software were written to solve the governing equations for the different spatial, geological and hydrological conditions, which are based on two numerical techniques usually: the finite difference method (FDM), e.g. MODFLOW, and finite element method (FEM), e.g. FEFLOW. The FDM approximates differential equations by a differential approach, while FEM approximates differential equations by an integral approach [12]. The FEM is more appropriate choice when a subject of a research is a calculation of a precise water table elevation. On the other hand, it is better to use FDM when a calculation of a precise groundwater flux is performed. In our study the finite element flow simulation system FEFLOW® v6.1 from MIKE Powered by Danish Hydraulic Institute (DHI) was used, which enables to use stable isotopes as conservative transport materials. FEFLOW is a software package for three dimensional (3D) visualization of groundwater flow, mass transfer and heat transport in porous media. [14] 3D FEM for fluid flow and solute transport was developed by Huyakorn and Wadssworth in 1985 [13], while it is based on the equations for conservation of mass, Darcy's law for the anisotropic medium and the time-dependent species transport equation [14, 15]. The FEFLOW has the following advantages: (1) an ability to treat irregular shapes of a complex model boundary; (2) to refine the density of the mesh around features (e.g. wells, rivers, border, etc); (3) to simulate a groundwater age and transport of stable O and H isotopes [14]. Numerical models are used in parallel with development of powerful and affordable computer systems. In our case a 64-bit workstation under Windows 7 with CPU INTEL Core i7 and graphics ATI FirePro W7000 was used as a platform. The coefficients of the governing equations are the parameters that described the properties, boundaries, and stress of the system [11]. The choice of parameters is known as parameterisation, which is generally the most time-consuming procedure during modelling. The following data were obtained from different institutions/companies: (1) the counter map of the hardpan depth from VO-KA; (2) depths of piezometric groundwater table and groundwater abstractions from 48 and 35 wells, respectively from the VO-KA and Pivovarna Union, d. d, brewery; (3) isotopic analyses from the Jožef Stefan Institute. The ArcMap and AutoCAD Map were used for data processing before using FEFLOW. All three software FEFLOW, ArcMap and AutoCAD Map are fully integrated allowing quite easy transfer of the data between them. ArcMap and AutoCAD Map are geographic information systems (GISs), which are used for the management, analysis, and display of geographic information. GIS supports several views for working with geo- graphic information: (1) the geodatabase view, (2) the geovisualization view, (3) the geoprocessing view [15]. Further, a coordinate system is extremely important when we are dealing with geographic information. It is recommended to define it already in the GIS. Our model is prepared in D48/GK coordinate system. However, the replacement of the current coordinate system (D48/GK) with the European Spatial Reference System (ESRS) is in progress in Slovenia. Its horizontal direction is defined with the European Terrestrial Reference System (ETRS89/TM) [16]. Therefore, an extremely caution is required, while using GIS data since they are archived in different coordinate systems. The direction of groundwater flow within the Ljubljansko polje aquifer system is toward southeast in general. However, groundwater mining has a large influence on the hydrody-namic in the aquifer [17]. Ljubljansko polje aquifer system recharges mainly through infiltration of precipitation and river water. A smaller proportion of groundwater is recharged by lateral underground inflow [17, 18] in the area where the hardpan is deep enough at Ljubljansko barje in south and Kamniško-Mengeško polje in the north [18]. However, the underground inflow is still not very well known therefore a study about it should be made in the future. Small water recharge from surrounding hills is present as well. The Sava river is one of the most important sources of groundwater in Ljubljansko polje aquifer system, although it does not recharge the aquifer along its entire length. The Sava river water inflows into the aquifer system between Tacen and Šentjakob [19], while the groundwater recharges the Sava downstream from Šentjakob. It is important to highlight that there are some locations where the Sava's bed lies on the hardpan, therefore no river-ground-water interactions are present. An important loss of groundwater is also water supply. There are four pumping stations in Ljubljansko polje: Kleče, Jarški prod, Hrastje and Šentvid. The largest is in Kleče, where 55 % of drinking water was pumped in 2003 [20]. Figure 1: Map of Ljubljansko polje. Information about geological composition of Ljubljansko polje was taken from previous studies and approximately 115 drilling reports from the Javno podjetje Vodovod-Kanalizaci-ja, d. o. o., (VO-KA) and Pivovarna Union, d. d. Different lithological terms in drilling reports were used in the last 70 years therefore it was first necessary to unify them. Golden Software Strater was used for this purpose. According to all lithological information, the aquifer was divided into five modelled layers presented in Table 1. The aquifer is composed of Holocene and Pleistocene fluvial sediments, while its basement is composed of impermeable rocks (a hardpan). The geometry of the hardpan is diverse, its outcrops are in Šentjakob, Črnuče and Tacen. On the other hand, it is located quite deep, even more than 100 m, in the central part of the aquifer. Three topographic data types, obtained from different institutions/companies, were used for generation the model surface: (1) Digital elevation model (DEM) from The Surveying and Mapping Authority of the Republic of Slovenia; (2) LiDAR data for Sava floodplain from Holding Slovenske elektrarne, d. o. o., (HSE, d. o. o.); and (3) 126 Sava's cross-sections from the City Municipality of Ljubljana (MOL). The precise surface generation of the aquifer is a demanding task and thus presented in the second part of the article. Table 1: Lithological composition of modelled layers (ML) Precise surface generation for modelling The most appropriate GIS data for FEFLOW were prepared in the SHP or DAT files. The model geometry: polygon of the border, lines of the left and right banks of the Sava river, and points of wells and piezometers was constructed using SHP files. The outcrops of the hardpan, the Ljubljanica and Kamniška Bistrica rivers were the main criteria when the border of the study area was defined. The lines were extrapolated from TIN file with the following procedure. Designing a grid is one of the critical steps in groundwater modelling. Although fine grids provide more accurate solutions, the grid should be generated wisely since the duration of calculation depends on number of finite elements. Therefore, it is reasonably to use fine grids where accurate solutions are expected and coarse girds where details are not important [12]. The aquifer domain was discretised by Triangle Mesh Generator, into 96 129 mesh elements per slice or 480 645 mesh element in total. Triangle is an extremely fast meshing tool for complex meshes. It has been developed by Jonathan Richard Shewchuk [21]. The mesh in floodplain domain was refined after the discretization. The reason is higher aquifer vulnerability to pollution at that area. The mesh structure for the Sava river bed has been refined until the mesh was small enough and consequently appropriate for interpolation of the river cross-sections data. The mesh is succinct around pumping wells and piezometers as well, since pumping wells have a great influence on the aquifer and the known elevation of the hydraulic head is important for the model calibration. Figure 2 shows the full map extent of the Ljubljansko No. ML Lithology Epoch 1 Gravel, sand, gravel with sand and silt Holocene 2 Conglomerate, clay, conglomerate with lens and clay 3 Gravel with clay, sand with clay, gravel with narrow lines of conglomerate, gravel with _sand and silt_ Pleistocene 4 Conglomerate is predominate 5 Gravel with sand and silt Figure 2: Two-dimensional view of the Ljubljansko polje model domain with a finite element mesh. polje with aquifer discretization using a finite element mesh. The model domain is greater than 70.85 km2 and it is bounded approximately by the latitudes 46.12° N and 46.08° N and the longitudes 14.43° E and 14.64° E. DEM, Light Detection and Ranging (LiDAR) and the river cross sections were needed for interpolation of a surface, which were obtained as bunch of XYZ, LAS and KOO files, respectively. A XYZ file stores x, y, and z coordinates as floating-point values, where each row represents a distinct point record. A LAS file is an industry-standard binary format for storing airborne LiDAR data [22]. A KOO file is used for storing coordinates in national coordinate system, this is Gauss-Kruger coordinate system in Slovenia. These ASCII files were transformed into SHP files. The transformation of the XYZ and the KOO files is easy, while transformation of LAS files is more demanding. The LAS files were transformed indirectly into the SHP file through Terrain and Raster files. The ArcCata-log is needed for creating the Terrain file and 3D Analysis Tools in ArcMap for transformation Terrain file into the raster file with its anxious resolution. Once the data were imported into the FEFLOW, the surface interpolation was provided using the following steps: (1) interpolation of the full model surface by DEM data (25 m x 25 m); (2) interpolation of the floodplain by LiDAR data (10 m x 10 m); and (3) interpolation of the river bed. AKIMA linear interpolation was used in the first and second steps, while Neighbourhood Relationship interpolation method was used in the third step. Some important properties in Parameter Association should be defined before using Neighbourhood Relationship in- terpolation method: (1) under Default link selection the line among which the interpolation will be done should be chosen, (2) the Snap distance should be at least similar to the largest distance between cross-sections. The distances between the Sava cross-sections range from approximately 20 m to 690 m along the flow path, thus the Snap distance 690 m was used in our case. To find the border between wet and dry region of the river bed is demanding and time-consuming process. For this purpose the 3D Analysis Tool was used for conversion of the Terrain file into TIN file (Triangle Irregular Network). Figure 3 shows the result of the river bed topography interpolation in FEFLOW. Further, nodes for single interpolation should be chosen carefully. To make it easier the SHP files (polygons of the floodplain and the wetted region of the river bed) were made, which enable fast selection of mesh items by map polygon. The Data Management Tool in ArcMap was used, which further enables a fast conversion of lines into a polygon. The counter map of the hardpan defined by hydrogeologists, was reconstructed using: (1) coordinates of later drilled wells that reached the hardpan, and (2) extrapolated outcrops of Figure 3: The Sava river bed terrain a) and in the model b). the hardpan around Ljubljansko polje from the digitalized geological map. Few interpretations of new datasets were made and finally the most reasonable geometry of the hardpan was selected. Conclusion The conceptual model of Ljubljansko polje aquifer system was constructed, which will be progressed into the regional groundwater flow model. The main advantage of conceptual model is the construction of very high resolution of model domain geometry, especially the Sava river topography, which was constructed with river cross-sections data. The hardpan is the basement of the five model layers. The existing counter hardpan map was updated according to later drilling reports and digitalized geological map. However, during a conceptualisation it is necessarily to accept the trade-off among realism, generality and precision [23]. The simplification of a natural system is another great problem for beginners, who do not understand the philosophy of modelling yet - less is more [2]. Extremely high resolution of our model geometry causes instability of the model as well. Therefore, the model geometry had to be reconstructed. When the model geometry was constructed the initial and boundary conditions were defined. In the next step the calibration of the steady state model of groundwater flow will be performed. Acknowledgment The model is under preparation as a part of the ongoing EU 7th Research Project - GLOBAQUA (Managing the effects of multiple stressors on aquatic ecosystem under water scarcity). The project is financially supported by the IAEA Project under Contact No. RER8016 TC entitled 'Using Environmental Isotopes for Evaluation of Streamwater/Groundwater Interactions in Selected Aquifer in the Danube Basin' and Environmental Social Found (KROP 2012). The authors would like to thank all those who provided the data or offered us any other help: Javno podjetje Vodovod-Kanalizacija, Ljubljana, The City Municipality of Ljubljana, Holding Slovenske elektrarne, d. o. o., Pivovarna Union, d. d., Slovenian Environmental Agency, The Surveying and Mapping Authority of the Republic of Slovenia, DHD, d. o. o., and Institute for Water of the Republic of Slovenia. Special thanks go to hydrogeological specialist of Ljubljansko polje aquifer system Mrs Branka Bračič Železnik, who helped us with her professional advices and experience. References [1] L4018/GWM I, Groundwater modelling. 2000. [2] Anderson, M. P., Woessner, W. W. (1992): Applied Groundwater Modelling: simulation of flow and ad-vective transport. San Diego: Elsevier. [3] Kmecl, V., Simončič, A., Sušin, J., Žnidaršič-Pongrac, V. (2005): Rodovitnost tal. Podtalnica Ljubljanskega polja; pp. 164-177. [4] Bračič Železnik, B., Kladnik, D., Rejec Brancelj, I., Smrekar, A. (2005): Mestna raba tal. Podtalnica Ljubljanskega polja; pp. 188-201. [5] Janža, M., Meglič, P., Šram, D. (2015): Numerical hydrological modelling, Final Report [online, cited 2015]. Available on:. [6] Vižintin, G., Souvent P., Veselič, M., Cencur Curk, B. (2009): Determination of urban groundwater pollution in alluvial aquifer using linked process models considering urban water cycle. Journal of Hydrology; 377 (3-4), pp. 261-273. [7] Mikulič, Z., Andjelov, M. (2008): Groundwater modelling. International Symposium on Groundwater Flow and Transport Modelling, Ljubljana. [8] Radinja, D. (1951): Sava na Ljubljanskem polju. Geografski vestnik; 23, pp. M. [9] Žlebnik, L. (1971): Pleistocen Kranjskega, Sorškega in Ljubljanskega polja. Geologija; 17, pp. 477-491. [10] Premru, U. (1983): Osnovna geološka karta SFRJ. Tolmač lista Ljubljana; L 33-66, pp. 111. [11] Konikow, L. F. (1999): Use of numerical models to simulate groundwater flow and transport. In: Yurtsever, Y. (Editor), Modelling. Environmental isotopes in the hydrological cycle: Principles and applications, IAEA. Vienna, pp. 75-115. [12] Faust, C. R., Mercer, J. W. (1980): Ground-Water Modeling: Numerical Models. Ground Water; 18 (4), pp. 395-409. [13] Huyakorn, P. S., Wadsworth, T. D. (1985): FLAMENCO: A three-dimensional finite element code for analyzing water flow and solute transport in saturated-unsatu-rated porous media. GeoTrans, Inc, Herdon. [14] MIKE Provided by DHI [online, cited 2015]. Available on:. [15] ESRI[online, cited 2015]. Available on:. [16] Kete, P., Berk, S. (2012): Stari in novi državni koordinatni system v Republiki Sloveniji ter koordinatni system zveze Nato. Geoprostorska podpora obrambnemu sistemu Republike Slovenije; pp. 259-279. [17] Auersperger, P., Čenčur Curk., B., Jamnik, B., Janža, M., Kus, J., Prestor, J., Urbanc, J. (2005): Dinamika podzemne vode. Podtalnica Ljubljanskega polja; pp. 39-61. [18] Cerar, S., Urbanc, J. (2013): Carbonate Chemistry and isotope characteristics of groundwater of Ljubljans- ko polje and Ljubljansko barje aquifers in Slovenia. The Scientific World Journal; 2013, pp. 1-11. [19] Andjelov, M., Bat, M., Frantar, P., Mikulič, Z., Savić, V., Uhan, J. (2005): Pregled elementov vodne bilance. Podtalnica Ljubljanskega polja, pp. 27-38. [20] Bračič Železnik, B., Jamnik, B. (2005): Javna oskrba s pitno vodo. Podtalnica Ljubljanskega polja; pp. 101-120. [21] FEFLOW Help [online, cited 2015]. Available on:. [22] ArcGIS Help 10.1[online, cited 2015]. Available on:. [23] Levins, R. (1966): The strategy of model-building in population biology. American Scintist; 54 (4), pp. 421-431.