**Application of Stable Isotopes of Water to Study Coupled Submarine Groundwater Discharge and Nutrient Delivery**

#### **Carlos Duque 1,2,\*, Søren Jessen 2, Joel Tirado-Conde 2, Sachin Karan <sup>3</sup> and Peter Engesgaard <sup>2</sup>**


Received: 12 August 2019; Accepted: 31 August 2019; Published: 4 September 2019

**Abstract:** Submarine groundwater discharge (SGD)—including terrestrial freshwater, density-driven flow at the saltwater–freshwater interface, and benthic exchange—can deliver nutrients to coastal areas, generating a negative effect in the quality of marine water bodies. It is recognized that water stable isotopes (18O and 2H) can be helpful tracers to identify different flow paths and origins of water. Here, we show that they can be also applied when assessing sources of nutrients to coastal areas. A field site near a lagoon (Ringkøbing Fjord, Denmark) has been monitored at a metric scale to test if stable isotopes of water can be used to achieve a better understanding of the hydrochemical processes taking place in coastal aquifers, where there is a transition from freshwater to saltwater. Results show that 18O and 2H differentiate the coastal aquifer into three zones: Freshwater, shallow, and deep saline zones, which corresponded well with zones having distinct concentrations of inorganic phosphorous. The explanation is associated with three mechanisms: (1) Differences in sediment composition, (2) chemical reactions triggered by mixing of different type of fluxes, and (3) biochemical and diffusive processes in the lagoon bed. The different behaviors of nutrients in Ringkøbing Fjord need to be considered in water quality management. PO4 underneath the lagoon exceeds the groundwater concentration inland, thus demonstrating an intra-lagoon origin, while NO3, higher inland due to anthropogenic activity, is denitrified in the study area before reaching the lagoon.

**Keywords:** submarine groundwater discharge; stable isotopes of water; nutrients; freshwater– saltwater interface; Ringkøbing Fjord

#### **1. Introduction**

Coastal regions are exposed to anthropogenic pressures frequently associated with high-density populations that require water supply [1]. Additionally, these regions are recipients for nutrients derived from farming, industrial activity, and sewage systems within the coastal region as well as their inland catchments [2]. Excess nutrient loads cause eutrophication [3–5] and algal blooms in coastal waters [6,7], and ultimately may have a global scale impact by altering marine food webs [2,8]. Another consequence is the effect on touristic activities such as fishing and leisure, [9] which are important drivers of economic activity.

The impact of groundwater-borne nutrients to the coastal regions was recognized by pioneering studies in the 1980s [10–12]. Subsequent investigations on submarine groundwater discharge have used seepage meters, radioactive tracers, or hydraulic methods to estimate the flux of groundwater [13–15]. Flux estimates combined with nutrient concentrations were used to calculate nutrient loads to coastal waters in several regions of the world [16–20].

Nevertheless, coastal aquifers present a series of hydrodynamic features that complicate the direct application of flux and nutrient concentration to estimate nutrient loads. The presence of a freshwater–saltwater interface separates two domains where flow is driven by different mechanisms, and where the hydrochemical properties of groundwater are distinct. Freshwater flows as a consequence of the higher hydraulic head inland after the recharge of aquifers. Saltwater fluxes are associated with density-driven flow [1,21]. The benthic fluxes, referred to here as the exchange between groundwater and surface water across the sea or lagoon bed [22], can be generated by different mechanisms including, among others, tides and waves [23]. Also, the mixing of water with different properties and ionic strengths can trigger multiple chemical transformations [24,25], making an accurate estimation of mass loads difficult. While the general trends are well-known, the complexity of coastal aquifers (i.e., spatial variability in hydraulic conductivity, flow paths, residence times) hinders the understanding of all the processes that are taking place, as well as their hydrochemical consequences. Recently, it was recognized that there are limited studies where there is a distinction between the different types of fluxes [26]. One reason for this is that it is required to differentiate between fluxes; however, that can be challenging due to the amount of field work required or the lack of appropriate methods to discern these fluxes. Despite the solid theoretical and modeling knowledge about hydrochemical processes [24,25], detailed field studies are required to link them to conditions that characterize aquifers, often with a high degree of complexity.

Using stable isotopes of water as a tracer is a well-established discipline in hydrogeological studies that allows the differentiation of water exposed to different physical processes relating to precipitation and evaporation [27]. Its use has been applied with multiple objectives [28]: To quantify recharge sources [29–31], for climate assessment [32], or to evaluate water origin in lakes [33] or hydrothermal systems [34,35]. The continuous improvement in analytical techniques has enabled researchers to quantify 18O and 2H with small volumes of water through a simple and direct sampling system. Since chemical processes do not have an effect on stable isotopes, they can be suitable tracers when combined with hydrochemical studies. If different sources of water are exposed to different physical processes, as naturally occurring in coastal areas, water stable isotopes can be an indicator of flow paths and the mixing of water, without being affected by chemical reactions.

