**Evaluation of Temperature Profiling and Seepage Meter Methods for Quantifying Submarine Groundwater Discharge to Coastal Lagoons: Impacts of Saltwater Intrusion and the Associated Thermal Regime**

#### **Joel Tirado-Conde 1,\*, Peter Engesgaard 1, Sachin Karan 2, Sascha Müller <sup>1</sup> and Carlos Duque <sup>3</sup>**


Received: 2 July 2019; Accepted: 7 August 2019; Published: 9 August 2019

**Abstract:** Surface water-groundwater interactions were studied in a coastal lagoon performing 180 seepage meter measurements and using heat as a tracer in 30 locations along a lagoon inlet. The direct seepage meter measurements were compared with the results from analytical solutions for the 1D heat transport equation in three different scenarios: (1) Homogeneous bulk thermal conductivity (Ke); (2) horizontal heterogeneity in Ke; and (3) horizontal and vertical heterogeneity in Ke. The proportion of fresh groundwater and saline recirculated lagoon water collected from the seepage experiment was used to infer the location of the saline wedge and its effect on both the seepage meter results and the thermal regime in the lagoon bed, conditioning the use of the thermal methods. The different scenarios provided the basis for a better understanding of the underlying processes in a coastal groundwater-discharging area, a key factor to apply the best-suited method to characterize such processes. The thermal methods were more reliable in areas with high fresh groundwater discharge than in areas with high recirculation of saline lagoon water. The seepage meter experiments highlighted the importance of geochemical water sampling to estimate the origin of the exchanged water through the lagoon bed.

**Keywords:** heat as a tracer; temperature profiles; seepage meter; lagoon; coastal areas; seawater intrusion; groundwater discharge; seawater-groundwater interactions

#### **1. Introduction**

Exchange fluxes between groundwater (GW) and surface water (SW) bodies are important due to the hydrological [1] and ecological [2,3] consequences for aquatic plants and animals. Understanding the temporal and spatial heterogeneity of exchange fluxes is important as it can affect water resource management plans. Although GW-SW interactions have been extensively studied [4], it remains a challenge to accurately measure the exchange fluxes in-situ and also decide on where and how often [5].

Lewandowski et al. [6] discussed the historical evolution of studies investigating Lacustrine groundwater discharge (LGD), i.e., exchange fluxes between lakes and groundwater versus that of submarine groundwater discharge (SGD), i.e., exchange fluxes between sea/lagoons and groundwater. The current study is interested in evaluating the methods for quantifying SGD and its inherent spatial heterogeneity. Duque et al. [7] recently showed that spatial heterogeneity in exchange fluxes was more

important than temporal heterogeneity and highlighted the impact of a correct characterization of the spatial heterogeneity on the discharge estimations, so this study focuses on capturing SGD as a snapshot in time.

Seepage meters, such as the ones designed by Lee [8], are the only tool to directly measure LGD and SGD [9]. The use of seepage meters to quantify SGD is often reported [7,10,11], but it is a labor intense methodology and it has certain limitations [9,12]. Other techniques (see review by Rosenberry et al. [5]) that have been extensively applied are the use of groundwater head measurements from piezometers [13] and the use of heat as a tracer [14–16], among others [17], that indirectly derive the horizontal/vertical flux by the application of Darcy's law or by the use of analytical solutions of the 1D heat transport equation. Most of these methods rely on prior information of hydraulic and thermal characteristics of the sediment, which may not be as well-known and predictable in coastal environments as they are in other hydrological systems [18]. In other cases (e.g., [19]), the time series of temperature data at different depths have been used to estimate exchange vertical fluxes. The analytical models for heat transport also require the specification of upper and lower temperature boundary conditions. The upper boundary condition is often set to what was measured at the sediment-water interface, while the lower boundary condition is assumed to be equal to the average groundwater temperature (equal to the average seasonal air temperature) at some depth below the measurement point. While the thermal method has found extensive use for LGD, it has not been used that much to study SGD [18,20,21].

In coastal areas, the GW-SW exchange, or SGD, is influenced by the density-driven flow produced by the differences between the low-density fresh groundwater (FGW) and the higher density saline surface water (SSW). These differences generate a mixing zone with a dynamic behavior that depends on the salinity of both FGW and SSW as well as on the magnitude of the groundwater discharge flux to the coast. Befus et al. [18] summarized a field study using heat as a tracer in a conceptual model of a groundwater-bay system that showed the well-known complex flow with a saltwater interface, a salt water plume at the high tide mark, and a freshwater tube with fresh groundwater discharge. However, they also superimposed a model of the thermal regime, which due to the ~3 ◦C difference in the average annual temperature of groundwater and water in the bay developed into distinct thermal regimes, where the subsurface near the bay was warmer than discharging groundwater. The complexity of the flow field and thermal regime therefore will affect the choice of how to specify boundary conditions for heat transport modeling and possibly also how well two different methods like the seepage meter and thermal methods compare.

In this study, 180 seepage meter measurements in 30 different locations along a coastal lagoon inlet were combined with indirect estimations of the vertical flux using heat as a tracer in those same 30 locations. The proportion of discharging fresh groundwater and recirculated saline lagoon water was used to infer the location of the saline wedge and the consequences that it had on the seepage meter results. The time series of the temperature in deeper wells were used to assess the boundary conditions for the thermal regime of the lagoon sediments that affected the applicability of the thermal methods.

