Development of a Hydrogeological Conceptual Model of the Varaždin Alluvial Aquifer

: The Varaždin aquifer represents the main source of water for public supply, agricultural, and industrial purposes in the Varaždin County in NW Croatia. In the last decades, this area has experienced contamination of groundwater with nitrates. This study describes the conceptualization of the Varaždin aquifer for the purpose of developing numerical model of groundwater ﬂow and nitrate transport. Within the study, three important elements are deﬁned: aquifer geometry, recharge from precipitation, and other boundary conditions. 3D aquifer model revealed that Varaždin aquifer consist of three layers: upper aquifer, semipermeable interlayer, and lower aquifer. The Wetspass-M model was used for the assessment of spatial and temporal distribution of water balance components for the period 2008–2017. Results of the model indicate that the average annual precipitation is distributed as 34% groundwater recharge, 21% surface runoff, and 45% actual evapotranspiration. The maps of equipotential lines show the behavior of the aquifer system and deﬁne boundary conditions, i.e., recharge and discharge areas of the aquifer: an inﬂow boundary from Drava River and accumulation lake Varaždin on the northwest and north, no ﬂow boundary on the west and south, and an outﬂow boundary on the east.


Introduction
The Varaždin aquifer is a vital source of water for public supply, agricultural, and industrial purposes in the Varaždin County in NW Croatia. Moreover, according to its hydrogeological characteristics, it represents one of the strategic groundwater resources in Croatia. In the last few decades, this area has experienced high nitrate concentrations caused by anthropogenic sources, such as manure, synthetic fertilizers, septic systems, and other wastewaters. The contamination of groundwater with nitrates have caused the shutting down of the pumping site Varaždin. This paper is part of a broader study being conducted in Varaždin alluvial aquifer with the aim of assessing the origin, fate, and the transport of nitrate within the study area.
Conceptual model is a simplified representation of a groundwater system and is based on geological, geophysical, hydrological, and hydrogeochemical information [1]. The development of an appropriate conceptual model is one of the most important steps in any successful modelling study [2,3]. A detailed survey of hydrodynamic characteristics of the aquifer was undertaken to develop hydrogeological conceptual model, which will be used as a foundation for setting up a numerical groundwater flow and nitrate transport model. In the present study, preparation of conceptual model involved identification of the study area, creation of 3D model of hydrogeological system, estimation of recharge from precipitation, and defining boundary conditions are shown.
The sustainability and efficient management of groundwater reserves in the aquifer relies on groundwater recharge. According to Healy [4] groundwater recharge can be classified into two categories: focused recharge from surface water bodies such as rivers, canals, and lakes, and the diffuse recharge from infiltration of precipitation through the unsaturated zone to the groundwater. The analysis of groundwater and surface water levels allows defining the boundary conditions as well as the relationship between surface water and groundwater. However, quantifying diffuse recharge from precipitation is of particular importance to this study. Over the past decade, Wetspass model has been used in different parts of the world to calculate water balance components, including groundwater recharge. Gebreyohannes et al. [5] developed Wetspass model to assess water resources in Geba basin in Ethiopia. Porretta-Brandyk et al. [6] evaluated and verified Wetspass model with focus on river runoff modeling in rural catchments in Poland. Zarei et al. [7] used Westpass model for assessment of groundwater recharge, surface runoff, and evapotranspiration in different land-use types in northeast Iran. Zhang et al. [8] addressed the effects of urbanization on water balance components in Beijing, China by using Wetspass model. Salem et al. [9] applied the Wetspass model to assess the water balance components in the Drava basin in Hungary. Previous assessments of groundwater recharge from precipitation and its spatial distribution in the study area are rare and poorly understood. Patrčević [10] used the experimental station to analyze the vertical water balance of groundwater in the Drava alluvium and estimated that recharge accounts for 38.3% of the total precipitation. Larva [11] developed a numerical model of the Varaždin aquifer to predict future nitrate concentration depending on the abstraction rates on pumping sites. The author assigned the recharge as a share of average annual precipitation, depending on covering layer thickness. A value of 35% was used in the area where covering layer does not exceed thickness of 2.5 m, while a value of 20% was applied in the area with covering layer thickness above 2.5 m.
In this study, the improved Wetspass-M model was used to explore the relationship between precipitation and recharge. The selection of the software was based on the data availability, and insights from previous investigations [9,12] where authors recommended using the Wetspass-M model for groundwater recharge assessment in developing groundwater flow models for the Drava basin. The presented research is the first study to evaluate the spatial distribution of long-term average groundwater recharge from precipitation in the Varaždin aquifer, and this information will serve as important input for developing numerical model, together with aquifer geometry and other boundary conditions.

