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

#### **Tamara Markovi´c 1,\*, Igor Karlovi´c 1, Melita Perˇcec Tadi´c <sup>2</sup> and Ozren Larva <sup>1</sup>**


Received: 9 December 2019; Accepted: 28 January 2020; Published: 30 January 2020

**Abstract:** 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 δ18O and δ2H in ground- and surface waters and precipitation samples were carried out. The relationship between ratios δ18O, δ2H, 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.

**Keywords:** stable isotopes; deuterium and oxygen-18; hydrogeological conceptual model; alluvial aquifer; Varaždin area

#### **1. 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 (δ18O) and hydrogen (δ2H) in atmospheric precipitation and groundwater is important for hydrological, hydrogeological, climatological, and meteorological applications [1–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 δ18O and δ2H in the ground- and 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 δ18O and δ2H values in precipitation and in the ground and surface waters, and this information will be used to determine recharge areas.

#### **2. Geological, Hydrogeological, and Climatic Setting of the Study Area**

#### *2.1. 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<sup>2</sup> (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š´cak (Figure 2). High concentrations of nitrate have caused the shutting down of the pumping site Varaždin, but Vinokovš´cak is still in operation.

**Figure 1.** Geographical position of the study area with locations of the sampling points; the transects A–A' and B–B' correspond to the representative hydrogeological cross-sections shown in Figure 2. The general groundwater flow direction is defined by the water heads and stable isotopes in the study area.

**Figure 2.** Schematic hydrogeological cross-sections across the study area (modified according to Reference [10]).

#### *2.2. 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.

**Table 1.** Monthly and annual average air temperature in Varaždin. Statistics in rows refer to maximum, mean, and minimum average monthly and annual temperature in the 1981–2010 period, and the standard deviation (sd).


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).


**Table 2.** The monthly and annual precipitation sum in Varaždin. Statistics in rows refer to maximum, mean, and minimum monthly and annual precipitation in the 1981–2010 period, the standard deviation (sd), and the coefficient of variation (cv).

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 p-value 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.


#### **3. Materials and Methods**

#### *3.1. Water Sampling*

Monthly composite precipitation was sampled in the Miko family courtyard in the village Hraš´cica (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 intoa1L 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š´cak 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.


**Table 5.** Depths of the observation wells.

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.

#### *3.2. Stable Isotope Analyses*

The δ18O and δ2H were determined using Picarro L2130i (Santa, Clara, USA) in the Hydrochemical Laboratory of the Croatian Geological Survey. The instrument uses CRDS (Cavity Ring-Down Spectroscopy) technology [14]. All measurements were checked with Picarro's standards (Depleted <sup>−</sup>29.6 <sup>±</sup> 0.2 <sup>δ</sup>18O; <sup>−</sup><sup>235</sup> <sup>±</sup> 1.8 <sup>δ</sup>2H; Mid <sup>−</sup>20.6 <sup>±</sup> 0.2 <sup>δ</sup>18O; <sup>−</sup><sup>159</sup> <sup>±</sup> 1.3 <sup>δ</sup>2H; Zero 0.3 <sup>±</sup> 0.2 <sup>δ</sup>18O; 1.8 <sup>±</sup> 0.9 <sup>δ</sup>2H), which were checked periodically against the International Atomic Energy Agency (IAEA) standards: Vienna Standard Mean Ocean Water 2 (VSMOW2) and Standard Light Antarctic Precipitation 2 (SLAP2). Measurement precision was <sup>±</sup> 0.3 <sup>o</sup>/oo for <sup>δ</sup>18O and <sup>±</sup> <sup>1</sup> <sup>o</sup>/oo for <sup>δ</sup>2H.

It is generally known that all isotopic results are expressed as per the international measurement standard, VSMOW2 [15,16].

A global relationship between δ2H and δ18O has been observed [16] and called the Global Meteoric Water Line (δ2H <sup>=</sup> <sup>8</sup>·δ18O <sup>+</sup> 10).

The deuterium excess (d-excess [17]) was calculated for each sample as follows:

$$\text{d-excess (\text{\textquotedblleft}\text{\textquotedblright})} = \text{\textquotedblleft}\text{H} - 8\text{ \textquotedblleft}\text{\textquotedblright}\text{O.}$$

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].

#### **4. Results and Discussion**

#### *4.1. 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).

**Figure 3.** Monthly anomaly of precipitation (**a**) and air temperature (**b**) in the period from June 2017 to June 2019 based on 1981–2010 climate normal.

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.

#### *4.2. Stable Isotopes in Precipitation*