The objectives of this work are therefore to: (1) Test and compare the in-situ direct measurements of groundwater discharge to a coastal lagoon using seepage meters with indirect estimates using heat as a tracer; (2) determine the effect of vertical and horizontal sediment heterogeneity in bulk thermal conductivity on the results obtained by the use of analytical solutions of the heat transport equation; (3) address the impact of lagoon water recirculation on the subsurface thermal regime and on the results of the GW-SW exchange for both the direct measurements and the estimation using heat as a tracer

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

#### *2.1. Study Area*

Ringkøbing Fjord is a coastal lagoon in Western Jutland, Denmark, bordering the North Sea (latitude 55◦58 40 N, longitude 8◦14 21 E), covering an area of 300 km<sup>2</sup> (Figure 1). The mean water depth is 1.9 m and although the deepest points are up to 5 m, more than 25% of the lagoon has water depths of less than 0.5 m. The shallow depths also dominate the study site from the eastern shore line to the elongated island offshore (Figure 1b). The lagoon contains brackish water (5–15%, salinity of 6–15 g L<sup>−</sup>1) due to the connection to the sea to the west (controlled by a sluice) and the seasonal variation of the fresh water discharge mainly from the Skjern River (annual average of 50 m3 s−1). The tides are not an issue [18], although the operation of the sluice and storms can generate small flooding of the shore line. The annual mean precipitation and evaporation are 1050 and 630 mm, respectively [22]. More details of the hydrology can be found in Haider et al. [23], Duque et al. [7,20], and Müller et al. [24].

The study area is one of several human-made inlets on the eastern side of the lagoon (Figure 1c). The geology is predominantly Pleistocene fluvio-glacial sand, including relatively thin layers of silt or lower permeable materials, as well as paleo-channels with higher hydraulic conductivity. At the eastern border, close to the shore, the dense coastal vegetation covers the sandy sediments, creating a layer of organic material that can extend into the fjord for tens of meters [20]. The lagoon bed is composed of homogeneous sand with small wind- and wave-controlled ripples. The orography both in the lagoon and the surrounding areas is mostly flat, with small variations (in the order of 10–20 cm) of the sediment bed topography offshore. The boreholes indicate the existence of an unconfined aquifer of at least 12 m thickness. The hydraulic gradient, obtained from the water level measurements in the boreholes on land and lagoon stage, is 0.005–0.010 (east to west direction) [24].

**Figure 1.** (**a**) The location of Ringkøbing Fjord, (**b**) the study area, the water depth in the lagoon is ~30–50 cm from the shore line and 2 km offshore to the elongated island, (**c**) The lagoon inlet. The black numbered dots indicate the location of the measuring points. The red stars indicate the location of the deep wells used to obtain the average annual groundwater temperature.

#### *2.2. Fluxes from Seepage Meters*

Thirty Lee-type seepage meters [8] were installed on 4 July, 2017 along a 43 m long transect, covering most of the inlet to the lagoon and the first few meters offshore (Figure 1c). Each seepage meter covered an area of 0.25–0.29 m2, and were installed by inserting them around 20 cm into the lagoon bed. The space between the lagoon bed and the top part of the seepage meter was between one and two centimeters to assure that the flow of water to or from the surface was possible. A valve located on one side of the seepage meter allowed connecting a plastic bag to measure the exchange fluxes with the lagoon. The seepage meters were left for several days without attached bags before sampling to allow the brackish lagoon water inside the seepage meters to be displaced by discharging water. The seepage meter experiments were performed six times over a period of two consecutive days, 9–10 July, 2017. These 180 seepage measurements had an average duration of 140 min, with the shortest and the longest experiments of 118 and 175 min, respectively. The plastic bags were pre-filled with 1 L of tap water and weighed before being attached to the seepage meters. After the experiment, the plastic bags were weighed again. The electrical conductivity (EC) of the water in the bags attached to the seepage meters was measured before and after the experiments. The EC was used to approximate

the total dissolved solids (TDS) and the density (ρ) [25] of the water in the bag in order to convert from mass to volume. The exchange flux was then calculated based on the difference in the volume divided by the experimental duration. The exchange fluxes are here reported as a Darcy flux, q = Q/A, where Q is the measured volumetric flux (as described above) and A is the area of the seepage meter.

The EC of the tap (initial water in the bags) and the lagoon water (end members) were used to discern between the fresh groundwater component of the exchange flux and the component corresponding to the brackish lagoon water recirculation in the exchange fluxes. This was based on a mass balance:

$$V\_{\theta}EC\_{\theta} + V\_{\text{ex}}EC\_{\text{ex}} = V\_{f}EC\_{f} \tag{1}$$

$$V\_{ex} = V\_f - V\_o \tag{2}$$

where *Vo* and *ECo* are the initial water volume (L) with its corresponding EC (mS cm<sup>−</sup>1) in the plastic bag prior to the experiment, *Vex* and *ECex* are the volume and the EC of water exchanged during the duration of the experiment, and *Vf* and *ECf* are the volume and the EC of water in the bag after running the experiment. The only unknown is *ECex*:

$$EC\_{c\chi} = \frac{V\_f EC\_f - V\_o EC\_o}{V\_{c\chi}} \tag{3}$$

Then, the mass balance of the exchanged water was calculated:

$$V\_{c\mathbf{x}} E C\_{c\mathbf{x}} = V\_{\mathcal{S}^{\mathbf{w}}} E C\_{\mathcal{S}^{\mathbf{w}}} + V\_{lw} E C\_{lw} \tag{4}$$

