Application of Stable Water Isotopes to Improve Conceptual Model of Alluvial Aquifer in the Varaždin Area

: To understand groundwater flow and geochemical processes within an aquifer, it is necessary to set up a conceptual model of the aquifer. To accomplish this, different methods are used, and one of them is an isotopic technique. The study area is located in the Varaždin area (NW Croatia). The aquifer represents the main source of potable water for the town of Varaždin and the surrounding settlements. The conceptual model of the alluvial aquifer has to be set up prior to creating a groundwater flow and transport model. Measurements of ratios δ 18 Ο and δ 2 H in ground ‐ and surface waters and precipitation samples were carried out. The relationship between ratios δ 18 Ο , δ 2 H, and d ‐ excess for local precipitation in the study area showed that precipitation originates from the Atlantic air masses, although during the colder periods of the year, influence of the Mediterranean air masses was not negligible. The monitored period was warmer and wetter than average. Evaporation was observed at all monitored surface waters, but the largest rate was at the location of a gravel pit in Šijanec. The isotopic composition of the precipitation and groundwater showed a good correlation due to the isotopic homogenization of groundwater along the flow path.


Introduction
To understand groundwater flow and geochemical processes within an aquifer, it is necessary to set up a conceptual model of it. Different methods are used, and one of them is an isotopic technique, which is often used successfully to help elucidate hydrological studies [1]. Knowledge about the isotopic ratios of oxygen (δ 18 Ο) and hydrogen (δ 2 H) in atmospheric precipitation and groundwater is important for hydrological, hydrogeological, climatological, and meteorological applications [1][2][3][4][5][6][7][8] because it can provide information on the mean recharge elevation of the aquifer, the mean residence time, water-rock interactions, etc.
The study area is located in the Varaždin area (NW Croatia). The aquifer represents the main source of potable water for the town of Varaždin and the surrounding settlements. The favorable climate, topography, and available groundwater have insured intensive agricultural practices involving the application of large amounts of synthetic fertilizers and manure that have subsequently led to high nitrate concentrations in the Varaždin aquifer. High concentrations of nitrate have caused the shutting down of the Varaždin pumping site. To determine the behavior of nitrates, a groundwater flow and transport model will be used. A conceptual model of the alluvial aquifer, therefore, has to be set up, and to improve this model, measurements of δ 18 Ο and δ 2 H in the groundand surface waters and precipitation samples were carried out.
The research is still ongoing and, in this paper, only the results of two-year measurements are elaborated in detail. The main goal of this paper is to provide an overview of the spatial and temporal variability in δ 18 Ο and δ 2 H values in precipitation and in the ground and surface waters, and this information will be used to determine recharge areas.