The utility of water stable isotopes for studying the impact of the hydrogeological dynamics in coastal areas on nutrient discharge has been tested in this study. To this end, a field campaign was designed to combine the spatial distribution of nutrients and stable isotopes of water, where groundwater transitions from fresh to saline. Hence, the objective was to evaluate the use of stable isotopes of water, in combination with the spatial distribution of nutrients (PO4 and NO3) for characterizing submarine groundwater discharge (SGD) at the metric scale, including the challenges associated with field sampling (i.e., spatial variability, scattered samples, and outliers). The field site was established in a Danish lagoon (Ringkøbing Fjord), as a benchmark in which there are also major local concerns about water quality and potential eutrophication issues.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Ringkøbing Fjord is a brackish coastal lagoon located adjacent to the North Sea near the west coast of Jutland in Denmark (Figure 1). The lagoon covers 300 km2 with a mean water depth of 1.9 m (Ringkøbing Amt., 2004). Lagoon salinity ranges from 5 to 12 g/L during the year, due to the seasonal discharge from the Skjern River and in connection with the sea managed by a sluice. The region is characterized by a flat topography and low altitude hills. Annual mean precipitation and evaporation are 1050 mm and 630 mm, respectively [36].

The local hydrogeological setting consists of a shallow 12-m thick unconfined sandy aquifer, limited by a clay layer that was 5–10 m thick [37]. Other aquifer units below this range are out of the scope of this study, since their input to the lagoon in the near-shore environment would not be

quantitatively relevant at the study area. Groundwater originates as recharge in agricultural and residential (summer housing) areas flowing from east to west, until discharging to the lagoon [37,38]. Nutrient transport could be generated by farming activity as NO3 [39], and is associated with septic tanks and PO4 [40].

Ringkøbing Fjord has been studied in the last decades, focusing on the status of the lagoon water and the in-fjord processes affecting water quality [41–43]. The management of the sluice changes the salinity in the lagoon, and the impact on the ecological systems was described by Petersen et al. [44]. However, the groundwater flow with the associated nutrient input to the lagoon has not been researched. Previous groundwater-related research focused exclusively on the physical processes that determine the flow of groundwater to the lagoon [22,37], and the development of methodologies to help quantify water fluxes [38,45].

**Figure 1.** Study area location in Denmark and in Ringkøbing Fjord. Sampling points (1–30), inland wells, and general characteristics of the study area.

#### *2.2. Field Campaigns*

Field campaigns were conducted between the 14–19 August 2016, (Inland wells A and B) and between 7–14 July 2017 (sampling points 1–30). The data collection aimed to characterize the hydrogeological properties, fluxes, and hydrochemistry, including: freshwater aquifer inland (Inland wells A and B), shallow part of the aquifer underneath the lagoon (sampling points 1–30), and groundwater–lagoon water fluxes. The area was investigated in previous studies [22,38], allowing the optimization of the placement of sampling points 1–30 and the requirements to complete field work. Sampling points 1–30 extended along a transect from the shoreline up to 40 m offshore, with small differences in the water depth of the lagoon (water depth of 30 to 50 cm). The water depth of the lagoon did not exceed 50 cm for hundreds of meters offshore. Part of the transect was surrounded by submerged autochthonous vegetation; however, the field equipment was deployed along a corridor without vegetation, which serves as access to the lagoon for leisure and fishing. The first few meters of the lagoon from the shoreline were covered with

organic material (degraded vegetation), followed by well-classified sand that extended to at least half a meter below the surface (tested by surficial core samples).

In 2016, groundwater samples were collected at two wells (Inland wells A and B) located inland 60 and 15 m from the lagoon shore, respectively (Figure 1). These wells were 8 and 9 m deep, respectively, and samples were collected at one-meter intervals.

In 2017, 30 measurement locations were established irregularly distributed along a transect (Figure 1). The distance between sampling points ranged between 1–3 m. This allowed for data to be obtained along the zone, with a transition from fresh to saline/brackish groundwater, gaining a representation of the spatial heterogeneity via multiple measurements at similar distances from the shore. At all locations, a set of measurements were completed both at the surface of the lagoon bed and at multiple depths to determine groundwater hydrochemistry (number of samples, n = 150), water isotopic characteristics (n = 150), aquifer hydraulic properties (n = 120), and groundwater seepage to the lagoon (n = 30). In four cases, it was not possible to pump up groundwater from piezometers due to local low hydraulic conductivity, or potential damage to the piezometer screens. Electrical conductivity (EC) data were presented as continuous in the study area, through the kriging interpolation method. The stage level in the lagoon and the water table in the Inland A and B wells were monitored regularly, showing steady conditions indicating that the results presented here correspond to a stable snapshot of the groundwater–lagoon water interaction. The sampling protocol and analytical methods (see below) were the same for both campaigns (August 2016 and July 2017).

#### *2.3. Measurement of Fluxes*