The mean stable isotope δ18O and δ2H values and the associated d-excess are shown in Table 6.

**Table 6.** Yearly averaged isotopic composition of precipitation at the rain gauge in Hrašˇcica (177 m a.s.l.).


The stable isotope <sup>δ</sup>18O values varied from <sup>−</sup>14.91 to <sup>−</sup>4.5 <sup>o</sup>/oo, and <sup>δ</sup>2H varied from <sup>−</sup>108.1 to <sup>−</sup>25.1 <sup>o</sup>/oo (Figure 4, Table S1). The lowest <sup>δ</sup>18O and <sup>δ</sup>2H values were observed in winter and highest were observed in summer. The d-excess varied from 4.1 to 12.9 <sup>o</sup>/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 δ18O and δ2H values were automatically more negative. In addition, the lowest δ18O and δ2H 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.

**Figure 4.** Monthly variations of (**a**) δ18O (o/oo) values, (**b**) δ2H (o/oo) values, and **(c)** d-excess (o/oo) in precipitation at the rain gauge in Hraš´cica.

The measured stable isotope δ18O and δ2H 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 δ18O and δ2H values and weighted by the amount of precipitation. Because of this, they are not discussed here.

The calculated LMWL for the period from June 2017 to June 2019 is:

$$\text{OL.SR } \delta^2 \text{H} = (7.54 \pm 0.12) \text{ } \delta^{18} \text{O} + (5.00 \pm 1.00) \text{, n} = 23$$

$$\text{RMA } \delta^2 \text{H} = (7.56 \pm 0.11) \text{ } \delta^{18} \text{O} + (5.17 \pm 1.04) \text{, n} = 23$$

$$\text{PWL.SR } \delta^2 \text{H} = (7.55 \pm 0.13) \text{ } \delta^{18} \text{O} + (4.85 \pm 1.13) \text{, n} = 23.$$

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 <sup>o</sup>/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 <sup>o</sup>/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).

**Figure 5.** Distribution of groundwater δ18O (o/oo) and δ2H (o/oo) values around the local meteoric water line (LMWL).

#### *4.3. 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 7.** Minimum, maximum, and averaged isotopic composition of the groundwater by sampled well.

**Table 8.** Minimum, maximum, and averaged isotopic composition of surface water.


The measured <sup>δ</sup>18O values in the groundwater varied from <sup>−</sup>11.47 to <sup>−</sup>8.26 <sup>o</sup>/oo, and the <sup>δ</sup>2H values varied from <sup>−</sup>81.7 to <sup>−</sup>58.5 <sup>o</sup>/oo (Table 7). The measured <sup>δ</sup>18O values in the surface water varied from <sup>−</sup>12.12 to <sup>−</sup>3.36 <sup>o</sup>/oo, and the <sup>δ</sup>2H values varied from <sup>−</sup>81.7 to <sup>−</sup>38.8 <sup>o</sup>/oo (Table 8).

The correlation between the δ18O and δ2H 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 δ18O and δ2H 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 δ18O and δ2H 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 <sup>o</sup>/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 <sup>o</sup>/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 δ18O and δ2H values of the surface waters distributed around the LMWL shown in Figure 6 indicate a relationship between δ18O and δ2H 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 δ2Hand δ18O 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.

**Figure 6.** Distribution of surface waters δ18O (o/oo) and δ2H (o/oo) values around the LMWL.

**Figure 7.** Distribution of δ18O (o/oo) values in surface waters over the monitored time.

To connect the measured values, a simplified statistical correlation method was used, and results are shown in the correlation matrix in Table 9.


**9.**Correlationmatrixbetweengroundandsurfacewatersforδ18O(o/oo)values.

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š´cak pumping site (Figure 1), and the influence of the local precipitation is minor.

#### **5. Conclusions**

Although δ18O and δ2H values of ground- and surface waters and precipitation were measured for only 24 months, the following conclusions are proposed:


**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/2/379/s1, Table S1: Determined isotopic composition of the precipitation by month for monitored period.

**Author Contributions:** T.M. contributed to conceptualization, isotopic data interpretation, and writing. I.K. contributed to the hydrogeology and geology description of the study area and graphical editing. M.P.T. contributed to statistical analysis and visualization of climatological data. O.L. contributed by polishing the article. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Croatian Scientific Foundation (HRZZ), grant number HRZZ-IP-2016-06-5356.

**Acknowledgments:** The authors would like to thank the Miko family for taking care of the rain gauge, as well as the Croatian Meteorological and Hydrological Service for providing climatological data for the Varaždin meteorological station.

**Conflicts of Interest:** The authors declare that there is no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