$$V\_{lw} = V\_{ex} - V\_{\text{gw}} \tag{5}$$

$$V\_{\mathcal{S}^w} = \frac{V\_{cx}(EC\_{cx} - EC\_{lw})}{EC\_{\mathcal{S}^w} - EC\_{lw}} \tag{6}$$

where *Vgw* and *ECgw* are the volume and the EC of the groundwater exchanged during the experiment and *Vlw* and *EClw* are the volume and the EC of recirculated lagoon water exchanged during the experiment. The relation between *Vgw* and *Vex* and between *Vlw* and *Vex* yield the groundwater component of the exchanged volume of water (%GW) and the recirculated brackish lagoon water component of the exchanged volume of water (%LW), respectively:

$$\% \text{GW} = \frac{V\_{\text{gra}}}{V\_{\text{cc}}} \times 100 \tag{7}$$

$$\% \text{LW} = \frac{V\_{lv}}{V\_{cc}} \times 100\tag{8}$$

Or, in terms of EC:

$$\% \text{GW} = \left(1 - \frac{EC\_{lv} - EC\_{vx}}{EC\_{lv}}\right) \times 100\tag{9}$$

$$\% \text{LW} = \frac{EC\_{lv} - EC\_{ex}}{EC\_{lv}} \times 100\tag{10}$$

#### *2.3. Fluxes from Temperature Profiling*

The vertical temperature distribution in the lagoon bed next to each seepage meter was measured by inserting a probe 0.5 m vertically into the lagoon bed. The probe had 10 temperature sensors (thermocouples with an accuracy of ±0.2 ◦C) located at different depths (0, 2.5, 5, 7.5, 10, 15, 20, 30, 40 and 50 cm from the sediment-water interface). The temperature was recorded after an equilibration period of 10 min based on previous field experiments [20]. Two analytical solutions were fitted to the vertical sediment temperature profiles. This process was performed manually to account for possible measurement errors or inconsistencies. The fitting of the analytical solutions to the temperatures obtained from the deepest temperature sensors was prioritized as they are located further from the lagoon water. The effects of the diurnal temperature variations near the lagoon bed were thus minimal.

Bredehoeft and Papadopulos [26] proposed a steady-state analytical solution (BP) for estimating vertical fluxes from the temperature-depth profiles:

$$T(\mathbf{z}) = T\_s + \left(T\_\mathcal{J} - T\_s\right) \frac{\exp\left[N\_{\rm pc}\left(\frac{\tilde{z}}{L}\right)\right] - 1}{\exp\left[N\_{\rm pc}\right] - 1} \tag{11}$$

where T(z) is the (steady-state) temperature at a given depth *z*(m) in degrees Celsius, *TS* is the (steady-state) lagoon water temperature in degrees Celsius, *Tg* is the temperature of the groundwater at a given depth *L*(m) in degrees Celsius, and *Npe* is the Peclet number representing the ratio of convection to conduction:

$$N\_{\rm pc} = \frac{q\_{\rm z} \rho\_f c\_f L}{K\_{\rm c}} \tag{12}$$

where *qz* is the vertical fluid flux or exchange flux (m s<sup>−</sup>1), ρ*<sup>f</sup>* is the density of the fluid (kg m<sup>−</sup>3), *cf* is the specific heat capacity of the fluid (J kg−<sup>1</sup> ◦C<sup>−</sup>1) and *Ke* is the bulk thermal conductivity (W m−<sup>1</sup> ◦C<sup>−</sup>1).

The BP was applied considering two different scenarios: (1) Homogeneous sediment composition with respect to *Ke* in the study area (Figure 2a); (2) heterogeneous *Ke* (Figure 2b) due to the different organic matter fractions in the lagoon bed, based on the observation of decreasing organic matter in the top sediment offshore. To reproduce these two scenarios, two different parameterizations were used. The values of the parameters considered for the homogeneous scenario were chosen assuming a sandy sediment (Table 1).


**Table 1.** The parameters used for the application of a steady-state analytical solution (BP).

The parameters used for the heterogeneous scenario were based on the work of Duque et al. [20] and are also listed in Table 1. Duque et al. [20] observed that the thermal conductivity in the top part of the lagoon bed was not homogeneous. It is distinguished between the different types of sediments along the fjord inlets related to the amount of vegetation present. Close to the shore, the vegetation was abundant and the thermal conductivity of the top part of the lagoon bed was low (0.72 W m−<sup>1</sup> ◦C<sup>−</sup>1). Further offshore, the amount of vegetation decreased until eventually disappearing and the thermal conductivity of the top part of the lagoon bed increased until reaching a value typical for sandy sediments (1.82 W m−<sup>1</sup> ◦C<sup>−</sup>1). The lagoon water temperature (*Ts*) was obtained from the top temperature sensor in the vertical probe, located at the lagoon water-sediment interface. Unlike the remaining parameters, the groundwater temperature (*Tg*) was based on the mean annual groundwater temperature collected in previous studies from the different wells within the study area (Figure 3). This value was chosen as it represents the overall thermal behavior of deep groundwater in the study area. As the BP solution assumes a steady state, it is necessary to choose a groundwater temperature

with minor variability during the measurement period, a condition that could not be met if using groundwater temperature at shallow depths. Furthermore, as temporal variation of exchange fluxes in the study area are small [7], the use of the mean annual groundwater temperature provides a robust lower boundary condition for this method, representing the average thermal characteristics of the system. The density of the fluid was assumed to be dominated by discharging fresh water as the area is close to the shore and the measured EC values were relatively small.