Seepage meters were deployed one week prior to the campaign in July 2017 to allow the flushing of the seepage meter chambers. Seepage meters had a classic Lee design [46], and were constructed from 200-L 58-cm diameter steel drums. To minimize the disturbance of the sediments in the study area, the fluxes of groundwater discharge were measured in the first two days of the campaign [47]. The seepage meter bags were thin-walled, durable plastic autoclave bags prefilled with approximately 1 L of freshwater, with known electrical conductivity (EC ~0.28 mS cm−1). The measurements were repeated with 2–3 h measuring intervals on day 1 (round 1: 10:00–13:00, round 2: 13:00–15:00, round 3: 15:00–17:00) and day 2 (round 1: 9:00–11:00, round 2: 11:00–14:00, round 3: 14:00–16:00). After each interval, the sample bag was removed, and another prepared bag was immediately attached, to prevent inflow of lagoon water to the seepage meter during the bag exchange. Each bag had a valve that was closed during removal and while unattached. To approximate the simultaneous sampling of the seepage meters, the bag exchange was completed by two to three operators. The total time for sample collection and bag replacement never exceeded 20 min. The EC of seepage meter samples, corrected to 25 ◦C, was measured using a HACH IntelliCAL CDC401 Standard Conductivity probe, connected to a HACH HQ40d instrument (expected error less than 0.5% of the reading). The EC of measured flux was estimated with a mass balance:

$$\text{EC flux} = \frac{\left[\left(\text{W}\_{\text{off}} - \text{W}\_{\text{bag}}\right) \times \left(\text{EC}\_{\text{off}}\right)\right] - \left[\left(\text{W}\_{\text{on}} - \text{W}\_{\text{bag}}\right) \times \left(\text{EC}\_{\text{on}}\right)\right]}{\text{W}\_{\text{off}} - \text{W}\_{\text{on}}} \tag{1}$$

where Woff: Weight of the seepage meter bag after measuring, ECoff: EC of water in the seepage meter bag after measuring, Won: Weight of the seepage meter bag before measuring, ECon: EC of water in the seepage meter bag before measuring, and Wbag: Weight of the seepage meter bag and attachment accessories without water.

The proportion of freshwater and saltwater flux from each discharge sample was calculated by assuming an EC of 0.28 mS cm−<sup>1</sup> for the freshwater endmember measured in Inland wells A and B, while the EC of the saline lagoon endmember was assumed to be equal to that measured for the lagoon water on the day of sampling (22 mS cm<sup>−</sup>1). The results presented are the arithmetic mean of the set of six tests.

#### *2.4. Wells*

Two types of drilling methods were used depending on the desired depth. For shallow depths (up to 20 cm), a mini-piezometer was used (screen length 3 cm), to allow the sampling of small volumes with minor disturbance to the system, as the proximity to the lagoon bed could induce the vertical downward flow of lagoon water to the piezometer screen. Additionally, this sample could be a better indicator of the chemical characteristics of water flowing to the lagoon than deeper samples, as sampling was very close to the discharge point. In the same locations, a borehole was drilled with a <sup>3</sup> <sup>4</sup> inch steel pipe (screen length 10 cm) with a pneumatic hammer to depths of 1 m, 2 m, 3 m, and 4 m below the lagoon bed (sampling points 1–30), and for Inland wells A and B, every meter after reaching the water table (3–8 m or 1–9 m below terrain, respectively). For each location, sample collection was followed by slug tests, prior to pushing the well to the next depth. This system was used to facilitate the completion of hydraulic conductivity tests, and because it was difficult to reach depths below 1 m with the mini-piezometers.

#### *2.5. Hydrochemistry Sampling and Chemical Analysis*

In July 2017, groundwater samples were collected from the piezometers immediately after installation and purging (3 times the water volume contained in the well or after drying it) using a peristaltic pump. Water samples were collected in 60-mL polypropylene syringes. In total, 146 samples were collected, of which 69 were immediately passed through a 0.21-μm filter (Sartorius Minisart cellulose acetate) with a syringe, and directly into the collection containers. Filtration of the remaining 77 samples, of which most were saline, proved difficult and these were collected unfiltered into 100-mL acid washed amber glass bottles, completely filled to avoid headspace. These were filtered in the laboratory using a Mikrolab Aarhus A/S ML 303810 pneumatic syringe filtration unit. Upon collection in the field, the samples were immediately refrigerated to 5 ◦C. In either case, the filtrate was passed into one 10-mL polyethylene (PE) vial for analysis of major cations (Na+, K+, Ca2+, Mg2+) and ammonium (NH4 <sup>+</sup>), another 10-mL PE vial for major anions (Cl−, SO4 <sup>2</sup>−, NO3 <sup>−</sup>, Br−), a 5-mL PE vial for PO4 3− and Fe2<sup>+</sup> analysis, a 100-mL amber glass flask for alkalinity titration, and a 2-mL septum glass for stable isotopes of water. In addition to filtration and refrigeration, samples for PO4 <sup>3</sup><sup>−</sup> and Fe2<sup>+</sup> were conserved by 2 vol%2MH2SO4, samples for major cations and ammonium were acidified by 1 vol% 1 M HNO3 (resulting in pH 2–2.1, depending on alkalinity), and samples for anions were stored frozen (−18 ◦C). In the August 2016 campaign, the procedure was similar; however, here a flow cell was used measuring O2, pH, and EC directly in the field. A sample was taken when the readings were stable; then, they were filtered as described above.

