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

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 (¹⁸O and ²H) 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 ¹⁸O and ²H 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. PO₄ underneath the lagoon exceeds the groundwater concentration inland, thus demonstrating an intra-lagoon origin, while NO₃, higher inland due to anthropogenic activity, is denitrified in the study area before reaching the lagoon.

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,4,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,11,12]. Subsequent investigations on submarine groundwater discharge have used seepage meters, radioactive tracers, or hydraulic methods to estimate the flux of groundwater [13,14,15]. Flux estimates combined with nutrient concentrations were used to calculate nutrient loads to coastal waters in several regions of the world [16,17,18,19,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,30,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 ¹⁸O and ²H 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 (PO₄ and NO₃) 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 km² 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 NO₃ [39], and is associated with septic tanks and PO₄ [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,42,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].

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⁻¹). 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:

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

(1)

where W_off: Weight of the seepage meter bag after measuring, EC_off: EC of water in the seepage meter bag after measuring, W_on: Weight of the seepage meter bag before measuring, EC_on: EC of water in the seepage meter bag before measuring, and W_bag: 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⁻¹ 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⁻¹). 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 ¾ 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⁺, Ca²⁺, Mg²⁺) and ammonium (NH₄⁺), another 10-mL PE vial for major anions (Cl⁻, SO₄²⁻, NO₃⁻, Br⁻), a 5-mL PE vial for PO₄³⁻ and Fe²⁺ 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 PO₄³⁻ and Fe²⁺ were conserved by 2 vol% 2 M H₂SO₄, samples for major cations and ammonium were acidified by 1 vol% 1 M HNO₃ (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 O₂, 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.

PO₄³⁻ were analyzed spectrophotometrically using the molybdate blue method of Murphy and Riley (1986). Total dissolved Fe, henceforth assumed to a direct proxy for the Fe²⁺ concentration, was measured on residual sample aliquots of the PO₄³⁻ samples after reduction to Fe²⁺ with hydroxylamine hydrochloride, using the spectrophotometric ferrozine method of Stookey (1970).

δ¹⁸O 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 δ¹⁸O 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 (O₂) 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/H₀), 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⁻¹. 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⁻¹ (Inland well A) and 0.215–0.274 mS cm⁻¹ (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.

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

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

$$\mathrm{D-excess}= \mathsf{\delta }\mathrm{D}-8{\mathsf{\delta }}^{18}\mathrm{O}$$

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

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

$$\mathrm{LMWL}=7.29\times {\mathsf{\delta }}^{18}\mathrm{O}+4.81‰$$

(3)

The δ¹⁸O values ranged between −7.69 and −2.75 ‰, while δD was −50.1 and −22.6‰ (Table S1). Many of the samples grouped around values near δ¹⁸O = −7.6‰ and δD = −47‰, corresponding to freshwater (Figure 5). These values were similar to groundwater in other regions of Denmark [50,51,52,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).

Hydraulic conductivity (K) along the transect ranged between 0.02–28.0 m d⁻¹. The average values for each depth were similar: 9.3 m d⁻¹ (1 m), 7.4 m d⁻¹ (2 m), 9.3 m d⁻¹ (3 m), and 9.8 m d⁻¹ (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⁻¹) than at deeper locations (standard deviation at 2 m = 8.3 m d⁻¹, 3 m = 7.5 m d⁻¹, and 4 m = 10.9 m d⁻¹). 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).

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, NO₃ 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 NO₃ 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 NO₃ 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 NO₃ 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 NO₃ except for five samples that were below the detection limit.

For PO₄, Inland wells A and B had average values of 0.8 and 3.9 µmol/L, respectively, which is a notable increase in PO₄ moving toward the shore of the lagoon. In Inland well A, all the sampling points, except at 3-m depth, did not show any PO₄. 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).

PO₄ 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, PO₄ 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, 1 m = 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 PO₄ 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 PO₄ 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 PO₄ concentrations in groundwater. Groups I and II showed significant differences (p = 1.97 × 10⁻²) as well as Groups II and III (p = 3.41 × 10⁻²) and Groups III and IV (p = 1.20 × 10⁻²).

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,56,57].

In the study area, the differentiation between fresh and saline groundwater allowed the identification of different hydrochemical processes. While in the freshwater sector NO₃ concentration was decreasing, PO₄ 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 (Figure 2 and Figure 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 δ¹⁸O 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 δ¹⁸O, but in fact, each per-mille increase in δ¹⁸O 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 PO₄, 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 PO₄ 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 PO₄ distribution (Figure 8). The coincidence between the chemical composition of PO₄ 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 PO₄:

(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 PO₄ 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 PO₄. This interpretation supports the observed changes in PO₄ 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 PO₄ sediment composition is required to confirm this hypothesis.

(2)

Hydrodynamics of the coastal zone observed from the isotopic signature. In this case, the groundwater flow paths and the mixing of water with different properties would be responsible for the PO₄ distribution. The freshwater zone would be characterized by groundwater oxic conditions that do not enhance the mobility of PO₄ [24,25]. The shallow saline zone would be dominated by exchange processes where lagoon water with variable salinity, oxygen conditions, and organic matter generated by primary production would mix with saline, reducing groundwater. These processes can generate the microbial remineralization of organic matter [62]. The increase of PO₄ due to a high exchange of water has been documented in sand flats [63]. The agreement between PO₄ distribution and groundwater characteristics based on δ¹⁸O and deuterium, as well as the measured seepage and the salinity changes, makes this possibility more feasible than sedimentological changes (Hypothesis 1). Groundwater hydrodynamics would control the hydrochemical distribution, rather than the sediment composition changes. The changes in hydraulic conductivity from the first meter to deeper zones can favor the sharp transition in isotopic signal associated to different hydrogeological processes.

(3)

Diffusive and biological processes in the lagoon bed. A well-documented source of PO₄ is the release from sediments by diffusion [64]. Nevertheless, this has been highlighted as a process that only takes place in the top 20 cm [64], while here, the changes were observed further than 1-m depth. Another alternative is associated with the higher temperature due to the season (summer) when it was measured. In the lagoon bed, the high temperature would increase the biological activity, and the releasing of PO_4, as has been documented in lakes [65,66].

Most likely, the distribution of PO₄ 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 NO₃ in Inland well A showed contaminant sources in the upland area (Figure 7). The origin of NO₃ 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 NO₃. 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 NO₃ 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 PO₄ 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 PO₄ content and high NO₃ content in the study area. The processes taking place in the proximity of the lagoon, as denitrification, reduce NO₃ concentrations before discharge to the lagoon, while PO₄ 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.