**Figure 2.** A schematic cross section of the study area with different zonation and parametrization for the lagoon bed temperature analysis. (**a**) BP homogeneous soils scenario with a single thermal conductivity, (**b**) BP heterogeneous soil scenario. Three zones are defined, each of them with a different thermal conductivity, (**c**) SB multi-layered. See text for the definition of the two BP and SB analytical solutions. Three different zones and two layers per zone are defined. Each zone has a different thermal conductivity for the upper layer. All zones have the same thermal conductivity for the lower layer.

The vertical heterogeneity in thermal conductivity (layering) was likely to be present on site [27] and its effect on the estimation of fluxes were also assessed. Previous thermal conductivity measurements mapped the top 10–15 cm of the lagoon bed [20]. It is the hypothesis that the thermal conductivity at greater depths is more homogeneous as sediments are less affected by organic matter depositions which originated from the root activity and plant debris, and the consolidation increases with depth, decreasing porosity [27]. Accounting for the layering of the thermal conductivity, the analytical solution by Shan and Bodvarsson [28] was applied.

Shan and Bodvarsson [28] modified the Bredehoeft and Papadopulos [26] method to allow for the introduction of different thermal conductivity values representing a multi-layer system (SB):

$$T\_i(z) = \mathbb{C}\_{i,1} \mathfrak{e}^{q\_i z/a\_i} + \mathbb{C}\_{i,2} \text{ ( $i = 1, 2, \dots, n$ )}\tag{13}$$

where *Ti(z)* is the temperature at a given depth z in layer *i*, *Ci*,1 and *Ci*,2 are two integral constants, *qz* is the vertical fluid flux across all the layers and α<sup>i</sup> is the thermal diffusivity of the *i*th-layer, defined as the ratio between the effective thermal conductivity (*Ke*) and the heat capacity (ρ*<sup>c</sup>* <sup>=</sup> 4.18 <sup>×</sup> 106 <sup>J</sup>/m3/ ◦K). It is important to note that, as described in Kurylyk et al. [29], Shan and Bodvarsson [28] incorrectly defined thermal diffusivity. Here, α*<sup>i</sup>* is obtained from the volumetric heat capacity of the water (ρ*c*) and should not be confused by the volumetric heat capacity of the medium. By setting the origin at the surface of the top layer, the depth z increases downward and the base of each layer at a depth di, the thickness of each layer is the difference of its two boundary coordinates *bi* = *di*−*di*−1, where *d*<sup>0</sup> = 0. *Water* **2019**, *11*, 1648

Assuming constant temperatures at the top and the bottom of the system (*Ts* and *Tg* respectively), the boundary conditions are:

$$T\_1(0) = T\_s; \ T\_n(d\_n) = T\_{\mathcal{R}} \tag{14}$$

$$T\_i(d\_i) = \ T\_{i+1}(d\_i) \ (i = 1, \ 2, \ \dots, \ n - 1) \tag{15}$$

$$a\_i \left(\frac{dT\_i}{dz}\right)\_{z=d\_i} = a\_{i+1} \left(\frac{dT\_{i+1}}{dz}\right)\_{z=d\_i} (i = 1, \ 2, \ \dots, \ n - 1) \tag{16}$$

Equation (14) are the top and the bottom temperature boundary conditions, (15) is the temperature continuity condition at the layer interfaces, and (16) is the thermal flux continuity condition.

The two integral constants are:

$$\mathbb{C}\_{i,2} = \frac{aT\_s - T\_{\mathcal{S}}}{a - 1} \text{ ( $i = 1, 2, \dots, n$ )}\tag{17}$$

$$\mathcal{C}\_{1,1} = \frac{T\_{\mathcal{S}} - T\_s}{a - 1} \tag{18}$$

$$\mathbb{C}\_{(i+1),1} = \mathcal{e}^{q\_i d\_i (\frac{1}{a\_i} - \frac{1}{a\_{i+1}})} \mathbb{C}\_{i,1} \ (i = 1, \ 2, \ \dots, n - 1) \tag{19}$$

where the parameter *a* is defined as:

$$a = e^{q\_z d\_n / a\_{eff}} \tag{20}$$

where *dn* is the total thickness of the *n* layers and αeff is the effective thermal diffusivity of the *n* layers:

$$a\_{eff} = d\_n / \sum\_{i=1}^{n} \left(\frac{b\_i}{\alpha\_i}\right) \tag{21}$$

**Figure 3.** A schematic (not to scale) cross-section of the study area. In summer, fresh cold groundwater discharges into the lagoon from the east, creating a mixing zone when it encounters the saline slightly warmer lagoon water (a saline wedge is generated due to the density difference). The EC measurements east of the mixing zone showed lower EC values than the measurements west of the mixing zone. The mean annual temperature (9.8 ◦C) measured in two deep wells (4 m depth, black vertical lines, see Figure 1c) was used as a lower boundary condition for the BP method.