Cations, anions, and ammonium were analyzed by ion chromatography using a Metrohm 820-IC Separation Center with 819-IC Detector. Alkalinity was measured by end point titration on a Metrohm autotitrator.

PO4 <sup>3</sup><sup>−</sup> were analyzed spectrophotometrically using the molybdate blue method of Murphy and Riley (1986). Total dissolved Fe, henceforth assumed to a direct proxy for the Fe2<sup>+</sup> concentration, was measured on residual sample aliquots of the PO4 <sup>3</sup><sup>−</sup> samples after reduction to Fe2<sup>+</sup> with hydroxylamine hydrochloride, using the spectrophotometric ferrozine method of Stookey (1970).

δ18O and δD values, reported in % units relative to the VSMOW/SLAP scale, were measured at the Geological Survey of Denmark and Greenland (GEUS) on a Picarro Cavity Ring-Down Spectrometer (CRDS) L2120-i. The standard deviation for each analysis was equal to or lower than 0.31% for δ18O and 0.74% for δD.

Field measurements of pH and temperature-specific electrical conductivity (EC; corrected to 25 ◦C) and temperature were conducted using an IntelliCAL PHC101 Standard gel filled pH probe and an IntelliCAL CDC401 Standard Conductivity probe, both connected to a HACH HQ40d instrument. Field measurements of dissolved oxygen (O2) were conducted using a WTW Optical IDS FDO® 925 sensor connected to a WTW oxi-3310 IDS Portable Dissolved Oxygen. The probes were in all cases

immersed in a 0.3-L beaker flushed and freshly filled with sample water, delivered from the piezometer by the peristaltic pump.

#### *2.6. Hydraulic Properties*

The horizontal hydraulic conductivity at the well screen locations was estimated by the Hvorslev straight-line method, using specific software (AQTESOLV/Pro). It assumes a quasi-steady-state flow by omitting the storativity of the aquifer. After conversion into normalized response data (H/H0), the Hvorslev solution for unconfined partially penetrating wells was used [48]. While developed for confined aquifers, the Hvorslev method can also be used for unconfined conditions, provided that the well screen is not too close to the water table [48]. The method assumes that the aquifer is homogeneous and of uniform thickness, and can be used in both cases of fully or partially penetrating wells. The test was repeated up to four times, depending on the time to get one slug test completed, and the results presented here are the arithmetic mean for each location.

#### **3. Results**

#### *3.1. Physical Tracers of Flow (EC, Stable Isotopes of Water and Groundwater Seepage)*

The data of sampling points 1 to 30 (Table S1) were projected into a cross-section (EC, D-excess, and hydraulic conductivity) to facilitate the interpretation and comparison. Some data are overlapped since they were located at the same distance from the shoreline; however, it is also a representation of the spatial variability in the study area.

Groundwater EC showed a consistent pattern starting from freshwater and transitioning to saltwater moving offshore (July 2017 campaign). The change from fresh to saltwater takes place in a mixing zone of 3–5 m width. Interpolation of the measurements at different depths and distances from the shore showed the shape of a classic saline wedge (in spite of representing only the top four meters of the aquifer), where higher density water intruded freshwater (Figure 2). The EC ranged from 0.15 to 20.20 mS cm<sup>−</sup>1. The lowest EC corresponded to the fresh groundwater inland, and the highest corresponded approximately to the maximum EC of the lagoon during a year [38]. The measured profiles of EC in Inland wells A and B tapped the freshwater part of the aquifer (August 2016 campaign). EC fell in the range of 0.115–0.386 mS cm−<sup>1</sup> (Inland well A) and 0.215–0.274 mS cm−<sup>1</sup> (Inland well B), confirming that the inland water was fresh. The high density of measurements and the absence of outliers indicated that the observed patterns are apt for the further interpretation of chemical and hydrogeological processes.

**Figure 2.** Interpolated electrical conductivity (EC) along the cross-section and sampling points (vertical scale × 2.5). The dashed line represents the transition from freshwater to saltwater. Sampling points are located with empty circles, and the locations of sampling points 1, 10, 20 and 30 are marked.

Seepage rates were spatially variable (Figure 3), due to changes in hydraulic properties and local hydrodynamics. There were a few points in the proximity to the shoreline and in the beginning of the saltwater area with low fluxes near zero. Still, in the majority of the locations, the seepage was from the aquifer to the lagoon. The observation of spatial variability of flux is frequent in seepage meters studies [49]. This was also confirmed in this study with fluxes ranging between ~0–10 cm/d in a 40-m transect. Within distances of 1–3 m, fluxes changed up to 10-fold. The EC that was estimated, with the mass balance of the seepage meter bags (Equation (1)) that were higher or lower than 10 mS/cm, correlated with measurements of saline/fresh groundwater below the seepage meters. Thus, a clear pattern of groundwater becoming more saline in the offshore direction was evident, from both seepage meter estimated EC and measured subsurface EC (Figure 3).