Geological and Hydrogeological Setting of the Study Area
The study area of the Varaždin aquifer system is located in NW Croatia, upstream of the town of Varaždin in the valley of the Drava River and it covers an area of approximately 200 km 2 ( Figure  1). The Varaždin alluvial aquifer is composed of sediments of the Quaternary age, deposited during the Pleistocene and Holocene eras as a result of accumulation processes of the Drava River [9]. The alluvial deposits consist primarily of gravel and sand with occasional lenses and interbeds of silt and clay ( Figure 2). The thickness of the aquifer varies from less than 5 m in the NW part to about 65 m in the SE part of the study area ( Figure 2, cross-section A-A'). The Varaždin aquifer is an unconfined aquifer, and the groundwater is in direct contact with the surface water: the Drava River and the Plitvica stream, the derivation channel, the accumulation lake Varaždin (Varaždin Lake), and the gravel pit Šijanec ( Figure 2). The direction of groundwater flow is generally NW-SE and is parallel to the direction of the Drava River flow. The covering layer of the aquifer is not continuously developed ( Figure 2), which makes the aquifer vulnerable to contamination from the surface. Favorable hydrogeological conditions enabled the development of two pumping sites-Varaždin and Vinokovšćak (Figure 2). High concentrations of nitrate have caused the shutting down of the pumping site Varaždin, but Vinokovšćak is still in operation.

Climatic Setting of the Study Area
The Varaždin meteorological station (46.3° N, 16.137° E, 167 m a.s.l.) and the surroundings have a characteristic precipitation regime with more precipitation during the summer [11]. Because of these features, the local climate is categorized in the Cfb group of the Köppen-Geiger classification system, which is known as "warm-temperate climate" or "marine west coast climate." The study area is the former [12].
The coldest month is January, with average temperatures of 0.0 °C, while July and August are the warmest months with average temperatures of 20.9 and 20.1 °C, respectively (Table 1). Mean annual temperature is 10.6 °C. The inter-annual variability is the largest in February and January (standard deviation (sd) in Table 1), meaning that average differences from the mean are largest at the end of the winter. According to the data from the last climate normal period (1981-2010), the lowest average precipitation amounts were during the cold part of the year, and in January the mean precipitation was 38.7 mm ( Table 2). From June to September, the precipitation amounts were on average larger than 80 mm, with the highest amount in September (98.3 mm) and the second-highest in June (96.1 mm). The average annual precipitation was 832 mm. Inter-annual variability was largest in January (coefficient of variation (cv) = 0.81). Based on the long-term data records for 1951-2018, the existence of seasonal trends was tested by a Mann-Kendall test at the 0.05 significance level, and Sen's slope was calculated to determine the trend value [13]. There is significant warming in all seasons, with Sen's slope ranging from 1.8 °C/100 years in autumn to 3.7 °C/100 years in summer (Table 3). Table 3. Trend analysis of long-term temperature data for the period 1951-2019. A Mann-Kendall pvalue of <0.05 indicates a significant trend. Sen's slope is a value of that trend. A Sen's slope p-value of <0.5 indicates a significant value of a trend. There are no statistically significant changes in seasonal monthly precipitation in the 1951-2019 period (Table 4), even though there is a slight indication of drying in summer and in spring, while autumn and winter show a slight tendency of becoming wetter. Table 4. Trend analysis of long-term precipitation data for the period 1951-2019. Parameters are the same as in Table 3.

Water Sampling
Monthly composite precipitation was sampled in the Miko family courtyard in the village Hrašćica (46.3° N, 16.292° E; 177 m a.s.l.) in the period from June 2017 until June 2019. In the field, sample was poured into a 1 L plastic bottle with a tight-fitting cap.
In addition, ground-and surface water sampling was conducted on a monthly basis in the period from June 2017 to June 2019 for chemical and isotopic analyses. Ten observation wells (nine of them are located in the recharge area of the Varaždin pumping site and one in the recharge area of the Vinokovšćak pumping site) and four surface waters (Drava River, the Plitvica stream, Varaždin Lake, which is an accumulation, and a gravel pit in Šijanec) were sampled. The water depth was measured prior to pumping. The observation wells were sampled after stabilization of the parameters EC (electrical conductivity), T (water temperature), pH, and O2 (dissolved oxygen content). Depths of the observation wells are given in Table 5. Samples were poured into a 50 mL plastic bottle with a tight-fitting cap. All samples were measured in the laboratory immediately upon returning from the field. Meteorological parameters (precipitation and air temperature) used here are from the main Varaždin meteorological station, located in the vicinity of the installed rain gauge.
It is generally known that all isotopic results are expressed as per the international measurement standard, VSMOW2 [15,16].
The deuterium excess (d-excess [17]) was calculated for each sample as follows: In 2003, researchers verified this value by using global maps derived by interpolation from more than 340 stations [18] and, in 2005, this was reconfirmed using IAEA-GNIP datasets at 410 stations [19]. Higher values of d-excess are caused by intense evaporation of seawater in conditions of moisture deficit [20].
The local meteoric water line (LMWL) has been calculated using three types of linear regression analysis, two of which are recommended by the IAEA [21]: ordinary least squares regression (OLSR) and reduced major axis (RMA) regression. The new one takes into account the amount of precipitation, using precipitation weighted least squares regression (PWLSR) [22].

Precipitation and Temperature
The climatological conditions from June 2017 to June 2019, the period of the isotope sampling, were compared to a climate normal 1981-2010 ( Figure 3).
The summer of 2017 was drier and warmer from the average of 1981-2010. Autumn 2017 was mostly wetter than normal, with September 2017 being the wettest month of the period, and it was also colder than average in September. This period, with larger precipitation than average, continued until May 2018, and it was also warmer, except in February and March 2018. From June 2018 to March 2019, precipitation was mostly around or smaller than average, and it was warmer. May 2019 was the second wettest month during this period, and it was also colder than average. June 2019 was again warmer than average.

Stable Isotopes in Precipitation
The mean stable isotope δ 18 Ο and δ 2 H values and the associated d-excess are shown in Table 6.   Table S1). The lowest δ 18 Ο and δ 2 H values were observed in winter and highest were observed in summer. The d-excess varied from 4.1 to 12.9 o /oo. The d-excess shows the influence of the Atlantic air masses. Nevertheless, the influence of the Mediterranean air masses in the study area was observed during the autumn and winter months ( Figure 4). The Mediterranean air masses (precipitation) are characterized by a higher d-excess than the Atlantic air masses [20]. This was observed in References [23,24] in the continental part. Atypical climatological conditions during the observed period had influenced variations of monthly isotopic composition. A sudden change in the air temperature and/or precipitation amount during the season influenced the variation of the monthly isotopic composition of the rain. For example, May 2019 was colder and wetter than average (even than May 2018), and δ 18 Ο and δ 2 H values were automatically more negative. In addition, the lowest δ 18 Ο and δ 2 H values were measured in the coldest month, which was February 2018 (Figures  3 and 4). It was observed that the isotopic composition of the precipitation in the study area reflects climatological conditions well.
The measured stable isotope δ 18 Ο and δ 2 H values were weighted by the amount of precipitation at the Varaždin meteorological station for the observed period. However, there were no large differences between the measured stable isotope δ 18 Ο and δ 2 H values and weighted by the amount of precipitation. Because of this, they are not discussed here.
It was observed that all three methods yielded a very similar slope value and axis intercept (b value) of the LMWL, which was supported by very similar measured and weighted values.  In addition, calculated data were compared with data published in Reference [25]. There is a difference between these two slopes values for 0.09 o /oo, and it can be concluded that the OLSR values calculated from the measured data and the published meteoric water line of the study area are not different in terms of slope. However, there is a large difference between these two lines in the axis intercept values for 2.6 o /oo, and the published LMWL is slightly below the measured one ( Figure 5). There are several reasons for that: shorter monitored period in our research; very untypical and extreme climatological conditions during our monitored period; and different measurement techniques (our samples were measured using CRDS technology and published were measured using Isotope Ratio Mass Spectrometry (IRMS) technology).