Following Duque et al. [20], three different zones were defined for the application of SB, each representing a 2-layer system with a lower thermal conductivity top layer and a higher thermal conductivity bottom layer (Figure 2c). The first zone, with a top layer thermal conductivity of 0.72 W m−<sup>1</sup> ◦C<sup>−</sup>1, corresponds to the initial 15 m from the shore (points 1 to 11, Figure 1), the second zone, with a top layer thermal conductivity of 1.46 W m−<sup>1</sup> ◦C<sup>−</sup>1, corresponds to the next 5 m (points 12 to 17) and the third zone, with a top layer thermal conductivity of 1.82 W m−<sup>1</sup> ◦C<sup>−</sup>1, corresponds to the rest of the points located more offshore (points 18 to 30). The thermal conductivity value chosen for the bottom layer was 1.84 W m−<sup>1</sup> ◦C−1, representing sand. The last zone therefore has a nearly uniform thermal conductivity. The temperatures measured by the top and bottom sensors of the vertical probe were used to fix *Ts* and *Tg* respectively. The location of the interface between the two layers (the base of the top layer) was based on the best possible fit of the SB solution.

#### **3. Results**

#### *3.1. Seepage Meter Fluxes*

The differences between estimated fluxes using the density of each sample compared to using a constant density (1000 kg m<sup>−</sup>3) were negligible (lower than <sup>±</sup> 0.1 cm day−<sup>1</sup> for 156 of the 180 samples). Thus, only the results obtained using the same density for all the samples are shown and discussed in the following. The fluxes measured during the second day were generally higher (25 of the 30 points for rounds 4 to 6). A trend in the variability of the measured fluxes at each point indicated an increase in the standard deviation when the distance to the shore increased (Figure 4). The maximum and minimum standard deviations (2.65 cm day−<sup>1</sup> and 0.11 cm day−1) correspond to locations 24 and 22 respectively, with a separation of only 5 m. The arithmetic mean seepage fluxes (red diamonds, Figure 4) are used to study the spatial variability of the exchange fluxes, Figure 5a.

**Figure 4.** The exchange fluxes corresponding to the seepage meter experiments for each round of measurements at each measuring point. The red diamond corresponds to the mean value. The black dashed line shows the mean flux for all locations and all rounds. The blue dashed line points to the location from which no fresh water was found in the water samples to the left of this line (see Section 3.4).

**Figure 5.** The spatial distribution of the measured and calculated exchange fluxes: (**a**) The mean from seepage meters, (**b**) calculated with the Bredehoeft and Papadopulos (1965) method for the homogeneous scenario, (**c**) calculated with the Bredehoeft and Papadopulos (1965) method for the heterogeneous scenario and (**d**) calculated with the Shan and Bodvarsson (2004) method.

The exchange fluxes ranged from <sup>−</sup>0.02 to 10.22 cm day−<sup>1</sup> with an overall mean of 3.44 cm day−<sup>1</sup> (Figure 4) and a standard deviation of 2.63 cm day−1. Further, 28 of the 30 points showed a positive groundwater–lagoon water exchange flux, indicating a predominant process of groundwater discharge to the lagoon. The two negative exchange fluxes, indicating groundwater recharge, were very small (−0.02 cm day−1) compared with the rest of the measured fluxes. There is not an obvious pattern in the distribution of the discharge fluxes (Figure 5a), but two different higher discharge zones are distinguished. The first is located 18 m offshore (points 12 to 14 and 17) and the second is located from 30 to 43 m offshore (points 21, 23, 24 and 26 to 30).

#### *3.2. Fluxes Estimated with BP Solution to Heat Transport Equation*

The calculation of fluxes using the BP solution was merely possible for 18 of the 30 points (points 1 to 18). In some locations (points 19 and 21), the temperature of the lagoon bed showed a small decrease for the first 15 cm and then dropped more than 4oC over 5 cm (from the sensor located at 15 cm depth to the sensor located at 20 cm depth) (Figure 6). In other locations (points 20 and 22 to 30), the temperature stabilized at 16–17 ◦C at 35-40 cm depth (Figure 6), far from the 9.8 ◦C expected for groundwater in the study area (Figure 5b,c). In all these locations (19–30), it was not possible to satisfactorily fit the BP to the data. It was noted that the vertical temperature data not fitted by the BP was located farthest offshore.

**Figure 6.** The temperature-depth profiles for points 5, 19 and 20. An example of reasonably fitted data is point 5 (empty black dots), fitted with both BP (black dashed line) and SB (magenta dashed line). An anomalous change in the temperature can be observed for point 19 (red dots) from depth 0.15 m to depth 0.2 m. The temperature at point 20 (blue dots) remained above 16 ◦C along the entire 0.5 m profile.

For the scenario of homogeneous sediment, the exchange fluxes ranged from 0.86 to 6.05 cm day−<sup>1</sup> with a mean of 3 cm day−<sup>1</sup> and a standard deviation of 1.56 cm day<sup>−</sup>1. A general trend was observed of higher exchange fluxes close to the shore and decreased with the distance offshore (Figure 5b). Thus, these data indicate that most of the groundwater discharge took place within the first meters of the lagoon.

For the scenario of heterogeneous sediment, the exchange fluxes ranged from 0.69 to 2.25 cm day−<sup>1</sup> with a mean of 1.44 cm day−<sup>1</sup> and a standard deviation of 0.47 cm day−1. In this case, the higher discharge values were also located close to the shore, but a second zone with high values can be observed at points 12 to 17 (Figure 5c).

#### *3.3. Fluxes Estimated with SB Solution to Heat Transport Equation*

The calculation of the fluxes using the SB solution was possible in 27 of the 30 points (19–21 showed anomalous temperatures, see explanation in Section 3.2). The main difference was that the

multi-layer approach allowed the authors to fit the data even if the temperature of the lowest sensor for each measurement was far from the expected groundwater temperature.