**Figure 3.** Spatial plot of discharging fluxes measured with seepage meters and an interpolation of measured EC in collected groundwater. The red/purple colors of the flux bars indicate if the mass balance of the seepage meters bags showed fresh (<10 mS/cm) or saline water (>10 mS/cm).

The deuterium excess (D-excess) (Equation (2)) was used as an indicator of evaporation:

$$\text{D-excess} = \text{ } \text{\ $D-8\$ }^{18} \text{O} \tag{2}$$

An independent spatial differentiation with three distinct zones was identified from the changes in D-excess that correlates with the EC distribution: a freshwater zone and a saline zone separated into shallow and deep zones (Figure 4). To probe the differences between the various zones, a two-tailed t-test was conducted between the D-excess measured in each zone. There was a significant difference (*p*-value < 0.1) in the D-excess content between freshwater and shallow saline samples (*p*-value < 0.001), and also between the samples from the two shallow and deep saline zones (*p*-value < 0.001). The freshwater zone was characterized by high D-excess and uniform mean and standard deviation values: −10.21 + 1.23%, maximum: −13.09%, n = 85). The samples from the shallow saline zone had a much lower D-excess (−4.76 + 1.61%, maximum: −9.37%, n = 32). In the deep saline zone, the D-excess was the lowest (−1.22% + 2.68%, maximum: −6.13%, n = 27). A dipping line marks the differences between the shallow and the deep saline zone (Figure 4). In the furthest offshore location, the shift took place at 2 m below the lagoon bed, while at the location of the saltwater–freshwater interface, this shift was at 4-m depth. The spatial distribution of the isotopic tracer was in good agreement with the EC distribution (Figure 2).

**Figure 4.** Deuterium excess along the cross-section with dashed lines dividing data into a freshwater zone, a shallow zone, and a deep saline zone, as explained in the text (vertical scale × 2.5).

Stable isotopes of water (δ18O and δD) in groundwater showed alignment with the Global Meteoric Water Line (GMWL) and the Local Meteoric Water Line (LMWL) [50].

$$\text{LIMWL} = 7.29 \times \delta^{18} \text{O} + 4.81\%\_{00} \tag{3}$$

The <sup>δ</sup>18O values ranged between <sup>−</sup>7.69 and <sup>−</sup>2.75 %, while <sup>δ</sup>D was <sup>−</sup>50.1 and <sup>−</sup>22.6% (Table S1). Many of the samples grouped around values near <sup>δ</sup>18O = <sup>−</sup>7.6% and <sup>δ</sup>D = <sup>−</sup>47%, corresponding to freshwater (Figure 5). These values were similar to groundwater in other regions of Denmark [50–53]. There was no evidence of fractionation of freshwater samples, as all grouped together with no outliers. A conservative mixing line was established between the groundwater endmember and oceanic seawater (i.e., VSMOW standard), in which part of the samples aligned well. Additionally, an evaporation line was added, where the start and the slope were selected based on the last freshwater sample presence, and following the trend of the samples with lower δD, respectively (Figure 5). The enrichment in heavy isotopes for samples with higher EC was indicated by the alignment between the conservative mixing line with seawater and the evaporation line (Figure 5).

**Figure 5.** δ18O and δD of groundwater samples with the Global Meteoric Water Line (GMWL), the Local Meteoric Water Line (LMWL), a mixing line between fresh groundwater in the study area, and the sea (δ18O and δD: 0%, 0%), and a potential evaporation line. Colors indicate the sample categories from Figure 4.

Hydraulic conductivity (K) along the transect ranged between 0.02–28.0 m d−1. The average values for each depth were similar: 9.3 m d−<sup>1</sup> (1 m), 7.4 m d−<sup>1</sup> (2 m), 9.3 m d−<sup>1</sup> (3 m), and 9.8 m d−<sup>1</sup> (4 m). In spite of the apparent homogeneity, there were spatial differences in K. At 1-m depth, K was relatively homogenous (standard deviation of 4.5 m d<sup>−</sup>1) than at deeper locations (standard deviation at 2 m = 8.3 m d−1,3m = 7.5 m d−1, and 4 m = 10.9 m d−1). At 2-m, 3-m, and 4-m depths, K values were more variable, especially in the zone dominated by saltwater (Figure 6). In some of the locations, measurements could not be completed, reflecting low K or field technical problems (Figure 6).

**Figure 6.** Hydraulic conductivity at the different depths measured, and location of the freshwater– saltwater transition in groundwater (vertical scale × 2.5).

#### *3.2. Hydrochemistry of Groundwater*