Stable Isotopes in Ground and Surface Waters
The minimum, maximum, and average isotopic composition of the surface-and groundwaters are given in Tables 7 and 8 together with their average d-excess.   Table 8).
The correlation between the δ 18 Ο and δ 2 H measured values of the groundwater is shown in Figure 5 and indicates that this relationship has a slope of 7.14. Using a Student's t-test according to Reference [26], a good relationship between groundwater and precipitation was observed. Generally, an isotope relationship between δ 18 Ο and δ 2 H with a slope of about 8 is normally observed for precipitation [16]. Since the relationship between the isotopic composition of precipitation and groundwater is good, it can be concluded that groundwater is recharged by precipitation. Values that are slightly more negative were measured in the SPV-11 well and the private well, while at observation wells PDS-5, PDS-6, PDS-7, and P-1529, values are almost identical ( Table 7). The highest δ 18 Ο and δ 2 H values were measured at observation wells P-4039 and P-2500 ( Table 7). The calculated average d-excess values varied from 9.6 to 10.7 o /oo, indicating the influence of recharge by precipitation with signatures of the Atlantic air masses and good homogenization of groundwater along the flow path. This was observed at SPV-11, PDS-6, PDS-7, the private well, P-1530, and P-1529. However, depending on hydrodynamic conditions (low/high water levels), the vicinity of the river or lake, and the depth of the observation well, it was observed that wells, especially the shallower ones and/or those closer to the river and lake, showed high variation in d-excess values. These values were higher than 11 o /oo, indicating recharge by surface waters and faster recharge by precipitation. This was observed at P-1556, PDS-5, P-2500, and P-4039.
The measured δ 18 Ο and δ 2 H values of the surface waters distributed around the LMWL shown in Figure 6 indicate a relationship between δ 18 Ο and δ 2 H for surface waters with slopes of 5.73 at Varaždin Lake, 6.32 at Drava River, and 4.77 at the gravel pit in Šijanec, indicating an influence of evaporation. A slope from 4 to 6 is attributed to waters with a significant rate of evaporation relative to the input [16]. It was observed that the evaporation process was strongest at the location gravel pit in Šijanec (Figure 7). Nevertheless, the gravel pit was used for fish farming. Because of this activity (resulting in an extra nutrient load due to fish feeding), a low water level, a high load of nutrients, high temperatures, and algae bloom occurred every summer, which had a significant influence on the isotopic and chemical features of this water. The winter-measured values of the δ 2 Hand δ 18 O of the Drava River and the summer-measured values of Varaždin Lake are above the LMWL (Figure 6). For the Drava River, this can be explained by the fact that the larger part of the recharge area of the Drava River is situated far upstream of the study area and is under the influence of a different climate, and the influence of the recharge area in the study area is small. Varaždin Lake is recharged by the Drava River, especially in the late winter and springtime when isotopic values are more negative in the river. Since the lake has a high volume, turnover in the lake takes some time. In addition, the Plitvica stream, like the Drava River, has its recharge area in a mountain area where the climate is different, and because of that, more negative values are measured in the wintertime. During the late spring/summer period, the discharge of the stream is low. In the watercourse of the stream, small connected ponds are formed where evapotranspiration is present and, because of that, some values are below the LMWL.  To connect the measured values, a simplified statistical correlation method was used, and results are shown in the correlation matrix in Table 9. No statistical connection was observed between either the groundwater, the waters of Drava River, Varaždin Lake, or the Plitvica stream with water from the gravel pit in Šijanec. This is partly because there was no observation immediately downstream from the gravel pit. Moreover, it has a very small volume, and its influence is thus limited to its immediate surroundings. Furthermore, the gravel pit is not connected to the river, stream, or lake; consequently, the isotopic composition differs (Figure 1). In addition, a very weak correlation was observed between the waters from the P-2500 and P-4039 wells, which are in the vicinity of the Plitvica stream. This weak correlation is attributed to the drainage roll of the Plitvica stream in this part of the aquifer. A higher correlation was observed between both Varaždin Lake and the Drava River waters and the groundwater of the observation wells. A high correlation between observation wells on the right side of the intake/drain channels indicates a homogenization of the groundwater source (a mixture of the precipitation, river, and lake waters). The left side is different because of the influence of Varaždin Lake. Namely, the Drava River flows from the lake and has the same isotopic composition as the lake water ( Figure 1). The river recharges the aquifer as a consequence of groundwater abstraction at the Vinokovšćak pumping site (Figure 1), and the influence of the local precipitation is minor.