Study Area
The study area is situated in the Drava River valley in the northwestern part of Croatia ( Figure 1). It covers the part of Varaždin alluvial aquifer upstream of the town of Varaždin where highest nitrate concentrations were observed, with an area of approximately 200 km 2 . The Varaždin alluvial aquifer is mostly unconfined and represented by Pleistocene and Holocene alluvial deposits [13], which are mainly composed of gravel and sand with occasional lenses and interbeds of silt and clay [14]. There are two pumping sites in the study area: active-Vinokovšćak; and inactive-Varaždin. The boundaries of the study area are represented by the surface water on the northwestern and northern edge, Haloze hills on the western and Varaždinsko-Topličko gorje hills on the southern edge of the aquifer. On the northwest, Drava River flows into accumulation lake of the hydroelectric power plant Varaždin (HPP Varaždin), from which it continues in two paths: (1) as the Drava River watercourse in the north, and (2) through an intake channel for hydroelectric power plant Varaždin. The derivation channel carries the water from the power plant downstream. The second watercourse that runs through the area is Plitvica stream, located close to the southern edge. The climate of the study area is classified as a warm temperate climate (Cfb in the Köppen-Geiger climate classification system) [15]. The long-term (1981-2010) average annual amount of precipitation measured at the Varaždin meteorological station is 832 mm, although annual precipitation can be highly variable ranging from 481 to 1312 mm, with more precipitation occurring during the summer [16]. Precipitation generally originates from the Atlantic air masses, but Mediterranean influence is also notable, especially du ring the colder periods of the year [14]. The average annual temperature in Varaždin is 10.6 °C. The warmest month is July, with an average temperature of 20.9 °C. The coldest month is January, with an average temperature of 0 °C.

Aquifer Geometry
To achieve better understanding of the aquifer system of the study area, the presented methodology was followed ( Figure 2). The first phase was to collect all available geological data in the study area, including existing maps, borehole logs, and cross-sections in order to define boundaries of the aquifer. Geological data used within this study are mainly related to construction of the hydroelectric power plant Varaždin and the development of pumping sites in the Varaždin area. The data were collected from different sources, including the database of the Croatian Geological Survey, but also numerous technical reports from water utility, geotechnical and civil engineering companies. Horizontal characterization of the aquifer, i.e., model domain was determined by the Basic Geological Map of the Republic of Croatia at a scale of 1:100,000 for the Varaždin [17] and Čakovec sheet [18], and is represented by Quaternary alluvial deposits. The second phase included construction of contour maps based on digital elevation model (DEM) and borehole logs, in order to achieve vertical characterization of the aquifer. DEM was used to create the top surface of the model in ArcGIS software. Total 60 geological boreholes were The climate of the study area is classified as a warm temperate climate (Cfb in the Köppen-Geiger climate classification system) [15]. The long-term (1981-2010) average annual amount of precipitation measured at the Varaždin meteorological station is 832 mm, although annual precipitation can be highly variable ranging from 481 to 1312 mm, with more precipitation occurring during the summer [16]. Precipitation generally originates from the Atlantic air masses, but Mediterranean influence is also notable, especially du ring the colder periods of the year [14]. The average annual temperature in Varaždin is 10.6 • C. The warmest month is July, with an average temperature of 20.9 • C. The coldest month is January, with an average temperature of 0 • C.