The exchange fluxes ranged from 0.09 to 34.6 cm day−<sup>1</sup> with a mean of 12.9 cm day−<sup>1</sup> and a standard deviation of 14.8 cm day<sup>−</sup>1. The high standard deviation indicates that there were extreme values. The fluxes above 8 cm day−<sup>1</sup> were obtained only for offshore points (22 to 30). At points closer to the shore (1 to 18), the maximum calculated flux was 7.78 cm day−<sup>1</sup> with a mean of 2.77 cm day<sup>−</sup>1. There were no obvious patterns in the distribution of the fluxes (Figure 5d), but two different high-discharge zones were distinguished: The first zone located 18 m offshore (points 12 to 18) and the second zone located from 30 to 43 m offshore (points 22 to 30).

#### *3.4. Fresh Groundwater and Lagoon Water Proportion in Exchange Fluxes*

The proportion of fresh groundwater in the seepage meter samples showed that after point 18, there was little to no fresh water in the exchange flux (Figure 7). Further, the highest percentage of groundwater collected was approximately 78%, indicating that even the points located close to the shore were influenced by brackish lagoon water.

**Figure 7.** The mean percentage of fresh groundwater calculated from the seepage meters experiment. The red dashed line indicates the approximate location of the brackish lagoon water–fresh groundwater interface.

#### *3.5. Method Comparison*

The comparison between the two analytical solutions and seepage meter results in points 1 to 18 are analyzed due to lack of data for all methods at the remaining measuring points. This showed that fluxes obtained with both BP and SB for points closer to shore (1 to 8) were generally much higher than those obtained with the seepage meters (Figure 8). In the remaining points, the fluxes obtained with the seepage meters were generally higher than those obtained with the temperature approach.

The mean difference between the fluxes calculated with BP and the fluxes measured with seepage meters was 3.55 cm day−<sup>1</sup> with a standard deviation of 1.96 cm day−<sup>1</sup> for the scenario with homogeneous soil and 2.66 cm day−<sup>1</sup> with a standard deviation of 2.33 cm day−<sup>1</sup> for the scenario with heterogeneous soil. The mean difference between the fluxes calculated with SB and the fluxes measured with seepage meters was 2.14 cm day−<sup>1</sup> with a standard deviation of 1.56 cm day<sup>−</sup>1.

**Figure 8.** The mean of exchange fluxes measured with seepage meters and calculated using analytical solutions for temperature profiles.

#### **4. Discussion**

#### *4.1. Comparison of Seepage Meter and Thermal Methods*

The results show the existence of an exchange of water between the subsurface and the lagoon. The measurements with seepage meters indicated that this exchange was generally from the subsurface to the lagoon (groundwater discharge) with only two locations showing groundwater recharge from the lagoon. However, seepage at those locations were negligible as the measured fluxes (−0.02 cm day<sup>−</sup>1) fell below the detection limit [11]. Only the discharge fluxes were inferred using the two analytical solutions for the heat transport in the lagoon bed. It is important to note that in some cases, it was not possible to reproduce the temperature-depth data with the analytical solutions.

The seepage meter results show a spatial distribution with two main zones of groundwater discharge. A similar behaviour was observed by Duque et al. [20] and corresponded to the position of the saline wedge showing a first peak in the discharge 18 m offshore and to a wave-induced pumping flux yielding a second peak, 30–43 m offshore [24]. A similar phenomenon was also observed by Shinn et al. [30]. The first main discharge zone matched the end of significant fresh groundwater component in the exchanged water from the seepage meters experiment (around point 18).

The BP method was able to reproduce the temperature profiles for the first 18 points. Those 18 locations corresponded to the area in which a portion of fresh groundwater was measured in the seepage meter samples. The fluxes at these locations decreased as the distance to the shore increased (Figure 5c). This agrees with the commonly assumed decreasing trend of groundwater discharge with the distance offshore in the surface water bodies such as lakes [31]. This indicates that the BP method is not applicable if the origin of the discharged groundwater is lagoon water.

The SB method was able to reproduce most of the temperature profiles and the calculated fluxes showed a pattern similar to the ones obtained from the seepage meters (Figure 5d). In this case, the method was applicable also when there was no fresh water in the exchanged water. Despite this, the fluxes calculated for the offshore points (22 to 30) were very high in comparison with the results obtained by other methods used at other points, and also in comparison with previous studies in the study area [7,20]. This may indicate that, although the method can be applied and the temperature-depth profiles are correctly fitted with the SB analytical solution, the results obtained are not realistic and should be considered with care.

The analytical solutions to reproduce the temperature-depth data suggested a higher discharging flux than the seepage meters for the first 6–8 measuring points. For the remaining points, the fluxes measured with the seepage meters were higher. This may be caused by the proximity of the saline wedge around points 10–12. On the one hand, in the fresher part of the saline wedge, the main component of the exchanged flux is fresh groundwater and it can be well captured by the thermal methods. On the other hand, closer to the saline wedge the proportion of the recirculated lagoon water increased (Figure 3), altering the thermal signal in the groundwater. Thus, the reliability of the thermal methods decreases, because it is difficult to precisely define an appropriate lower boundary condition. A common temperature of 9.8 ◦C was used at all points at a depth of 5 m, but this may not necessarily be representative. Although the groundwater discharge fluxes obtained from the different methods were in the same order of magnitude, it is clear that both the location and the range of fluxes differ considerably depending on the methodology. This has been reported in other studies (e.g., [32]) and highlights the importance of complementing methods for obtaining groundwater exchange fluxes with other hydrological measurements.