Conclusions
Although δ 18 Ο and δ 2 H values of ground-and surface waters and precipitation were measured for only 24 months, the following conclusions are proposed: (a) Meteorological conditions during the observed period were quite different from the average. The summer of 2017 was drier and warmer, while the autumn was wet and cold. This period had more precipitation than average, and this continued until May 2018. It was also warmer, except in February and March 2018. From June 2018 to March 2019, precipitation was mostly around or smaller than average, and it was warmer. May 2019 was the second wettest month during this period and was also colder than average. June 2019 was again warmer than average. (b) The isotopic composition of the local precipitation in the study area varied according to the variable meteorological conditions and relationship between δ 18 Ο and δ 2 H, and the d-excess showed that precipitation originated from the Atlantic air masses, although during the colder parts of the year, the influence of the Mediterranean air masses was not negligible. (c) Evaporation was observed at all surface waters, but the largest rate was at the location of the gravel pit in Šijanec. Moreover, it has a very small volume; therefore, its influence is limited to the immediate surroundings of the aquifer. (d) The isotopic composition of precipitation and groundwater showed a good correlation due to the isotopic homogenization of groundwater along the flow path. (e) The isotopic composition of the groundwater source on the right side of the intake/drain channels indicates a homogenization of the groundwater source (a mixture of precipitation, river, and lake waters), whereas the left side is different because of the higher influence of Varaždin Lake.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Table S1: Determined isotopic composition of the precipitation by month for monitored period.