Inland wells A and B present different hydrochemical characteristics despite the short distance between them (50 m). In Inland well A, NO3 concentrations reached a maximum of 772 μmol/L or 47.8 mg/L at a 6-m depth (Figure 7), indicating the presence of a NO3 plume. The EC profile follows this trend; thus, the uncontaminated freshwater is probably reflected by the most surficial measurement (~0.115 mS/cm), showing infiltrating rain water near the coast. In Inland well B, no NO3 is observed, except at the depth of 7 m (7 μmol/L). This indicates that between Inland well A to Inland well B, mixing and geochemical reactions occur, giving a relative uniform NO3 profile at Inland well B. Alkalinity showed an increase from Inland well A to Inland well B, suggesting that denitrification with organic matter took place (notice that Inland well A only has three measurements). The alkalinity in Inland well B compared with the concentration observed at locations 1–30 (data not shown). Groundwater samples from the 2017 campaign (sampling points 1–30 at multiple depths) also showed no presence of NO3 except for five samples that were below the detection limit.

For PO4, Inland wells A and B had average values of 0.8 and 3.9 μmol/L, respectively, which is a notable increase in PO4 moving toward the shore of the lagoon. In Inland well A, all the sampling points, except at 3-m depth, did not show any PO4. In Inland well B, the highest concentrations were close to the topographic surface, decreasing to a minimum at 5–6 m, followed by an increase at deeper locations (Figure 7).

PO4 distribution in groundwater along the cross-section under the lagoon showed differences that also allowed dividing the transect into three zones, by grouping the samples plus a fourth category (Figure 8) (Table S1). In the freshwater zone, PO4 had an irregular distribution with the highest values isolated between much lower concentrations (Group I). The mean concentration for this zone (n = 85) was 1.33 μmol/L with a maximum of 33.28 μmol/L, and a standard deviation of 3.89 μmol/L. At the saltwater–freshwater interface, four samples showed a higher concentration, which were assigned

to a different category (Group II) (arithmetic mean = 14.44 μmol/L, standard deviation = 5.86 for n = 4). In the saline zone, the arithmetic mean at variable depths showed marked differences (depth 0.2 m = 5.95 μmol/L,1m = 1.97 μmol/L, 2 m = 0.11 μmol/L, 3 m = 0.13 μmol/L, and 4 m = 0.09 μmol/L). Thus, higher PO4 concentrations were measured at 0–1 m depth, and supported a subdivision of the saline zone into shallow and deep zones. The shallow saline zone (Group III) (Figure 8) had higher PO4 concentrations (arithmetic mean 3.87 μmol/L for n = 21) with a maximum of 8.13 μmol/L, and a standard deviation of 2.52 μmol/L). For the deeper saline zone (Group IV) (Figure 8), the arithmetic mean, maximum, minimum, and standard deviation were much lower in comparison (0.11 μmol/L for n = 31, maximum 0.38 μmol/L, and standard deviation 0.10 μmol/L). The differences between groups were verified with a two-tailed t-test probing the differences in PO4 concentrations in groundwater. Groups I and II showed significant differences (*p* = 1.97 <sup>×</sup> 10−2) as well as Groups II and III (*p* = 3.41 <sup>×</sup> 10−2) and Groups III and IV (*<sup>p</sup>* <sup>=</sup> 1.20 <sup>×</sup> <sup>10</sup><sup>−</sup>2).

**Figure 7.** Profiles of measured EC, alkalinity, PO4, and NO3 for Inland wells A and B.

**Figure 8.** Bubble diagram of PO4 concentrations measured in groundwater along the transect, together with interpolated EC. PO4 concentrations have been plotted with different colors for the categories, as explained in the text (vertical scale × 2.5).

#### **4. Discussion**

#### *4.1. Salinity and Drivers of Fluxes*

Groundwater discharge to the lagoon takes place not only in areas dominated by freshwater at the subsurface, but also in areas dominated by saltwater at the subsurface. This was inferred from the patterns of the subsurface EC measurements (Figure 2) that corroborated with the seepage meter calculated EC (Figure 3), as well as groundwater stable isotopes' signals (Figure 4). Often, the conceptualization of groundwater discharge to lagoons includes only the freshwater component, as this can be estimated based on hydraulic gradient and estimated K value. In addition, it is frequently anticipated that the contamination sources originated inland [26]. However, these assumptions, based on the dataset presented in this work, can not alone explain our observations. There are fluxes that take place in the zone occupied by saltwater, which can have an important contribution to the chemical budgets of coastal areas and lagoons. The drivers of these fluxes are related to well-known variable density flow [21], and are also due to waves and water level oscillations of the surface water [23], as was demonstrated indirectly in several marine studies from the 1990s [15,54], and more specifically in recent studies [55–57].

In the study area, the differentiation between fresh and saline groundwater allowed the identification of different hydrochemical processes. While in the freshwater sector NO3 concentration was decreasing, PO4 was generated at different ratios, depending on the depth. Hence, the fresh SGD would have a minor impact on the nutrients discharge, while the saline SGD would be more relevant in this study area.