When analysing the proportion of fresh groundwater and lagoon water in the seepage meter samples, it was observed that many samples did not have any fresh groundwater. This indicates that the water flowing to the lagoon at these points is recirculated brackish lagoon water (Figure 3).

#### *4.2. Direct Measurements (Seepage Meters): Pros and Cons*

The use of seepage meters for obtaining the groundwater discharge fluxes to Ringkøbing fjord is the most straightforward methodology used in this study. The fact that the flux is obtained by a simple mass balance makes it a robust and reliable method. Furthermore, as the exchanged water is stored in the collecting plastic bags, a geochemical analysis can be performed, complementing the dataset with insights regarding the origin of the discharged groundwater. Care must be taken, however, when interpreting the water quality of seepage water as the seepage meter can by itself change the oxidation state of the lagoon bed [33]. Performing several rounds of measurements, as described in this work, allows for a better estimation of the robustness of the method. The measured fluxes had an average standard deviation of less than 1 cm day−<sup>1</sup> for all six rounds, which indicates that the method is reliable.

However, this method has some noticeable drawbacks that may make it not applicable in some situations. An extended description of such situations can be found in [34]. In this particular case, the most important ones were: First, the installation of the seepage meters has to be performed in advance before the experiment takes place [8]. In the case of geochemical sampling of the exchanged groundwater, the amount of time needed for the seepage meters to replace the trapped surface water with discharged groundwater (if any) is inversely proportional to the exchange flux, which is unknown before the experiment takes place. This implies that the discharge flux needs to be assessed prior to the experiment based on previous studies or by the use of other methods (for example, the thermal method), thus increasing the uncertainty of the results. In this case, the equipment was installed approximately one week before the data collection. Once the measurements were performed and the exchange fluxes known, this study could confirm that the time from the installation until the experiment was long enough, but it could not have been the case. Second, the amount of human resources needed for an intense sampling campaign, as the one described here, is extensive [35]. Decreasing the number of measuring points is not always a good option, as it has been observed that to obtain reliable results, many seepage meters are needed [36,37]. Third, the methodology has some associated uncertainty. Shinn et al. [30] observed that the volumetric data obtained from seepage meters in areas exposed to currents, waves and swells, like the studied area and almost any costal area, must be used with caution as the induced flux from wave-induced pumping affects the measurements. Further, Lee [8] highlighted the potential issues in the measurements due to the improper placement of the cylinder. The seepage meter installation is a manual activity and it may not be possible to ensure that all requirements for reliable data acquisition are met. Finally, the use of correction coefficients to better estimate the exchange flux has been studied for many years (e.g., [8,38]) and it may be needed. This study did not apply any of the correction coefficients from the literature. This parameter seems to be rather equipment dependent and for the comparison of methods carried out in this work, the authors considered it appropriate to use the original value obtained from the mass balance.

#### *4.3. Analytical Solutions for Temperature-Depth Profiles (BP and SB): Pros and Cons*

The temperature-depth data collection was performed by one person during one of the field campaign days. The two temperature probes and two logger boxes were used to obtain the streambed temperature profiles. Each measurement takes ten minutes which implies that in approximately three hours, all measurements were done. Both the temperature probes and the logger boxes can be brought to the field in a small van and can be carried by a single person. This is a noticeably smaller amount of time and human resources needed in comparison with the seepage meter measurements. Therefore, the field campaigns are cheaper and more measurements can be done in the same amount of time.

The use of heat as a tracer, and specifically analytical solutions to study temperature profiles, has several limitations that need to be taken into account.

Both the BP and SB rely on the assumption of vertical and steady water and heat flow [26,28]. These conditions are rarely found outside laboratory experiments and are unknown a priori, unless there is some knowledge of the specific site. In this case, a location near the saline wedge is expected to have exchange fluxes (if any) with a strong vertical component due to the density driven flow. However, the locations influenced by a recirculating flux may be biased by a horizontal or transverse flux, which may be responsible for a weaker performance of the analytical methods within our study. Furthermore, the data collection took place in mid-July assuring that the temperature fluctuation of the surface is as small as possible (higher fluctuations with changing weather conditions are expected in spring and autumn). Thus, the assumption of steady-state is more likely to be appropriate.

The BP assumes that the media in which the temperature is measured is homogeneous, but this is the exception rather than the rule in natural systems and previous studies indicated that this was not the case in this study site. Furthermore, several parameters need to be known a priori: The average steady groundwater temperature; the depth at which this temperature is reached; the thermal conductivity of the medium. Here, several temperature-depth profiles were measured in boreholes around the study area, thus allowing the assessment of the groundwater temperature regime (9.8 ◦C at 5 m). The subsurface is on average 1.8 ◦C warmer in the saline sector compared with the fresh groundwater, similarly to what was observed by Befus et al. [18].