Aquifer Geometry
To achieve better understanding of the aquifer system of the study area, the presented methodology was followed ( Figure 2). The first phase was to collect all available geological data in the study area, including existing maps, borehole logs, and cross-sections in order to define boundaries of the aquifer. Geological data used within this study are mainly related to construction of the hydroelectric power plant Varaždin and the development of pumping sites in the Varaždin area. The data were collected from different sources, including the database of the Croatian Geological Survey, but also numerous technical reports from water utility, geotechnical and civil engineering companies. Horizontal characterization of the aquifer, i.e., model domain was determined by the Basic Geological Map of the Republic of Croatia at a scale of 1:100,000 for the Varaždin [17] andČakovec sheet [18], and is represented by Quaternary alluvial deposits. The second phase included construction of contour maps based on digital elevation model (DEM) and borehole logs, in order to achieve vertical characterization of the aquifer. DEM was used to create the top surface of the model in ArcGIS software. Total 60 geological boreholes were identified within the study area ( Figure 1). Analysis of the borehole logs and cross-sections indicates that subsurface hydrostratigraphy consists of covering layer, upper aquifer, semipermeable interlayer, lower aquifer, and aquifer bottom. The borehole coordinates and information about depths of the covering layer, semipermeable interlayer top and bottom, and aquifer bottom were prepared in spreadsheet, which was used as an input data for construction of contour maps in Surfer software, using Kriging interpolation method. In the third phase, the contour maps were imported into Visual MODFLOW Flex software and 3D aquifer model of the study area was built. identified within the study area ( Figure 1). Analysis of the borehole logs and cross-sections indicates that subsurface hydrostratigraphy consists of covering layer, upper aquifer, semipermeable interlayer, lower aquifer, and aquifer bottom. The borehole coordinates and information about depths of the covering layer, semipermeable interlayer top and bottom, and aquifer bottom were prepared in spreadsheet, which was used as an input data for construction of contour maps in Surfer software, using Kriging interpolation method. In the third phase, the contour maps were imported into Visual MODFLOW Flex software and 3D aquifer model of the study area was built.

Groundwater Recharge from Precipitation
WetSpass software is a GIS-based quasi steady state spatially distributed water balance model, which was developed for estimation of the long-term average spatial patterns of surface runoff, actual evapotranspiration, and groundwater recharge [19,20]. The newer version of the model, Wetspass-M, has ability to use monthly input data and compute monthly water balance components. This way the user can assess water balance components on monthly, seasonal or yearly basis. The Wetspass-M model uses spatially distributed input parameters, including DEM, slope, land use, soil type, groundwater level, and meteorological data (precipitation, number of rainy days per month, air temperature, wind speed, and potential evapotranspiration) [20][21][22]. All input data were prepared in the form of grid maps in ArcGIS software and are presented in the official coordinate system of the Republic of Croatia (HTRS96/TM). DEM of the study area is displayed with 20 m resolution. The highest point of the study area is 266 m a.s.l. at Haloze hills in the west and the lowest point is 166 m a.s.l. in the east ( Figure 3a). The resolution of all other maps were based on DEM resolution with cell size 20 × 20 m, 948 columns, and 875 rows. The slope map was created from DEM using spatial analyst tool in ArcGIS. Land use and soil types are connected to respective maps through lookup tables. Land use data for the Varaždin area were obtained from the CORINE database for Land Cover (CLC2018), GIS vector layer available online at https://land.copernicus.eu/pan-european/corine-landcover/clc2018. Total 14 CLC classes were identified in the study area, which have been reclassified into 11 land use classes for Wetspass-M input data purpose (Figure 3b), to be suitable with land use lookup table. Around 77% of the total study area is agricultural land. The other 23% is divided between build up, including city center and open build up (10%), forest, including deciduous, coniferous and mixed forest (9%), shrub (2%), and water bodies, including lake and navigable river (2%). The soil map was constructed using the combination of Thiessen polygon method in ArcGIS for 60 geological boreholes ( Figure 1) and pedological map of the Republic of Croatia, especially around the Plitvica stream in the southern part of the study area where there are not many boreholes. Pedological map is obtained from ENVI environmental atlas, available online at http://envi.azo.hr/. The most common types of soil in the study area are loam, silty clay, sand, sandy clay, and clay, covering 49%, 27%, 12%, 6%, and 6% of the area, respectively (Figure 3c).