One major question is whether these results can be used to estimate nutrient budgets for the full lagoon by extrapolation. In natural environments, the spatial variability and the differences in nutrient sources require caution with this type of approach. In this study area, it can be assumed that the freshwater discharge is limited to a fringe of 10–20 m from the shoreline of the lagoon along all of the east coast, due to the similar hydrogeological setting (similar surrounding topography, hydrogeology, and precipitation regime). In terms of the saline discharge, it was measured in this study up to 40 m from the shore, and in principle, could be active even further offshore considering that: (i) The depth of the lagoon from the further offshore locations does not exceed 1.9 m [58], so similar conditions regarding depth can be assumed; (ii) the sediment characteristics at other offshore locations of the lagoon are unknown, though in small exploration campaigns a few hundred meters offshore, no changes in the sediment composition were observed; and, (iii) the presence of waves and lagoon level changes are common to the entire lagoon area. Under these assumptions, it is reasonable to assume that most of the lagoon, except for the 10–20 m wide shoreline fringe, is dominated by saline fluxes generated by processes not related to the inland hydraulic conditions. Other considerations, such as windows in deep confined aquifers [37], would be exceptions to this conceptualization.

The measured flux in seepage meters and estimated EC of groundwater that were used to characterize the lagoon and groundwater interactions agreed well, showing that one of the methods would be enough to inform about the differentiation between fresh–salt groundwater areas. Still, there are uncertainties when using seepage meters to estimate the EC of groundwater, due to the assumptions that are required: The seepage meter chamber has to be totally flushed out of groundwater before the first measurement, with no exchange between the chamber and the seepage meter bag during the measuring time, and where the groundwater seeping into the seepage meter mixes immediately in the chamber and can be collected in the seepage meter bag. Therefore, the variability in fluxes estimated with seepage meters, combined with the concentration of chemical components immediately below the lagoon bed (to estimate the mass load to Ringkøbing Fjord), should be considered with a range of variability that can only be provided if multiple seepage meters are used in field campaigns. The direct sampling of seepage meters can generate changes in the redox conditions of the samples [59].

#### *4.2. Stable Isotopes of Water*

Groundwater in coastal areas is exposed to different processes that affect its isotopic composition. Salinity (Figures 2 and 4) and isotopic signal (Figure 3) mapping along the cross-section supported the differentiation of fluxes, flow paths, and mixing of waters with different origin and physicochemical characteristics (i.e., salinity, oxygen, and content of nutrients). Overall, the LMWL matched the isotopic composition of groundwater. This is more evident for the freshwater samples, as the saline samples are likely affected by evaporation (Figure 5). Thus, the mixing line with saltwater fits well with part of the group of samples located in the shallow saline zone of the aquifer, while the other part indicated mixing between the lagoon water and groundwater. This mixing process can be explained by the benthic exchange occurring in the saline zone, as it affects only the shallow part of the aquifer where, for example, pressure oscillations (waves and surface water level changes) will have a reduced magnitude and would not encroach deeper into the aquifer [23].

The deviation in the plot of the samples for δ18O and δD (Figure 5) can be explained by an evaporation process. Under the assumption that evaporation explains the isotopic composition of these samples, the reason for being located at a deeper location (deep saline water in Figure 4), can be associated to seasonal flow dynamics, due to the lagoon density changes (i.e., salinity in the lagoon is higher during the summer and would move deeper into the aquifer due to its higher density). In principle, evaporation would both increase the salinity as well as the concentrations of deuterium and δ18O, but in fact, each per-mille increase in δ18O would result in only a few-percent increase in salinity.

Therefore, the use of water stable isotopes is an excellent method for identification of the different type of fluxes driving SGD. The main challenge is the high density of data required to determine the isotopic signature of each flux. At our field site, sampling at a meter scale was required, but in other natural systems, it might require an even higher density of data or a higher resolution at variable depths. The main advantage of water stable isotopes is the simplicity in the sampling, and the reliability of the results.

#### *4.3. Nutrients*

The distribution of PO4, along the cross-section (Figure 8), allowed the distinction of four groundwater categories with different chemical composition. The spatial variability in the freshwater zone (category I), with no clear pattern, could be due to the natural spatial variability of the sediment composition by the dissolution of carbonates, which may contain adsorbed or co-precipitated phosphate. Unfortunately, there are no data on sediment composition to confirm this hypothesis. In the saline zone, the three categories seemed to be supported by the hydrogeological and isotopic measurements. The transition from freshwater to saltwater was characterized by an increase in the PO4 content. This can be associated with the displacement by anions in seawater, and the release of P and N sorbed to Fe oxides from organic matter, as has been reported in oxic aquifers [60] and lab experiments [61].