The SB allows for the definition of a multi-layered system. This study differentiated three different zones in this study area, providing the basis for studying the effect of the horizontal and vertical heterogeneity in the effective bulk thermal conductivity. The bulk thermal conductivity values [20] were used as a framework to consider a more complex hydrogeological system. Fitting the data with the SB solution was considerably easier than with the BP due to the flexibility on the lower boundary condition (lower sensor temperature versus steady groundwater temperature). This decreases the amount of uncertainty derived from the assumption of a steady groundwater temperature at unknown depth. At the same time, it also forces the solution to satisfy this condition, which does not allow taking into account the bigger scale thermal behaviour of the system (i.e., the thermal regime at greater depths). With this solution, the estimated fluxes deviated substantially from the results of previous studies [7,24]. Fixing the lower boundary condition to the temperature of the sensor located at a depth of 50 cm (instead of using a deep groundwater steady temperature as in the BP method) focuses the calculation of the flux to only the upper part of the lagoon bed, disregarding the rest of the subsurface. This may lead to under/overestimations of the flux if the processes occurring in that upper section are not representative of the overall behaviour of the system. Also, the fact that these fluxes were obtained for the more offshore points indicates that the thermal methods usually applied to study GW-SW

interactions in fresh water environments seem less reliable, when there is a strong component of saline water recirculation in the flux exchange. Therefore, obtaining unrealistic results in some locations and the use of more parameters that may be unknown (number of layers and thickness, thermal conductivity for each layer) are important limitations to take into account before using this method.

Some of the temperature profiles were very uniform (temperature not decreasing much with depth). This could actually indicate that there was a dominant vertical downward flow in the top 50 cm (contrary to other measurements) or that the flow was not strictly vertical upwards having a significant horizontal or transversal component. The complex flow field near a salt water wedge may be the explanation, so the use of 1D solutions is then not warranted or, at least, it is only the vertical component that is estimated in many cases. The comparison of these data and the seepage meter results helps to distinguish in which cases this may occur. That is, the differences between the fluxes estimated by the thermal and seepage meter methods could be explained by the fact that the seepage meter, located right at the lagoon bed, is a composite measure of the discharging flux, whereas, it is only the average vertical flux component that is measured over 50 cm by the thermal method. However, the fact that this phenomenon (nearly uniform temperature profiles) is observed at the more offshore points seems to confirm that the applicability of the thermal methods is not well suited for these environments.

#### **5. Conclusions**

The use of different methodologies to study the groundwater discharge to a lagoon with saltwater intrusion and distinct thermal regimes reveal that there are several challenges to overcome in order to obtain reliable estimates of submarine groundwater discharge. It points out the difficulties in choosing which tools to apply in the field, but it also gives insights to what method to use depending on which information of the study site is available beforehand.

The implementation of the two analytical solutions was only satisfactory for 18 points located closer to the shore and with a significant fresh groundwater component in the exchange flux. This has two implications: (1) Though the seepage meters allowed the measurement of the amount of exchanged water through the lagoon bed, it is not possible to infer the origin of this water unless other measurements are performed (e.g., EC); (2) it indicates that the thermal methods may not be well suited to study groundwater–surface water interactions if surface water recirculation processes are of importance due to the subsurface thermal regime (essentially providing unknown temperature boundary conditions).

Characterizing the saltwater-freshwater interaction in coastal areas is a key issue in the application of temperature measurements to estimate fluxes. The establishment of the areas dominated by saltwater or freshwater can be indicators of the limits of the applicability for the usual assumption of stable groundwater temperature. Also, the detection of the saltwater-freshwater interface is highly relevant due to the dominance of vertical-horizontal fluxes in the study area.

Choosing between the Bredehoeft and Papadopulos (BP) [26] or the Shan and Bodvarson (SB) [28] analytical solution to evaluate the subsurface temperature profiles depends mainly on the knowledge of the site. On the one hand, if groundwater is monitored in the surrounding area and thermal data is available, SB can be applied, allowing for a better constraint of the lower boundary condition. Also, if stratification of the subsurface is known or expected and it is possible to measure or infer the different thermal parameters for each layer, SB provides a good solution to better fit the thermal profiles. BP is better suited for environments with few or no prior data acquisition, as only the average groundwater temperature in the area and one value for the thermal conductivity are needed. SB reproduces more complex systems that are frequent in nature and BP only homogeneous scenarios. In both cases, it is necessary to address the implications that surface water recirculation processes have on the calculated exchange fluxes.

Nevertheless, the seepage meter and thermal methods have their own strengths and based on our study, it can be concluded that the thermal methods represent a fast way of getting the right order of magnitude of the exchange fluxes. The results (data/model) can be generated within one day and can be used to design an efficient seepage meter campaign (where and how much), which would likely give the most representative estimate of the exchange fluxes. Additionally, seepage meter measurements may allow observing the origin of the flux when taking samples for water chemistry.

**Author Contributions:** Conceptualization, J.T.-C., C.D and P.E.; methodology, J.T.-C., C.D. and P.E.; validation, J.T.-C., C.D. and P.E.; formal analysis, J.T.-C. and C.D.; investigation, J.T.-C., C.D., S.K. and P.E.; resources, P.E., C.D. and S.M.; writing—original draft preparation, J.T.-C.; writing—review and editing, J.T.-C., C.D., S.K., S.M. and P.E.; visualization, J.T.-C.; supervision, P.E. and C.D.; project administration, C.D. and P.E.; funding acquisition, P.E. and C.D.

**Funding:** This research has received funding from the European Union's Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie Grant Agreement No 722028 and the European Union Seventh Framework Programme FP7/2007–2013 under the Grant Agreement No 624496.

**Acknowledgments:** Comments from the editor Jörg Lewandoswki and from the two reviewers are gratefully acknowledged. We thank Iris Tobelaim and Christian Hansen for their help in the field.

**Conflicts of Interest:** The authors declare no conflicts 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/).