Groundwater Recharge from Precipitation
WetSpass software is a GIS-based quasi steady state spatially distributed water balance model, which was developed for estimation of the long-term average spatial patterns of surface runoff, actual evapotranspiration, and groundwater recharge [19,20]. The newer version of the model, Wetspass-M, has ability to use monthly input data and compute monthly water balance components. This way the user can assess water balance components on monthly, seasonal or yearly basis. The Wetspass-M model uses spatially distributed input parameters, including DEM, slope, land use, soil type, groundwater level, and meteorological data (precipitation, number of rainy days per month, air temperature, wind speed, and potential evapotranspiration) [20][21][22]. All input data were prepared in the form of grid maps in ArcGIS software and are presented in the official coordinate system of the Republic of Croatia (HTRS96/TM). DEM of the study area is displayed with 20 m resolution. The highest point of the study area is 266 m a.s.l. at Haloze hills in the west and the lowest point is 166 m a.s.l. in the east (Figure 3a). The resolution of all other maps were based on DEM resolution with cell size 20 × 20 m, 948 columns, and 875 rows. The slope map was created from DEM using spatial analyst tool in ArcGIS. Land use and soil types are connected to respective maps through lookup tables. Land use data for the Varaždin area were obtained from the CORINE database for Land Cover (CLC2018), GIS vector layer available online at https://land.copernicus.eu/pan-european/corine-land-cover/clc2018. Total 14 CLC classes were identified in the study area, which have been reclassified into 11 land use classes for Wetspass-M input data purpose (Figure 3b), to be suitable with land use lookup table. Around 77% of the total study area is agricultural land. The other 23% is divided between build up, including city center and open build up (10%), forest, including deciduous, coniferous and mixed forest (9%), shrub (2%), and water bodies, including lake and navigable river (2%). The soil map was constructed using the combination of Thiessen polygon method in ArcGIS for 60 geological boreholes ( Figure 1) and pedological map of the Republic of Croatia, especially around the Plitvica stream in the southern part of the study area where there are not many boreholes. Pedological map is obtained from ENVI environmental atlas, available online at http://envi.azo.hr/. The most common types of soil in the study area are loam, silty clay, sand, sandy clay, and clay, covering 49%, 27%, 12%, 6%, and 6% of the area, respectively (Figure 3c).  Groundwater level data and meteorological data have been provided by the Croatian Meteorological and Hydrological Service (DHMZ) for the study period 2008-2017. In total, 34 observation wells ( Figure 1) were used for construction of groundwater level map. Groundwater level data were analyzed in detail for the study period, then hydrological condition of medium groundwater level was selected, and finally groundwater level map was produced using Kriging interpolation method in Surfer software (Figure 3d).
For the purpose of this study, daily meteorological data from Varaždin meteorological station were used as it is located in vicinity of the study area (46.28278 N, 16.36389 E, 167 m a.s.l.) ( Figure 1) and has all the necessary data records. The study area is fairly flat and Varaždin meteorological station is in vicinity so we assumed that station as a representative for the area. Meteorological input data (precipitation, number of rainy days per month, air temperature, wind speed, potential evapotranspiration) were prepared as monthly average values for the study period. The average annual precipitation in the study period is 912 mm, with about 70% of it being concentrated from May to November. The average number of rainy days per month varies between 7 and 11 days. The mean annual temperature is 11.5 °C, with July as the warmest month (average temperature 21.9 °C) and January as the coldest month (average temperature 0.7 °C). Average annual wind speed measured at height of 2 meters above soil surface is 2.4 m/s. Evapotranspiration represents the sum of water loss through the process of rainfall interception and transpiration from plants, and evaporation from soil surfaces. It is commonly calculated from climatological data because of the difficulty to obtain accurate field measurements. The evaporation power of the atmosphere is expressed by the evapotranspiration from the reference surface not short of water. That is so-called reference crop evapotranspiration or reference evapotranspiration, denoted as ET0. The standardized reference surface is a hypothetical grass reference crop with specific characteristics [23]. The Groundwater level data and meteorological data have been provided by the Croatian Meteorological and Hydrological Service (DHMZ) for the study period 2008-2017. In total, 34 observation wells ( Figure 1) were used for construction of groundwater level map. Groundwater level data were analyzed in detail for the study period, then hydrological condition of medium groundwater level was selected, and finally groundwater level map was produced using Kriging interpolation method in Surfer software (Figure 3d).
For the purpose of this study, daily meteorological data from Varaždin meteorological station were used as it is located in vicinity of the study area (46.28278 N, 16.36389 E, 167 m a.s.l.) ( Figure 1) and has all the necessary data records. The study area is fairly flat and Varaždin meteorological station is in vicinity so we assumed that station as a representative for the area. Meteorological input data (precipitation, number of rainy days per month, air temperature, wind speed, potential evapotranspiration) were prepared as monthly average values for the study period. The average annual precipitation in the study period is 912 mm, with about 70% of it being concentrated from May to November. The average number of rainy days per month varies between 7 and 11 days. The mean annual temperature is 11.5 • C, with July as the warmest month (average temperature 21.9 • C) and January as the coldest month (average temperature 0.7 • C). Average annual wind speed measured at height of 2 meters above soil surface is 2.4 m/s. Evapotranspiration represents the sum of water loss through the process of rainfall interception and transpiration from plants, and evaporation from soil surfaces. It is commonly calculated from climatological data because of the difficulty to obtain accurate field measurements. The evaporation power of the atmosphere is expressed by the evapotranspiration from the reference surface not short of water. That is so-called reference crop evapotranspiration or reference evapotranspiration, denoted as ET 0 . The standardized reference surface is a hypothetical grass reference crop with specific characteristics [23].
The crop evapotranspiration under standard conditions (ET C ) refers to the evapotranspiration from excellently managed, large, well-watered crop fields that achieve full production under the given climatic conditions. The ET C is determined by crop coefficients (K C ) that relate ET C to ET 0 [23].
A large number of empirical or semi-empirical equations have been developed for assessing crop or reference crop evapotranspiration from meteorological records. The Food & Agriculture Organization (FAO) of the United Nations Penman-Monteith (FAO-56 PM) method [23] is recommended as the international standard method for the definition and computation of the potential reference evapotranspiration (ET 0 ). The FAO-56 PM method requires radiation, air temperature, air humidity, and wind speed data. The FAO-56 PM enjoys worldwide adoption as the most accurate [24], but the number of requested climatic variables usually makes its application questionable. As a result, many articles deal with its comparison to other proposed methods in order to avoid that much meteorological variables [24][25][26][27]. Our study had proper meteorological records and the FAO-56 PM was used as a reference method to obtain ET 0 according to the equation where ET 0 is the reference crop evapotranspiration (mm day −1 ); R n is the net radiation (MJ m −2 day −1 ); G is the soil heat flux (MJ m −2 day −1 ), which is regarded as null for daily periods; T a is the average daily air temperature at a height of 2 m ( • C); u is the wind speed at a height of 2 m (m s −1 ); e s is the saturation vapor pressure (kPa); e a is the actual vapor pressure (kPa); e s − e a is the vapor pressure deficit (kPa); ∆ is the slope of the saturation vapor pressure-temperature curve (kPa • C −1 ); and γ is the psychrometric constant (kPa • C −1 ).
The daily records of minimum and maximum temperature ( • C), mean relative humidity (%), wind speed (m/s) and actual duration of sunshine in a day (h/day) were collected for the study period from the DHMZ database. The data were transformed into appropriate format to be used within "ET 0 calculator" software [28] in order to calculate daily reference evapotranspiration for Varaždin meteorological station. Calculated daily data of ET 0 were summed into monthly data for the future steps of the calculation.
Nistor et al. [29][30][31] assessed the relationship of land cover data to the crop evapotranspiration based on seasonal potential evapotranspiration and standard seasonal crop coefficients (K c ) presented in the FAO Paper no. 56 [23]. The authors used seasonal land cover coefficients K CLC to the reference ET 0 to carry out the crop evapotranspiration ET C as showed by the equation where ET C is crop evapotranspiration, ET 0 is potential reference evapotranspiration, and K CLC is land cover coefficient. Nistor & Porumb-Ghiurco [29] and Nistor et al. [30] distinguished four stages of functionality of crops: initial (March-May), mid-season (June-August), end-season or late season (September and October), and cold season (January, February, November, and December) and they assigned K CLC values to each land cover class depending on the season. In order to assess the most possible realistic values of ET C from previously calculated ET 0 over the Varaždin aquifer area from Equation 2 was used. The CORINE Land Cover 2018 database (CLC2018) was used to obtain a land cover map of the study area where 14 classes of land cover classes were identified. The seasonal coefficients of KCLC [29,30] were assigned to each present CLC2018 class and used to calculate their monthly crop i.e., land cover evapotranspiration in the study period.