The differentiation between a shallow and a deeper saline zone is also accompanied by changes in the PO4 distribution (Figure 8). The coincidence between the chemical composition of PO4 in groundwater and the isotopic/EC distributions could be an indicator of geochemical processes, which are triggered by the mixing of water with different oxic characteristics. Following this reasoning, the saltwater–freshwater interface would be characterized by the mixing of anoxic waters with different salinities, where the shallow part of the aquifer will be in contact with lagoon water enriched in oxygen, while the deeper zone will have more stable anoxic conditions. Additionally, changes in pH can also play a role. This is a complex environment, where the observations reported can be generated by different processes. We discuss three potential hypotheses that may explain the distribution of PO4:

(1) Sediment composition differences. The horizontal layering of the changes observed in groundwater is concordant with the deposition of sediments. The reason that the most surficial sediments have a higher content of PO4 would be associated with the changes that have taken place in Denmark for the last century. The increase in the use of fertilizers and the sediment transport by the Skjern River that drains an agricultural catchment support the development of recent layers of sediments enriched in PO4. This interpretation supports the observed changes in PO4 in the top 1 m of the aquifer. Additionally, the hydraulic tests indicate a change in the characteristics of the sediments between 1–2 m depths (Figure 8). Nevertheless, this difference was not observed in the freshwater zone. Additional data about PO4 sediment composition is required to confirm this hypothesis.


Most likely, the distribution of PO4 monitored in the study area is a coalescence of the different processes discussed above. The contribution of each of these processes would require further specific research and sampling campaigns.

The presence of NO3 in Inland well A showed contaminant sources in the upland area (Figure 7). The origin of NO3 is probably associated to the septic tanks used in the summer houses (between 200–300), and the diffuse input due to fertilizer in agricultural fields. The difference with Inland well B (Figure 7) indicated that geochemical reactions can be responsible for the fast degradation of NO3. Denitrification in the study area can be favored by the presence of organic matter in the sediments and the mixing of water with different redox conditions, as it has been shown in intertidal beaches [67]. The two wells are aligned along what would be expected to be a flow line, but an alternative hypothesis is that Inland wells A and B do not represent the same water, but different flow paths. Still, it is clear that NO3 transported by groundwater is not reaching the lagoon, as only minor concentrations were measured (albeit below the detection limit) in a few of the 150 sampling points collected under the lagoon. Under the lagoon bed, the presence of organic matter that can contribute to the denitrification processes has been documented, especially in the proximity of the shoreline of the lagoon [68].

#### **5. Conclusions**

SGD at the field site shows a clear distinction between the hydrodynamic characteristics and processes taking place in the freshwater-dominated zone and the saline zone. The drivers of fluxes, the origin of discharging water, the chemical reactions that can take place, and the type of mixing of water must be considered in the analysis and methods for studying nutrients borne by SGD.

Stable isotopes of water can be used as a tracer of hydrogeological flow paths in coastal areas. The simplicity in sampling stable isotopes of water, and the relative low cost of the analysis, point to

them as an efficient and useful tool for the differentiation of fluxes in SGD, which poses a persistent challenge in the study of coastal areas.

The spatial distribution of PO4 detected in Ringkøbing Fjord can be explained by: (1) Changes in the sediment composition of the study area, (2) mixing and chemical reaction of water with different origin and composition, due to fluxes driven by variable sources (freshwater, density driven flow, benthic exchange), or (3) diffusive flux and biological activity in the lagoon bed. It may be that all processes are taking place simultaneously, but with different contributions.

Groundwater originating inland has low PO4 content and high NO3 content in the study area. The processes taking place in the proximity of the lagoon, as denitrification, reduce NO3 concentrations before discharge to the lagoon, while PO4 can be generated in the aquifer below the lagoon, even if it is not detected in the inland wells. These two processes revert the input of nutrients that could be inferred by simple measurements of groundwater inland.

Mixing of water with different salinities and origin can trigger reactions generating and/or degrading chemical components at multiple spatial and temporal scales in coastal aquifers. Addressing salinity distribution and the distinction of flow paths will help to develop better monitoring strategies, and more accurate models to improve the estimation of nutrient budgets due to SGD in coastal areas.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/11/9/1842/s1, Table S1: In situ measurements, stable isotopes, nutrient and ion analysis for sampling points 1 to 30.

**Author Contributions:** Conceptualization, C.D. and P.E.; Formal analysis, C.D. and S.J.; Funding acquisition, C.D. and P.E.; Investigation, C.D., J.T.-C., S.K. and P.E.; Methodology, C.D., S.J. and P.E.; Project administration, C.D. and P.E.; Visualization, C.D. and S.J.; Writing—original draft, C.D.; Writing—review and editing, C.D., S.J., J.T.-C., S.K. and P.E.

**Funding:** The research leading to these results has received funding from the European Union Seventh Framework Programme FP7/2007-2013 under grant agreement number 624496, the European Union's Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie Grant Agreement number 722028 and the HOBE project (Hobe.dk) funded by the Villum Fundation. This work was carried out as part of the activities of the Aarhus University Centre for Water Technology, WATEC.

**Acknowledgments:** We thank Iris Tobelaim, Christian Hansen and students in the Hydrogeological field course (2016) for field and laboratory assistance.

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

#### **References**


© 2019 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/).