Boundary Conditions
The general behavior of the aquifer system and its boundary conditions can be well described by constructing map of water table contours lines or equipotential lines. Data on the groundwater levels and the surface water levels for the period 2008-2017 were used to construct the equipotential lines. The available time series data of groundwater levels for the study period included measurements from 34 observation wells in the study area (Figure 1), which are part of a national monitoring network. The measurements are performed by DHMZ every 3-4 days. The available time series data of surface water levels included measurements on four hydrological stations. The daily measurements of water level at inflow and outflow of the accumulation lake and downstream of the hydroelectric power plant Varaždin are provided by Croatian National Power Company (HEP), while daily measurements of water level for hydrological station Varaždin (where Drava River meets the derivation channel) are provided by DHMZ.
Due to the insufficient number of hydrological stations on the Drava River, it was necessary to create virtual hydrological stations between actual hydrological stations ( Figure 1). The water level at virtual hydrological stations were calculated by linear interpolation method between two hydrological stations with measurements of water level. Virtual hydrological stations were created at three sections at a distance of 1 km: (1) on the northwestern boundary of the aquifer between station at inflow to the accumulation lake and hydrological station Borl I that is located on the Drava River in Slovenia outside the study area; (2) on the Drava River between station at outflow of the accumulation lake and hydrological station Varaždin; (3) on the derivation channel between station downstream of the hydroelectric power plant Varaždin and hydrological station Varaždin. Water levels for hydrological station Borl I are available online at ARSO database (http://vode.arso.gov.si/ hidarhiv/pov_arhiv_tab.php?p_vodotok=Drava&p_postaja=2150). For the study period, the water level data for all monitoring stations, both measured and virtual, were analyzed in detail in Microsoft Excel and the hydrological conditions of low water levels and high water levels were selected. The results of the analysis were used as an input data for construction of equipotential lines for low and high water levels in Surfer software, using Kriging interpolation method.

Aquifer Geometry
The resulting 3D Varaždin aquifer model based on the contour maps constructed from borehole data is shown in the Figure 4. The model consists of three layers with different hydrogeological characteristics: upper aquifer, semipermeable interlayer, and lower aquifer. The covering layer of the aquifer is represented by low permeable silty-clay deposits. However, it is not continuously developed (Figure 4), with thickness from 0 to 5 m, which contributes to generally high vulnerability of the aquifer. Although heterogeneous, alluvial aquifer is composed mainly of gravel and sand with lenses and interlayers of silt and clay. The thickness of the aquifer increases from less than 5 m in the NW part to about 65 m in the SE part of the study area. A more significant semipermeable silty-clay interlayer appears in the east of the study area, at a depth of about 35 m. Its thickness is up to 5 m, but borehole data reveal that it is not continuously deposited, which indicates a hydraulic connection between aquifer sediments above and below. The aquifer bottom is composed of marl and sandstone in the west, while clay, silt, and marl are present in central and the eastern part below the aquifer. It is considered as being impermeable or having a no flow boundary in the vertical direction.

Groundwater Recharge from Precipitation
The results of WetSpass-M model are raster maps of mean monthly water balance components: actual evapotranspiration, interception, surface runoff, and groundwater recharge for the 2008-2017 period. Output raster maps of calculated monthly values were summed in ArcGIS software (total 12 maps for each water balance component) to obtain the spatial distribution of average annual values ( Figure 5).

Groundwater Recharge from Precipitation
The results of WetSpass-M model are raster maps of mean monthly water balance components: actual evapotranspiration, interception, surface runoff, and groundwater recharge for the 2008-2017 period. Output raster maps of calculated monthly values were summed in ArcGIS software (total 12 maps for each water balance component) to obtain the spatial distribution of average annual values ( Figure 5).
The actual evapotranspiration (AET) in Wetspass-M model is calculated as a sum of evaporation from bare soil, open water and impervious surface area, as well as transpiration and interception of vegetated area [22,32]. The simulated average monthly AET ranges from 9 to 80 mm/month, with average monthly interception between 0 and 21 mm/month. The average annual AET (Figure 5a) ranges from 142 to 2591 mm/year, with an average value of 414 mm/year. About 80% of the average annual AET occurs during the rainy and warmer period, from May to November. The lowest AET values are in the built-up area, higher values are represented by agricultural land, and the highest AET values are attributed to evaporation from water bodies. The average annual interception (Figure 5b) varies between 0 and 296 mm/year, with an average value of 111 mm/year. The highest values of average annual interception are observed at the area covered by forest and shrub. The results of AET and interception confirm that evapotranspiration depends greatly on land use. The average annual AET accounts for 45% of the average annual precipitation, meaning that evapotranspiration presents the major process by which water leaves the system. The actual evapotranspiration (AET) in Wetspass-M model is calculated as a sum of evaporation from bare soil, open water and impervious surface area, as well as transpiration and interception of vegetated area [22,32]. The simulated average monthly AET ranges from 9 to 80 mm/month, with average monthly interception between 0 and 21 mm/month. The average annual AET (Figure 5a) ranges from 142 to 2591 mm/year, with an average value of 414 mm/year. About 80% of the average annual AET occurs during the rainy and warmer period, from May to November. The lowest AET values are in the built-up area, higher values are represented by agricultural land, and the highest AET values are attributed to evaporation from water bodies. The average annual interception (Figure 5b) varies between 0 and 296 mm/year, with an average value of 111 mm/year. The highest values of average annual interception are observed at the area covered by forest and shrub. The results of AET and interception confirm that evapotranspiration depends greatly on land use. The average annual AET accounts for 45% of the average annual precipitation, meaning that evapotranspiration presents the major process by which water leaves the system.
The WetSpass-M model estimates monthly surface runoff in relation to precipitation The WetSpass-M model estimates monthly surface runoff in relation to precipitation amount, precipitation intensity, interception, and soil infiltration capacity [20]. The simulated average monthly surface runoff varies between 4 and 38 mm/month. The average annual surface runoff (Figure 5c) ranges from 37 to 987 mm/year, with an average value of 186 mm/year. About 50% of the average annual surface runoff occurs during the colder period from November to February. This is the period with low interception from vegetated surface and possible freezing of the soil occurs, resulting in lower infiltration of precipitation to the groundwater and higher surface runoff. High surface runoff values are observed in the water bodies, in the built-up area, and areas covered with soil of low permeability. The estimated average annual surface runoff represents 21% of the average annual precipitation.
The WetSpass-M model calculates monthly groundwater recharge as a residual term of the water balance R = P − S − AET (mm/month), where R is groundwater recharge, P is precipitation, S is surface runoff, and AET is actual evapotranspiration. The simulated average monthly groundwater recharge ranges from 14 to 57 mm/month. The average annual groundwater recharge (Figure 5d) varies between 0 and 511 mm/year, with an average of 312 mm/year. The spatial distribution of groundwater recharge of the study area greatly depends on soil type and land use. Lower recharge rates are generally observed at the area with low permeable soil (clay, silty clay), while higher recharge rates are associated with more permeable soil (sand, loam). In addition, higher values are attributed to agricultural areas, and especially forests. The highest recharge rates belong to the areas covered by forest on sandy soil. The average annual groundwater recharge constitutes about 34% of the average annual precipitation, which corresponds with previously used values of effective infiltration in the study area [10,11].

Boundary Conditions
The maps of equipotential lines show the behavior of the aquifer system ( Figure 6). The regional direction of groundwater flow in both hydrological conditions is from NW to SE, with local changes. The equipotential lines show that aquifer has an inflow boundary from Drava River and accumulation lake Varaždin on the northwestern and northern edge, no flow boundary on the western and southern edge, and an outflow boundary on the eastern edge. The lake water level variations are generally within 1 m, between 189 and 190 m a.s.l.
of 186 mm/year. About 50% of the average annual surface runoff occurs during the colder period from November to February. This is the period with low interception from vegetated surface and possible freezing of the soil occurs, resulting in lower infiltration of precipitation to the groundwater and higher surface runoff. High surface runoff values are observed in the water bodies, in the built-up area, and areas covered with soil of low permeability. The estimated average annual surface runoff represents 21% of the average annual precipitation.
The WetSpass-M model calculates monthly groundwater recharge as a residual term of the water balance where R is groundwater recharge, P is precipitation, S is surface runoff, and AET is actual evapotranspiration. The simulated average monthly groundwater recharge ranges from 14 to 57 mm/month. The average annual groundwater recharge (Figure 5d) varies between 0 and 511 mm/year, with an average of 312 mm/year. The spatial distribution of groundwater recharge of the study area greatly depends on soil type and land use. Lower recharge rates are generally observed at the area with low permeable soil (clay, silty clay), while higher recharge rates are associated with more permeable soil (sand, loam). In addition, higher values are attributed to agricultural areas, and especially forests. The highest recharge rates belong to the areas covered by forest on sandy soil. The average annual groundwater recharge constitutes about 34% of the average annual precipitation, which corresponds with previously used values of effective infiltration in the study area [10,11].

Boundary Conditions
The maps of equipotential lines show the behavior of the aquifer system ( Figure 6). The regional direction of groundwater flow in both hydrological conditions is from NW to SE, with local changes. The equipotential lines show that aquifer has an inflow boundary from Drava River and accumulation lake Varaždin on the northwestern and northern edge, no flow boundary on the western and southern edge, and an outflow boundary on the eastern edge. The lake water level variations are generally within 1 m, between 189 and 190 m a.s.l. Nearly all surface water features are in interaction with groundwater, except intake channel of HPP Varaždin, which is lined with concrete and has no impact on groundwater flow net. There is a clear bending of the equipotential lines towards derivation channel of Nearly all surface water features are in interaction with groundwater, except intake channel of HPP Varaždin, which is lined with concrete and has no impact on groundwater flow net. There is a clear bending of the equipotential lines towards derivation channel of HPP Varaždin, which suggest its drainage role. The impacts of the pumping site Vinokovšćak and Plitvica stream on groundwater flow net are not noticeable. The abstraction rate of groundwater at the well site Vinokovšćak (7315 m 3 /day on 3 September 2012 and 9350 m 3 /day on 15 September 2014) is clearly not enough to cause the groundwater level to drop significantly. Because there is not a sufficiently dense network of observation wells along the Plitvica stream (Figure 1), it is not possible to make a detailed interpolation along the stream to determine its contribution to the groundwater flow. Marković et al. [14] used stable water isotope measurements in Plitvica stream and adjacent observation wells, and concluded that Plitvica stream drains the aquifer. Comparing the maps of equipotential lines for low water levels ( Figure 6a) and high water levels (Figure 6b), it is evident that there is no significant change in the groundwater flow net. Oscillation in groundwater levels for a 10-year period are generally within 1-2 meters, which suggest that groundwater levels are strongly affected by the accumulation lake Varaždin and Drava River, which keep the aquifer in the quasi-steady state.

Conclusions
A combination of geological maps, borehole logs, cross-sections, DEM, land use, soil type, groundwater and surface water levels, and meteorological data was used to develop a hydrogeological conceptual model of the Varaždin alluvial aquifer. The hydrogeological conceptual model will be used for setting up a numerical groundwater flow and nitrate transport model.
The aquifer geometry is presented through a 3D model consisting of three layers: upper aquifer, semipermeable interlayer, and lower aquifer.
The groundwater recharge from precipitation was determined using WetSpass-M model for the period 2008-2017. The average annual actual evapotranspiration varies between 142 and 2591 mm/year, with about 80% of it occurring during the rainy period, from May to November. Lower values are observed in the built-up area, while higher values are attributed to agriculture and evaporation from water bodies. The average annual actual evapotranspiration represents 45% (414 mm/year) of the average annual precipitation. Estimated average annual surface runoff ranges from 37 to 987 mm/year, with an average value of 186 mm/year, which constitutes 21% of the average annual precipitation. About 50% of the average annual surface runoff occurs during the colder period from November to February. About 34% (312 mm/year) of the average annual precipitation recharges the aquifer, with 0 and 511 mm/year as a minimum and maximum average values, respectively. According to the results, permeability of the soil and land use control the spatial distribution of groundwater recharge. Higher values are associated with permeable soil types and agriculture or forest as a land cover. The results of groundwater recharge are consistent with previous researches, but with more detailed spatial distribution, which will serve as an input for a future numerical model of the Varaždin alluvial aquifer.
The general direction of groundwater flow is from NW to SE. The aquifer has an inflow boundary from Drava River and accumulation lake Varaždin on the northwestern and northern edge, and outflow boundary on the eastern edge. Western and southern edge of the aquifer are considered as no flow boundary. The equipotential lines show that the derivation channel of HPP Varaždin drains the aquifer, while the pumping site Vinokovšćak and Plitvica stream do not have a visible impact on groundwater flow net.
Author Contributions: I.K. conceptualization, methodology, software, writing-original draft preparation, review and editing; T.M. conceptualization, funding acquisition, supervision, writing-review and editing; T.V. data collection, visualization of meteorological data, writing-review and editing; O.L. data collection, writing-review and editing. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.