# Nutrient fluxes associated with submarine groundwater discharge from karstic coastal aquifers (Côte Bleue, French Mediterranean Coastline)

Simon Bejannin, Joseph James Tamborski, Pieter van Beek, Marc Souhaut, Thomas Stieglitz, Olivier Radakovitch, Christelle Claude, Pascal Conan, Mireille Pujo-Pay, Olivier Crispi, et al.

► To cite this version:

Simon Bejannin, Joseph James Tamborski, Pieter van Beek, Marc Souhaut, Thomas Stieglitz, et al. Nutrient fluxes associated with submarine groundwater discharge from karstic coastal aquifers (Côte Bleue, French Mediterranean Coastline). Frontiers in Environmental Science, 2020, 7, 10.3389/fenvs.2019.00205. ird-02483073

HAL Id: ird-02483073

https://ird.hal.science/ird-02483073

Submitted on 18 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

# Nutrient Fluxes Associated With Submarine Groundwater Discharge From Karstic Coastal Aquifers (Côte Bleue, French Mediterranean Coastline)

Simon Bejannin¹, Joseph James Tamborski¹*†, Pieter van Beek¹*, Marc Souhaut¹, Thomas Stieglitz², Olivier Radakovitch²,³, Christelle Claude², Pascal Conan⁴, Mireille Pujo-Pay⁴, Olivier Crispi⁴, Emilie Le Roy¹† and Claude Estournel¹

¹LEGOS, Laboratoire d'Etudes en Géophysique et Océanographie Spatiales UMR 5566 (Université de Toulouse, CNES, CNRS, IRD, UPS), Observatoire Midi Pyrénées, Toulouse, France, ²Aix-Marseille Univ, CNRS, IRD, INRA, Coll France, CEREGE, Aix-en-Provence, France, ³IRSN (Institut de Radioprotection et de Sûreté Nucléaire), PSE-ENV/SRTE/LRTA, Saint-Paul-Lez-Durance, France, ⁴Laboratoire d'Océanographie Biologique, Sorbonne Université, UPMC, CNRS, UMR 7621, Laboratoire Arago, Banyuls-sur-Mer, France

Determination of submarine groundwater discharge (SGD) from karstic coastal aquifers is important to constrain hydrological and biogeochemical cycles. However, SGD quantification using commonly employed geochemical methods can be difficult to constrain under the presence of large riverine inputs, and is further complicated by the determination of the karstic groundwater endmember. Here, we investigated a coastal region where groundwater discharges from a karstic aquifer system using airborne thermal infrared mapping and geochemical sampling conducted along offshore transects. We report radium data ($^{223}$Ra, $^{224}$Ra, $^{228}$Ra) that we used to derive fluxes (water, nutrients) associated with terrestrial groundwater discharge and/or seawater circulation through the nearshore permeable sediments and coastal aquifer. Field work was conducted at different periods of the year to study the temporal variability of the chemical fluxes. Offshore transects of $^{223}$Ra and $^{224}$Ra were used to derive horizontal eddy diffusivity coefficients that were subsequently combined with surface water nutrient gradients ($NO_2^-$ + $NO_3^-$, DSi) to determine the net nutrient fluxes from SGD. The estimated DSi and ($NO_2^-$ + $NO_3^-$) fluxes are 6.2 ± 5.0 *10³ and 4.0 ± 2.0 *10³ mol d⁻¹ per km of coastline, respectively. We attempted to further constrain these SGD fluxes by combining horizontal eddy diffusivity and $^{228}$Ra gradients. However, SGD endmember selection in this area (terrestrial groundwater discharge vs. porewater) adds further uncertainty to the flux calculation and thus prevented us from calculating a reliable flux using this latter method. Additionally, the relatively long half-life of $^{228}$Ra (5.75 y) makes it sensitive to specific circulation patterns in this coastal region, including sporadic intrusions of Rhône river waters that impact both the $^{228}$Ra and nutrient surface water distributions. We show that SGD nutrient fluxes locally reach up to 20 times the nutrient fluxes from a small river (Huveaune River). On a regional

OPEN ACCESS

**Edited by:**

Henrietta Dulai,
University of Hawaii at Manoa,
United States

**Reviewed by:**

Aaron Beck,
Geomar Helmholtz Center for Ocean
Research Kiel, Germany

Kimberly Null,
Moss Landing Marine Laboratories,
United States

*Correspondence:

Joseph James Tamborski
jtamborski@whoi.edu
Pieter van Beek
vanbeek@legosobs-mip.fr

†Present address:

Joseph James Tamborski and
Emilie Le Roy,
Department of Marine Chemistry and
Geochemistry, Woods Hole
Oceanographic Institution, Woods
Hole, MA, United States

**Specialty section:**

This article was submitted to
Water and Wastewater Management,
a section of the journal
Frontiers in Environmental Science

Received: 31 May 2019

Accepted: 20 December 2019

Published: 18 February 2020

**Citation:**

Bejannin S, Tamborski JJ, van Beek P,
Souhaut M, Stieglitz T, Radakovitch O,
Claude C, Conan P. Pujo-Pay M,
Crispi O, Le Roy E and Estournel C
(2020) *Nutrient Fluxes Associated
With Submarine Groundwater
Discharge From Karstic Coastal
Aquifers* (Côte Bleue, French
Mediterranean Coastline).
Front. Environ. Sci. 7:205.
doi: 10.3389/fenvs.2019.00205
scale, DSi fluxes driven by SGD vary between 0.1 and 1.4% of the DSi inputs of the Rhône River, while the (NO₂⁻ + NO₃⁻) fluxes driven by SGD on this 22 km long coastline are between 0.1 and 0.3% of the Rhône River inputs, the largest river that discharges into the Mediterranean Sea. Interestingly, the nutrient fluxes reported here are similar in magnitude compared with the fluxes quantified along the sandy beach of La Franqui, in the western Gulf of Lions (Tamborski et al., 2018), despite the different lithology of the two areas (karst systems vs. unconsolidated sediment).

**Keywords:** submarine groundwater discharge, Mediterranean sea, radium isotopes, thermal infrared remote sensing, nutrient fluxes, GEOTRACES

INTRODUCTION

Submarine Groundwater Discharge (SGD) is now recognized
as a vector for many chemical elements that impact the
biogeochemistry and ecology of the coastal seas (Moore, 1996;
Charette and Buesseler, 2004; Slomp and Van Cappellen,
2004; Burnett et al., 2006). SGD can facilitate phytoplankton
development in coastal areas (Paytan et al., 2006) or may play
a role in eutrophication of coastal seas, including in some cases
harmful algal blooms (Gobler and Sañudo-Wilhelmy, 2001). SGD
includes (i) the discharge of terrestrial groundwater to the sea and
(ii) the circulation of seawater through the coastal aquifer and
permeable coastal sediments (Burnett et al., 2003; Moore, 2010).
The mixing zone between groundwater and seawater is defined as
the subterranean estuary (Moore, 1999), a geochemically reactive
subsurface zone where many chemical species are modified by
biogeochemical reactions and water-rock interactions. Although,
SGD has been shown to supply essential nutrients (Charette et al.,
2001; Slomp and Van Cappellen, 2004; Hwang et al., 2005; Knee
et al., 2010; Beusen et al., 2013; Tamborski et al., 2017) and trace
elements to the coastal sea (Jeong et al., 2012; Trezzi et al., 2017),
relatively few studies have been conducted on nutrient fluxes
driven by SGD in the Mediterranean Sea (Garcia-Solsona et al.,
2010; Weinstein et al., 2011; Tovar-Sánchez et al., 2014; Rodellas
et al., 2015, 2017; Tamborski et al., 2018). In particular, little
information is available on the nutrient fluxes associated with
SGD along the French Mediterranean coastline (Rodellas et al.,
2017; Tamborski et al., 2018), despite the presence of several
well-known karstic springs (Arfib et al., 2006; Fleury et al., 2007;
Stieglitz et al., 2013; Bejannin et al., 2017).

Karstic springs may contribute groundwater to the coastal
ocean from different geological sources over different temporal
and spatial scales; thus, evaluation of SGD from coastal karst
aquifers remains a significant challenge (Montiel et al., 2018).
Karst aquifers may respond rapidly to even minor precipitation
events, and this rapid flushing can result in significant export
events (Fleury et al., 2007). Certain karstic springs may
continuously flow in the absence of rainfall, as is the case of
the brackish Port-Miou spring, one of the largest groundwater
springs that enter the Mediterranean Sea (Arfib and Charlier,
2016). Evaluation of nutrient loads via karst SGD is therefore
critical in semi-arid, oligotrophic regions like the Mediterranean
Sea (Tovar-Sánchez et al., 2014). Ra isotopes have proved to be
useful tracers of SGD from coastal karst aquifers, particularly in

the Mediterranean Sea (Ollivier et al., 2008; Garcia-Solsona et al.,
2010; Mejías et al., 2012; Tovar-Sánchez et al., 2014; Baudron
et al., 2015; Rodellas et al., 2015; Montiel et al., 2018; Tamborski
et al., 2018).

A recent study in the Mediterranean Sea reported SGD-
driven nutrient estimates based on a basin-wide ²²⁸Ra mass
balance (Rodellas et al., 2015). The ²²⁸Ra inventory, after
considering known ²²⁸Ra sources and ²²⁸Ra sinks, is divided
by the ²²⁸Ra submarine groundwater endmember activity to
determine a SGD flux. Importantly, Rodellas et al. (2015)
note that ²²⁸Ra SGD endmembers reported throughout the
Mediterranean Basin span a wide-range, from 10² to 10⁵ dpm
m⁻³. This example illustrates the complexity in determining
the submarine groundwater endmember, particularly in karst
environments. It is thus important to also conduct studies at
regional- and local-scales (i) to better constrain the processes
involved in the transfer of chemical species between land and the
coastal seas and (ii) to elucidate spatial and temporal variability
associated with these fluxes.

In this study, we investigated the karstic coastline of the Côte
Bleue region located west of the city of Marseille, France. This
coastline is located east of the Rhône river plume, where riverine
nutrient inputs support 23–69% of the primary production
of the Gulf of Lions (Ludwig et al., 2009). This coastline is
of specific interest since it may exhibit significant transfer of
water and associated solutes toward the sea through terrestrial
groundwater discharge, as is the case in other karstic systems
worldwide (Garcia-Solsona et al., 2010; Pavlidou et al., 2014;
Tovar-Sánchez et al., 2014). SGD may constitute an additional
input of nutrients that contributes to sustain the primary
production in the coastal sea. Here, we investigated the entire
Côte Bleue coastline (22 km), in order to evaluate if coastal
sections without springs could also contribute as a nutrient
source. In addition, we focused on the temporal variability of
these fluxes. We thus conducted offshore transects of radium
isotopes that were used to derive horizontal eddy diffusivity
coefficients (K<sub>h</sub>), that were subsequently combined with surface
water nutrient gradients (DSi, NO<sub>2</sub><sup>-</sup> + NO<sub>3</sub><sup>-</sup>) in order to
determine the net nutrient flux from SGD. Offshore transects
of ²²⁸Ra were also investigated with the aim to derive SGD-
driven nutrient fluxes by combining K<sub>h</sub> and ²²⁸Ra gradients
with the SGD endmember. Fluxes were investigated over a
1-year period (April 2016, October 2016, December 2016,
and March 2017).
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | Location of the study site, Côte Bleue, a karstic coastline that is situated in south-eastern France. Previous hydrogeological characterizations (based on Carte Hydrogéologique des Bouches du Rhône, 1972, BRGM) are reported. The location of the offshore transects (T1 to T7) are shown, with diamonds indicating the surface water samples. Arrows indicate the presence of known springs along this coastline.</figcaption></figure>

# MATERIALS AND METHODS

## Study Site

Côte Bleue is a 22 km long stretch of karstic coastline situated along the French Mediterranean Sea, within the Gulf of Lions. It is located 15 km east of the Rhône River, the largest river discharging into the French Mediterranean Sea, and is 15 km west of the Huveaune River, a small river discharging in Marseille (**Figure 1**). The steep rocky coast hosts a series of coves (Calanques), which shelter small harbors and beaches. Submarine flow paths and springs are known to exist from its three watersheds (**Figure 1**; Carte Hydrogéologique des Bouches du Rhône, 1972, BRGM). Little quantitative information exists on groundwater discharge magnitude, its seasonality and associated solute fluxes. This area receives little to no rain during the summer (from June to September) with low precipitation rates throughout the year (**Figure S1**). The total precipitation recorded during the 2 years of this study (2016–2017) is 632 mm (meteociel.fr). The water depth in this area is 6–20 m at 100 m offshore and 70–80 m depth at 8 km offshore (**Figure 1**). There are no permanent riverine inputs along Côte Bleue. The land and sea adjacent to the coast are protected areas (Natura 2000) and are subject to biologic monitoring (Jouvenel et al., 2004); fishing and diving is banned in the area. Artificial reefs have been installed to improve fish populations (Jensen et al., 2012).

The hydrodynamical regime of the Gulf of Lions is very complex (Millot, 1990). The main general circulation feature that influences the Gulf of Lions is the Northern Current (NC) that flows along the continental slope from east to west (Millot, 1990; see e.g., Figure 1 of Gatti et al., 2006). The NC was shown

to occasionally intrude on the shelf (Petrenko, 2003). Maps of horizontal currents in the eastern part of the Gulf of Lions indicate a predominant westward/northwestward flow offshore Côte Bleue (Pairaud et al., 2011). However, these patterns were found further offshore compared to the coastal region that was investigated in the present study. No information is thus available on the circulation patterns in the first 8 km along the coastline. West of Côte Bleue, the circulation pattern generally drives the Rhône river plume westward and therefore the Rhône river plume does not impact the investigated region. In some rare cases, however, the shelf waters in the area of Marseille were shown to be influenced by Rhône river intrusion that may impact the biogeochemical patterns in the area (Gatti et al., 2006; Fraysse et al., 2014).

## Field Methods

For reconnaissance of SGD locations, airborne thermal infrared (temperature) and in-water radon mapping were carried out. Airborne thermal infrared images were acquired on 20 September 2012 using a FLIR Systems ThermaCAM SC 3000 along the coast. Detailed information on the flight acquisition is presented in Bejannin et al. (2017). Groundwater temperature tends to follow the average annual air temperature; thus, surface water temperature may be qualitatively used to identify coastal areas potentially impacted by groundwater discharge. Images were acquired in the morning (7:30 a.m.) to provide the largest relative temperature difference between suspected discharging groundwater and ambient surface waters. Thermal infrared images were mosaicked and manually georeferenced against satellite imagery (Google Earth) for visualization.
<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | Thermal infrared imagery acquired along Côte Bleue on 20/09/2012 at 7:40 a.m. in the insets that show different zooms of the coastline (a–e), the sampling locations of the nearshore water and porewater samples are displayed on the TIR images. Black diamonds denote nearshore water samples and white diamonds denote porewater samples. It is important to note that the water sample collection and the thermal infrared overflight were not conducted at the same time.</figcaption></figure>

Radon mapping of surface waters was measured *in situ* on 27–28 October 2016 using two Rn-in-air detectors (RAD7-Durridge, Co. Inc.) routed simultaneously through an air-water exchanger in order to locate potential sites of terrestrial groundwater discharge (Burnett and Dulaiova, 2003; Dulaiova et al., 2005; Stieglitz et al., 2013; Cockenpot et al., 2015). Surface water was pumped from 0.5 m depth at a constant flow rate of 2.5 L min⁻¹ and filtered through an 80 µm cartridge while continuously moving along the coastline. Data integration of the run was fixed to 15 min, and the system was equilibrated at least for 30 min before the start. We considered that one integrated measurement corresponds to the activity of the water collected 15 min earlier. The distance between each point ranged between 200 and 400 m. The western transect was done from east to west (day 1), and the eastern one from west to east (day 2). In both cases, this direction affected the measurement due to the delay necessary for the equilibration time from high to low activity (Stieglitz et al., 2010). All activities were corrected from temperature, humidity (using the Durridge software) and salinity.

Data from the airborne thermal infrared imagery and continuous ²²²Rn survey were used to select nearshore and offshore transect sampling stations (Figures 1, 2). Nearshore surface water samples (~20 L; within 5 m of the shoreline) and porewater samples (~2–5 L) for Ra isotopes were collected using either a submersible pump or a manual hand pump. These samples were collected in the bay heading transect 1 (Laurons Bay) and in the areas where thermal infrared temperature anomalies exist (Figure 2). Three brackish karstic springs were sampled in May 2017 in a small bay (Laurons bay) at transect 1 by scuba divers. Additional porewater samples were collected

from hand-dug bore holes which intersected the saturated zone of the beach for areas where there was permeable sediment (Figure 2). Seawater samples (~100 L; from 100 m to 8.5 km offshore) were collected along offshore transects (Figure 1) using a submersible pump on board the R/V Antedon II. Various transects and coastal water samples were collected on April 28, October 26–27, December 7–8 in 2016, and on March 20–22 in 2017. Seven shore-perpendicular transects were investigated over these four campaigns (Figure 1); only transect T4 was sampled during all four seasonal campaigns (chiefly due to weather constraints). This includes transects conducted offshore sites of known or potential karstic springs and sites where no groundwater discharge exists. Salinity was recorded *in situ* during the sampling operations using the shipboard CTD for the seawater samples, and a handheld WTWTM probe (Xylem) for the nearshore samples. Water samples for nutrients were collected during each cruise using a submersible pump and filtered in the field (0.45 µm).

## Laboratory Methods

Water samples for Ra isotopes were passed through MnO₂-coated acrylic fibers (“Mn-fiber”) at a flow rate <1 L min⁻¹ to ensure a 100% yield of Ra extraction onto the fiber (Moore and Reid, 1973). Once back in the laboratory, the fibers were rinsed three times with Ra-free deionized water and dried with compressed air until a water:fiber ratio of 1:1 was reached (Sun and Torgersen, 1998). Short-lived ²²³Ra and ²²⁴Ra activities were determined using a Radium Delayed Coincidence Counter (RaDeCC; Moore and Arnold, 1996). Detector efficiencies and error propagation calculations were calculated after Moore
(2008) and Garcia-Solsona et al. (2008). A first counting session was run after sample collection to determine the total $^{224}$Ra and $^{223}$Ra activities. For high activity samples, another counting session was run after 1 week to determine the total $^{223}$Ra activity. The Mn-fibers were analyzed again 3 weeks after sampling to determine the $^{224}$Ra activities supported by $^{228}$Th. These supported activities were then subtracted to the total $^{224}$Ra activities to determine excess $^{224}$Ra (denoted $^{224}$Raex). Select samples were counted at least 1 year after sampling to determine $^{228}$Ra activities using RaDeCC via $^{228}$Th ingrowth, following Moore (2008). Select samples (T1 and T7 for all cruises) were additionally analyzed for $^{228}$Ra using the low-background gamma detectors at the LAFARA underground laboratory in Ferrières, French Pyrénées (van Beek et al., 2010, 2013). Prior to gamma analysis, the Mn-fibers were dried, pressed using a hydraulic press at 20 metric tons, placed into plastic boxes and vacuum sealed into bags to prevent from any Rn loss. A semi-planar detector (ORTEC/AMETEK; van Beek et al., 2013) was used to determine $^{228}$Ra activities from an average of the $^{228}$Ac photo-peaks (338, 911, and 969 keV).

Nutrient samples were filtered at 0.45 µm, immediately poisoned with HgCl2 (10 µg L⁻¹) and stored at 4°C in the dark. In the laboratory, nitrate (NO₃⁻), nitrite (NO₂⁻), phosphate (PO₄³⁻) and dissolved silica (DSi) concentrations were measured on a continuous flow autoanalyzer Technicon® AutoAnalyzer II at LOMIC, Banyuls-sur-Mer (Aminot and Kérouel, 2004). The analytical precision of NO₃⁻, NO₂⁻, PO₄³⁻, and DSi is ±0.02, ±0.01, ±0.02, and ±0.05 µM, respectively.

## Water and Solute Flux Estimations

Several methods based on the use of geochemical tracers (i.e., radium isotopes and radon) are used to quantify SGD fluxes (Burnett et al., 2006), including isotope mass balances and gradient flux calculations. The radium quartet ($^{223}$Ra, t₁/₂ = 11.4 d; $^{224}$Ra, t₁/₂ = 3.66 d; $^{226}$Ra, t₁/₂ = 1,600 y; $^{228}$Ra, t₁/₂ = 5.75 y), as well as $^{222}$Rn (t₁/₂ = 3.83 d), have been widely used to study SGD worldwide (Moore, 1996; Charette et al., 2001; Paytan et al., 2006; Rodellas et al., 2015; Tamborski et al., 2015). These radionuclides are produced within an aquifer by the decay of their sediment-bound U/Th series parent nuclide. Production near the mineral surface provides sufficient energy to recoil the daughter isotope into the surrounding pore fluid (Swarzenski, 2007). Ra isotopes tend to be adsorbed onto sediments at low ionic strengths (i.e., freshwater); however, Ra isotopes are desorbed and released into the surrounding pore fluid under brackish conditions (Webster et al., 1995). Ra isotope activities are typically 2–3 orders of magnitude greater in brackish groundwaters than in surface waters; thus, Ra isotopes are powerful tracers of SGD inputs to the sea. The range of half-lives (from days to thousands of years) of these isotopes allows for the quantification of SGD flow-paths which may occur over a wide-range of time scales (Moore, 1996; Charette et al., 2001).

SGD fluxes of conservative chemical elements can be quantified by combining horizontal eddy diffusivity coefficients $K_h$ (derived from short-lived Ra isotopes) and offshore gradients of the elements under consideration (Moore, 2000a; Windom et al., 2006). This method provides a direct estimate of the

chemical flux, and the knowledge of the SGD endmember is not required, which is an advantage since the endmember is not always easy to determine. However, the calculation is only valid when the chemical element under consideration behaves conservatively in seawater, that is, it is not significantly impacted by biological or geochemical processes within the time-frame of the coastal water residence time. Further, this assumes that the chemical element gradient is derived solely from SGD. Alternatively, $K_h$ can be combined to the offshore gradient of $^{228}$Ra (or $^{226}$Ra) and the calculated $^{228}$Ra SGD flux can then be converted into a chemical flux by multiplying the $^{228}$Ra flux by the chemical element/$^{228}$Ra ratio in the SGD endmember. Such a calculation implies that the SGD endmember is well-constrained for $^{228}$Ra and the chemical element of interest. Both methods have been widely used to derive SGD fluxes (Charette et al., 2001; Hancock et al., 2006; Niencheski et al., 2007; Knee et al., 2011; Li and Cai, 2011).

The principal method used here to derive SGD-driven nutrient fluxes combines horizontal eddy diffusivity coefficients ($K_h$; determined by short-lived Ra surface water gradients) with surface water nutrient gradients (section Nutrient Mixing Model). As a comparison, SGD-driven nutrient fluxes were determined by combining horizontal eddy diffusivity coefficients with surface water $^{228}$Ra gradients and nutrient endmember concentrations (when it was applicable; section $^{228}$Ra Mixing Model).

## Nutrient Mixing Model

Surface water $^{223,224}$Ra activities are used to estimate the exchange rate between the Côte Bleue coastline and the coastal sea. If horizontal dispersion can be approximated by a diffusive process, then a one-dimensional model can be applied, assuming that advection is negligible and conditions are in steady-state (Moore, 2000a):

$$ \frac{dA}{dt} = K_h \frac{\partial^2 A}{\partial d^2} - \lambda A \quad (1) $$

where A is the radium isotope activity, $K_h$ is the horizontal eddy diffusivity coefficient (m² s⁻¹), d is the distance offshore and λ is the decay constant of the Ra isotope used. Assuming steady state, the solution of Equation (1) is:

$$ A_d = A_0 \exp(-d \sqrt{\frac{\lambda}{K_h}}) \quad (2) $$

where $A_d$ is the Ra isotope activity at the distance d from the coast and $A_0$ is the radium activity at the boundary (d = 0). The horizontal eddy diffusivity coefficient ($K_h$) can thus be estimated from a plot of ln($^{224}$Raex) as a function of offshore distance (Moore, 2000a). $K_h$ is then calculated from the estimate of the slope of the linear regression ($K_h = \lambda/m^2$ where m is the slope of the linear regression). Horizontal eddy diffusivity coefficients (m² s⁻¹) were calculated for each campaign. The horizontal eddy diffusivity coefficient uncertainty was determined from the error associated with the slope of the ln(Ra) vs. offshore distance relationship, which includes the analytical uncertainty of the Ra
measurement (Garcia-Solsona et al., 2008). The uncertainty of K_h is thus:

$$ \Delta(K_h) = 2K_h \frac{\Delta m}{m} \qquad (3) $$

where Δ(K_h) is the uncertainty associated with K_h, m is the slope of the linear relationship ln(Ra) vs. offshore distance and Δm is the error determined using Origin Software.

The horizontal eddy diffusivity coefficients obtained for each campaign are then multiplied by the observed surface water nutrient gradient to calculate the SGD-driven nutrient flux. The offshore nutrient gradient (µmol L⁻¹ km⁻¹) is defined as the slope of the plot of a nutrient concentration as a function of offshore distance. A nutrient flux (µmol s⁻¹ km⁻¹) is calculated by multiplying the horizontal eddy diffusivity coefficient with the nutrient gradient by the depth of the water layer impacted by SGD (i.e., terrestrial groundwater inputs and also potentially circulation of seawater through the sediment). This thickness was determined from the CTD vertical profiles acquired at each sampling station and from several vertical profiles of Ra isotopes (Table S1). A range of 5–10 m depth was used in the calculations. This calculation assumes that nutrient uptake and utilization is negligible over the time-scale of offshore water transport, and that the surface water Ra and nutrient gradients are solely derived from SGD. Uncertainty in the nutrient flux was propagated as the uncertainty of the horizontal eddy diffusivity coefficients and the uncertainty of the slope of the nutrient concentration surface water gradient. The assumptions used in the model are the following: (i) we neglect advection; (ii) the source term for Ra and nutrients is the coastline (d = 0); and (iii) nutrient consumption is assumed to be negligible. These assumptions are further discussed in section Mixing Model: K_h * Nutrient Offshore Gradient.

## $^{228}$Ra Mixing Model

SGD-driven nutrient fluxes were also calculated using the observed surface water $^{228}$Ra gradient for each transect, as a comparison to the above approach (Moore, 2000a; Charette et al., 2007). The horizontal eddy diffusivity coefficient is multiplied by the $^{228}$Ra surface water gradient and the layer depth impacted by SGD to determine a $^{228}$Ra flux (dpm d⁻¹ km⁻¹), assuming that the observed $^{228}$Ra is derived from SGD. The volumetric SGD flux is estimated by dividing the $^{228}$Ra flux by the $^{228}$Ra activity of the SGD endmember. The nutrient flux can finally be estimated by multiplying this water flux with the nutrient concentration of the endmember. We choose the average $^{228}$Ra activity and nutrient concentrations of the three brackish springs sampled in Laurons bay as the SGD endmember. The SGD endmember must be well-characterized in this approach, which can introduce additional uncertainty; this is further discussed in section $^{228}$Ra Gradient Method.

## SYMPHONIE Model

Simulations of oceanic circulation in the Gulf of Lions have been conducted using the SYMPHONIE model (Marsaleix et al., 2008, 2019). This model has been widely used to study the Rhône plume in situations of realistic forcing by wind and river discharge

(Estournel et al., 2001; Reffray et al., 2004) and more broadly the circulation over the entire Gulf of Lions (Estournel et al., 2003; Petrenko et al., 2008). The model configuration is based on a bipolar grid with a resolution of about 375 m on the Gulf of Lions shelf gradually increasing further offshore. This simulation is embedded in a configuration of the whole Mediterranean. The meteorological forcing is calculated with bulk formulas based on the hourly outputs of the ECMWF weather forecast model. The daily flow of the Rhône and the main rivers of the Gulf of Lions is extracted from the HYDRO database (www.hydro.eaufrance.fr). The simulation runs from July 2011 to April 2017. Daily averages of salinity and current are stored and presented here.

# RESULTS

## Reconnaissance of SGD Sites

There are several cooler surface water temperature plumes (~14°C) with respect to the warmer, ambient seawater (>15°C) identified by the previous airborne thermal infrared survey (Figure 2). Locations of cooler surface water temperatures may be influenced by groundwater inputs, as groundwater temperatures reflect the mean annual air temperature (Anderson, 2005). However, it cannot be excluded that other processes may also impact the temperature of surface waters in such coastal environments at this time of the day (e.g., cooling of shallow nearshore waters during the night; impact of waves, etc.). There are several small, low temperature plumes in bays in front of transects 1, 2, and 3 (Figures 2a,b,e). Finally, several plumes are visible at two locations where the coastline transitions to a beach near transect 4 (Sausset-les-Pins and Carry-le-Rouet, Figures 2c,d, respectively). For Sausset-les-Pins (Figure 2d), several plumes are visible between jetties, while no temperature differences are visible on other parts of the coastline. For Carry-le-Rouet (Figure 2c), the plumes are located throughout the coastline. These two locations were further investigated with salinity and radium isotope measurements (see below) to validate that the surface water temperature signatures were impacted by SGD. It is important to note that the airborne TIR flight (2012) was not conducted at the same time as the water sampling campaigns (2016–2017).

Two zones of surface water enrichments (>1 km of shoreline length) of $^{222}$Rn were observed (Figure 3). The local $^{222}$Rn enrichments reach activities up to two orders of magnitude higher (370–2,150 dpm 100 L⁻¹) than other in situ measurements taken along the coast (<370 dpm 100L⁻¹), indicating two prominent locations of SGD that include the Laurons bay west of Côte Bleue (T1) and Niolon bay (T7) east of Côte Bleue. These two areas correspond to sites where springs are known to discharge into the coastal seas and where large fractures exist, as noted by geological maps (Figure 1). A slight increase in $^{222}$Rn is also observed near transect 3, consistent with the previous airborne temperature signatures (Figure 2).

$^{223}$Ra and $^{224}$Raex activities were highest closest to the shoreline and decreased in activity with increasing distance offshore, with significant temporal variability (Figure 4; $^{223}$Ra not shown). Transect 4 was sampled during all campaigns; the sample collected closest to the shoreline shows large differences
<figure><img src="image_5.png" /></figure>

**FIGURE 3** | Radon activities (dpm 100 L⁻¹) acquired from a continuous shoreline survey on 27–28 October 2016. The dotted lines represent the location of the offshore surface water transects. Note that the color scale is not linear.

<figure><img src="image_6.png" /></figure>

**FIGURE 4** | <sup>224</sup>Ra<sub>ex</sub> activities (dpm 100 L<sup>-1</sup>) determined in surface waters along the transects off Côte Bleue in April 2016, October 2016, December 2016, and March 2017. The same color scale is used for the four panels.

in activity, with a maximum value in October (4.8 dpm 100 L⁻¹) and a minimum value in March (1.0 dpm 100 L⁻¹). During October 2016 and March 2017, the highest observed Ra activities were observed for the western most and the eastern most transects (Figure 4), which coincides with the shoreline areas where the <sup>222</sup>Rn signal was the highest (Figure 3). However, T1 was not sampled in April and December 2016 and T7 was not sampled in April 2016. The sampling stations farthest from the coast (~6–8 km) have <sup>224</sup>Ra<sub>ex</sub> activities ranging from below detection limit (i.e., no excess <sup>224</sup>Ra) to 3.6 dpm 100 L⁻¹, with

an average (±STD) value of 1.3 (±0.9) dpm 100 L⁻¹ (n = 16; Table 1). The <sup>228</sup>Ra activities in the coastal seas (<8 km) range between 1.7 and 24 dpm 100 L⁻¹, with a mean value of 4.2 ± 2.3 dpm 100 L⁻¹ (n = 98; Table 1), indicating that these coastal waters are enriched in <sup>228</sup>Ra with respect to open Mediterranean seawaters (typically 1–2.5 dpm 100 L⁻¹; Schmidt and Reyss, 1996; van Beek et al., 2009; Rodellas et al., 2015).

Surface water samples collected in several locations of cooler surface water temperatures within the first 5 m of the shoreline, as identified by the previous airborne TIR flight (Figure 2), display
<table><caption>**TABLE 1** | Average ± standard deviation of salinity and Ra isotopes in seawater (from 100 m to 8 km offshore), nearshore water (within 5 m offshore), karstic springs and minimum—maximum values of salinity, and Ra isotopes for shallow porewater (0.5 m depth) sampled from April 2016, October 2016, December 2016, and March 2017.</caption><thead><tr><th>Sample type</th><th>Sampling month</th><th>n</th><th>Salinity</th><th><sup>223</sup>Ra<br>dpm 100 L<sup>-1</sup></th><th><sup>224</sup>Ra<sub>ex</sub><br>dpm 100 L<sup>-1</sup></th><th><sup>228</sup>Ra<br>dpm 100 L<sup>-1</sup></th></tr></thead><tbody><tr><td rowspan="5">Seawater</td><td>Apr-16</td><td>21</td><td>38.0 ± 0.3</td><td>0.3 ± 0.1</td><td>2.3 ± 1.3</td><td>4.0 ± 4.0</td></tr><tr><td>Oct-16</td><td>36 (32)</td><td>38.2 ± 0.1</td><td>0.4 ± 0.2</td><td>3.4 ± 1.9</td><td>3.3 ± 1.6</td></tr><tr><td>Dec-16</td><td>14</td><td>–</td><td>0.2 ± 0.1</td><td>1.5 ± 0.7</td><td>3.6 ± 0.8</td></tr><tr><td>Mar-17</td><td>46 (31)</td><td>38.1 ± 0.4</td><td>0.3 ± 0.3</td><td>2.9 ± 3.5</td><td>4.0 ± 1.2</td></tr><tr><td>Average</td><td>117 (98)</td><td>38.1 ± 0.3</td><td>0.3 ± 0.2</td><td>2.9 ± 2.5</td><td>4.2 ± 2.3</td></tr><tr><td rowspan="5">Coastal water</td><td>Oct-16</td><td>9</td><td>37.1 ± 0.7</td><td>2.2 ± 1.7</td><td>31 ± 51</td><td>5.4 ± 2.6</td></tr><tr><td>Dec-16</td><td>14</td><td>21.8 ± 0.4</td><td>0.7 ± 0.5</td><td>8.1 ± 3.3</td><td>4.5 ± 0.8</td></tr><tr><td>Mar-17</td><td>18 (15)</td><td>37.7 ± 1.3</td><td>2.6 ± 3.4</td><td>45 ± 83</td><td>4.7 ± 0.9</td></tr><tr><td>Average</td><td>41 (38)</td><td>37.7 ± 0.9</td><td>1.1 ± 0.8</td><td>11 ± 5.0</td><td>4.8 ± 1.4</td></tr><tr><td>Oct-16</td><td>2</td><td>17.8–37.5</td><td>2.4–13.1</td><td>21–278</td><td>5.9–82</td></tr><tr><td rowspan="5">Porewater</td><td>Dec-16</td><td>5</td><td>18.2–38.1</td><td>6–39</td><td>28–212</td><td>2.1–16</td></tr><tr><td>Mar-17</td><td>5 (3)</td><td>13.7–38.1</td><td>3.6–18.4</td><td>39–324</td><td>2.4–18</td></tr><tr><td>Average</td><td>13 (11)</td><td>13.7–38.1</td><td>2–39</td><td>21–324</td><td>2.1–82</td></tr><tr><td>Average (May-17)</td><td>3 (3)</td><td>27.0 ± 0.9</td><td>37 ± 25</td><td>1,485 ± 658</td><td>826 ± 45</td></tr></tbody></table>

Shallow porewater and coastal water sample locations are presented in **Figure 2**; seawater sample locations are presented in **Figure 1**. Numbers within brackets are the number of samples analyzed for <sup>228</sup>Ra.

a salinity of 37.7 ± 0.9, which is not significantly different from offshore seawater (salinity = 38.1 ± 0.3; **Table 1**). The salinity measurements conducted in nearshore waters, therefore, do not confirm that the TIR plumes could be related to terrestrial groundwater inputs. These samples display relatively higher short-lived Ra activities (**Table 1**). <sup>228</sup>Ra activities in nearshore water samples (along the beach) and seawater samples (up to 8.5 km) were similar, with activities of 4.8 (±1.4) and 4.2 (±2.3) dpm 100 L<sup>-1</sup>, respectively. Shallow porewater samples (0.5 m) taken along the coastline were reduced in salinity (29.3 ± 9.7; minimum = 13.7) and higher in short-lived Ra activities, reflecting influence of a terrestrial groundwater endmember (**Figure 5**). <sup>223</sup>Ra and <sup>224</sup>Ra<sub>ex</sub> activities were between 2–39 and 21–324 dpm 100 L<sup>-1</sup> in porewaters while <sup>228</sup>Ra activities were 2.1–203 dpm 100 L<sup>-1</sup> (**Figure 5**). Finally, three karstic springs that discharge in Laurons bay at transect T1 showed average <sup>223</sup>Ra, <sup>224</sup>Ra<sub>ex</sub>, and <sup>228</sup>Ra activities of 37 (±25), 1,485 (±658), and 826 (±45) dpm 100 L<sup>-1</sup>, respectively (n = 3; **Table 1**).

Taken together, these lines of evidence suggest that the chemical enrichments observed along this 22 km coastline (<sup>223</sup>Ra, <sup>224</sup>Ra<sub>ex</sub>, <sup>228</sup>Ra, and <sup>222</sup>Rn; **Table 1**) are driven by SGD. This includes the discharge of terrestrial groundwater (e.g., via <sup>222</sup>Rn) and, although it is predominantly a rocky coast, a component of seawater circulation through permeable sediment, as evidenced by the porewaters collected in several beaches (**Figure 5**).

## Nutrient Concentrations

Average (±STD) nutrient concentrations for nearshore waters (<5 m offshore) and seawater samples are reported in **Table 2**, as well as minimum and maximum values of nutrient concentrations for shallow porewaters. During October 2016, nearshore waters (taken in areas of previously identified thermal infrared anomalies) had an average (±STD) concentration of

0.1 (±0.1) µM PO<sub>4</sub><sup>3-</sup>, 5.8 (±11.5) µM DSi, and 3.8 (±6.9) µM NO<sub>2</sub><sup>-</sup> + NO<sub>3</sub><sup>-</sup> (sum of NO<sub>2</sub><sup>-</sup> and NO<sub>3</sub><sup>-</sup>, hereafter NO<sub>3</sub><sup>-</sup> as this sum is >80% of NO<sub>3</sub><sup>-</sup> on average). Nearshore waters sampled during December 2016 had lower nutrient concentrations than the samples collected in October, with average values of 0.03 (±0.02), 1.1 (±0.4), and 0.8 (±0.3) µM for PO<sub>4</sub><sup>3-</sup>, DSi and NO<sub>3</sub><sup>-</sup>, respectively. In general, offshore seawater samples collected were greatest for PO<sub>4</sub><sup>3-</sup>, DSi, and NO<sub>3</sub><sup>-</sup> during April and March, while average concentrations were lower in December and October. In general, the concentrations of DSi and NO<sub>3</sub><sup>-</sup> decreased with increasing salinity, and thus distance offshore, similar to <sup>223</sup>Ra and <sup>224</sup>Ra<sub>ex</sub> (**Figures 4–6**). Most PO<sub>4</sub><sup>3-</sup> concentrations in seawater were below the method detection limit (0.01 µM).

Porewater concentrations were approximately one order of magnitude higher than coastal surface waters, with values between 0.01 and 4.0 µM for PO<sub>4</sub><sup>3-</sup>, 1.9 and 133 µM for DSi, and 22 and 194 µM for NO<sub>3</sub><sup>-</sup> over the sampled salinity gradient (**Figure 6**). The variable average nutrient concentrations between sampling months (**Table 2**) may reflect changes in porewater salinity, rather than a seasonal endmember (**Figure 6**). The average NO<sub>3</sub><sup>-</sup> and DSi concentrations of the springs sampled in Laurons bay in front of transect T1 are 2.3 (±3.3) 113 (±40) µM, respectively (n = 3; mean salinity of 27). PO<sub>4</sub><sup>3-</sup> was not analyzed for these springs.

## DISCUSSION

### Estimate of the SGD Fluxes

#### Determination of the Horizontal Eddy Diffusivity Coefficients K<sub>h</sub>

All of the surface water transects were used to estimate horizontal eddy diffusivity coefficients, K<sub>h</sub> (**Figure 7**). We do not observe any significant difference between the <sup>223</sup>Ra and <sup>224</sup>Ra<sub>ex</sub> derived
<figure><img src="image_7.png" /></figure>

**FIGURE 5** | (A) $^{224}$Ra$_{ex}$ (dpm 100 L$^{-1}$), (B) $^{223}$Ra (dpm 100 L$^{-1}$), and (C) $^{228}$Ra (dpm 100 L$^{-1}$) along Côte Bleue as a function of salinity, for all sampling dates. Samples are categorized by water types (karstic spring, porewater, nearshore water within 5 m from the shoreline and seawater from 100 m to 8 km offshore). Right-hand side panels are zooms between salinity 35 and 40.

$K_h$. Here, we choose to report the $K_h$ values derived from the $^{224}$Ra$_{ex}$ activities that display lower uncertainties (Figure 7). The slopes (together with their associated uncertainty) of the linear relationships between ln($^{224}$Ra) and offshore distance are reported for each campaign on Figure 7, together with the correlation coefficients *r* and *p*-values. For a given season, the $^{224}$Ra$_{ex}$ activities reported for the different transects are similar and decrease with increasing offshore distance. Therefore, we report a single $K_h$ value for each season that is deduced from the slope of the linear relationship (following Equation 2). In doing so, the number of data points is increased and the significance of the correlation coefficient is thus improved. We calculate the error of the slope, which we use to determine the uncertainty on the $K_h$ estimate, unlike many studies that report $K_h$ values without considering any associated uncertainty (Moore, 2000a; Tamborski et al., 2018). Note that in March 2017, two different trends were observed (transects T1, T2, T7 showing higher Ra activities than transects T4, T5, T6). In this latter case, two $K_h$ estimates were determined from the two trends (Figure 7). However, the slopes obtained are the same (within error bars)

and therefore the $K_h$ estimates are not significantly different. This suggests that the two trends observed in March 2017 were not the result of different offshore mixing characteristics but are rather explained by differences in the absolute Ra activities at the coast (i.e., higher $^{224}$Ra activities in March 2017 for T1 and T7 that are located offshore of springs, but also for T2). The significance of the correlation coefficients suggests that the offshore dispersion of the Ra activities may indeed be approximated by a 1D diffusive mixing model. However, we have no *in situ* information that would support the assumption of negligible advection.

The mean $K_h$ for the two transects sampled in April 2016 is 39 (±22) m² s⁻¹ (Table 3). In October 2016, 5 transects were investigated and the mean $K_h$ is 96 (±44) m² s⁻¹. The highest $K_h$ was estimated for December 2016, with a $K_h$ of 184 (±112) m² s⁻¹, which could be related to the winter conditions (increased mixing). In March 2017, the $K_h$ values associated with the two groups of transects were 56 (±42) and 52 (±35) m² s⁻¹, for transects T1, T2, T7, and transects T4, T5, T6, respectively. The $K_h$ values do not exhibit a significant temporal variability when considering their associated uncertainties.
<table><caption>**TABLE 2** | Average ± standard deviation of nutrient concentrations in seawater (from 100 m to 8 km offshore), nearshore water (within 5 m offshore), karstic springs and minimum—maximum values for shallow porewater (0.5 m depth) sampled from April 2016, October 2016, December 2016, and March 2017, and karstic springs (Laurons bay).</caption><thead><tr><th>Sample type</th><th>Sampling month</th><th>n</th><th>PO<sub>4</sub><sup>3-</sup> µM</th><th>Si(OH)<sub>4</sub> µM</th><th>NO<sub>2</sub><sup>-</sup> + NO<sub>3</sub><sup>-</sup> µM</th></tr></thead><tbody><tr><td rowspan="5">Seawater</td><td>Apr-16</td><td>21</td><td>0.04 ± 0.02</td><td>2.1 ± 0.5</td><td>1.3 ± 0.5</td></tr><tr><td>Oct-16</td><td>36</td><td>0.02 ± 0.02</td><td>1.4 ± 0.3</td><td>0.3 ± 0.2</td></tr><tr><td>Dec-16</td><td>14</td><td>0.02 ± 0.02</td><td>1.5 ± 0.4</td><td>0.7 ± 0.3</td></tr><tr><td>Mar-17</td><td>46</td><td>0.03 ± 0.02</td><td>1.6 ± 0.5</td><td>1.5 ± 0.5</td></tr><tr><td>Average</td><td>117</td><td>0.02 ± 0.02</td><td>1.6 ± 0.5</td><td>1.0 ± 0.7</td></tr><tr><td rowspan="4">Coastal water</td><td>Oct-16</td><td>9</td><td>0.09 ± 0.10</td><td>5.8 ± 11.5</td><td>3.8 ± 6.9</td></tr><tr><td>Dec-16</td><td>14</td><td>0.03 ± 0.02</td><td>1.1 ± 0.4</td><td>0.8 ± 0.3</td></tr><tr><td>Mar-17</td><td>18</td><td>0.08 ± 0.13</td><td>7.0 ± 15.4</td><td>1.7 ± 2.5</td></tr><tr><td>Average</td><td>41</td><td>0.06 ± 0.09</td><td>4.2 ± 10.8</td><td>1.8 ± 3.8</td></tr><tr><td rowspan="4">Porewater</td><td>Oct-16</td><td>2</td><td>1.1–2.0</td><td>87–109</td><td>44–194</td></tr><tr><td>Dec-16</td><td>5</td><td>0.01–4.0</td><td>2–58</td><td>2–56</td></tr><tr><td>Mar-17</td><td>5</td><td>0.3–1.2</td><td>6–133</td><td>22–41</td></tr><tr><td>Average</td><td>12</td><td>0.01–4.0</td><td>2–133</td><td>22–194</td></tr><tr><td rowspan="2">Karstic spring</td><td>May-17</td><td>3</td><td>–</td><td>82–157</td><td>0.02–4.26</td></tr><tr><td>Average</td><td>3</td><td>–</td><td>113 ± 40</td><td></td></tr></tbody></table>

## Mixing Model: K<sub>n</sub> * Nutrient Offshore Gradient

For numerous transects, there is no significant nutrient gradient (for example, no DSi gradient for T5 in April 2016; **Figure 8A**); in this case, a slope is not reported and the nutrient flux is not estimated with this method. Statistically significant surface water gradients (and gradient uncertainties) are reported in **Figure 8** for DSi, and in **Figure 9** for NO₃⁻. Surface water DSi and NO₃⁻ gradients, with estimated uncertainties, are summarized by sampling season in **Table 3**. It bears mention that NO₃⁻ gradients may not capture the complete dissolved inorganic nitrogen gradient, which may comprise a significant pool of NH₄⁺. We do not report PO₃⁻ gradients, as the concentration of PO₃⁻ is very low and uniform in the Mediterranean seawater samples (**Table 2**).

²²⁴Ra<sub>ex</sub>-derived horizontal eddy diffusivity coefficients (K<sub>h</sub>; **Figure 7**) were thus multiplied by (statistically significant) surface water nutrient gradients (**Figures 8, 9; Table 3**) and the water depth impacted by SGD (5–10 m). Corresponding DSi fluxes varied between 0.9 (±0.9)–1.7 (±1.8) *10³ mol d⁻¹ km⁻¹ of shoreline and 6.3 (±4.5)–13 (±9.4) *10³ mol d⁻¹ km⁻¹ of shoreline, estimated in April 2016 and December 2016, respectively (**Figure 10**). NO₃⁻ fluxes were between 0.5 (±0.4)–1.1 (±0.7) *10³ and 3.4 (±3.5)–6.9 (±6.9) *10³ mol d⁻¹ km⁻¹ of shoreline for the same seasons (**Figure 10**).

April 2016 exhibited the lowest mean DSi flux, with only one statistically significant seawater DSi transect (T4; **Figure 10**). In comparison, the monthly average DSi fluxes were maximum in October and December 2016. We do observe significant temporal variability among individual transects that were repeatedly sampled (transects 4 and 7). DSi fluxes were four times higher in October than in April for transect 4, and up to 8 times higher in

December for that same transect. Temporal variability in nutrient flux for these repeated transects may be driven by temporally variable precipitation. The total precipitation registered (during 4 days) 12 days before the December sampling equaled 88 mm while 33 mm of rainfall was registered (during 4 days) 11 days prior to the October sampling. Monthly averaged NO₃⁻ fluxes are similar between the different sampling periods (**Figure 10; Table 3**). NO₃⁻ fluxes estimated for transect 4 were similar in April and October. The average (±STD) DSi and NO₃⁻ fluxes, for all the transects in which a nutrient gradient was estimated, equal to 6.2 (±5.0) *10³ and 4.0 (±2.0) *10³ mol d⁻¹ km⁻¹, respectively (**Figure 10; Table 3**).

It remains to be seen if SGD is the sole source of the observed nutrient gradients. Nutrient inputs from the Rhône River, the largest river that enters the Mediterranean Sea and whose river mouth is 15 km from the western shoreline of Côte Bleue, may be deflected eastward toward Côte Bleue on rare occasions (Gatti et al., 2006; Fraysse et al., 2014). Surface salinity maps have been generated using the SYMPHONIE model to study the fate of the Rhône river plume into the Gulf of Lions at the dates of the sampling campaigns. The plume was usually deflected toward the West (April 2016, October 2016, December 2016). In contrast, intrusion of Rhône river water (eastward transport) was observed on 15–16 March 2017 (that is, 1 week before the March 2017 sampling campaign), with diluted river waters reaching the entire Côte Bleue coastline (**Figure 11**). On 21 and 22 March, the impact of river waters was less important, with the exception of the western Côte Bleue that was still potentially impacted by these waters (**Figure 11**). The intrusion of Rhône river waters may thus have impacted the nutrient and Ra distributions offshore Côte Bleue during the March 2017 campaign. There are no other known nutrient sources along Côte Bleue. On the other hand, interpretation of offshore nutrient gradients (or lack thereof) is subject to a sound understanding of the uptake rate of inorganic nutrients in the water column, in addition to other known nutrient sources and sinks. As shown in Monterey Bay, nutrient loads driven by SGD can be quickly utilized by phytoplankton (Lecher et al., 2015). Along Côte Bleue, the apparent surface water ages derived from ²²⁴Ra/²²⁸Ra activity ratios (Moore, 2000b) range from several days (open coastline) to ca. ten days in semi-enclosed bays like Laurons bay located in front of T1. Nutrient uptake cannot be excluded over this time scale and could impact the offshore nutrient gradient. Thus, we rely on several different methods (e.g., Radon survey; ²²⁸Ra offshore gradients), in addition to the nutrient gradient method, to constrain SGD along Côte Bleue.

## ²²⁸Ra Gradient Method

The regional view of the entire dataset is that the ²²⁸Ra activities do not decrease with increasing distance offshore (**Figure 12**). Only two transects among the 16 transects investigated here display statistically significant offshore ²²⁸Ra gradients (*p* < 0.01) and we thus derived SGD fluxes from these data only. The ²²⁸Ra gradient for T7 in December 2016 is −0.164 (±0.065) dpm 100 L⁻¹ km⁻¹ while the ²²⁸Ra gradient for T5 in March 2017 is −0.225 (±0.058) dpm 100 L⁻¹ km⁻¹. The ²²⁸Ra fluxes are estimated by combining these gradients with the ²²⁴Ra-derived
<figure><img src="image_8.png" /><figcaption>**FIGURE 6** | (A) DSi (µM), (B) NO₂⁻ + NO₃⁻ (µM), and (C) PO₄³⁻ (µM) along Côte Bleue as a function of salinity, for all sampling dates. Samples are categorized by water types (karstic spring, porewater, nearshore water within 5 m from the shoreline, and seawater from 100 m to 8 km offshore). Right-hand side panels are zooms between salinity 35 and 40.</figcaption></figure>

horizontal eddy diffusivity coefficient corresponding to the appropriate month (section Determination of the Horizontal Eddy Diffusivity Coefficients K<sub>h</sub>). The ²²⁸Ra fluxes are thus 2.6 (±1.9) *10⁸ and 1.8 (±0.8) *10⁸ dpm d⁻¹ km⁻¹ for December (T7) and March (T5), respectively.

This method requires the use of the SGD endmember to determine the volumetric SGD flux (obtained by dividing the ²²⁸Ra fluxes by the ²²⁸Ra activity of the endmember). The nutrient flux can finally be estimated by multiplying this water flux with the nutrient concentration of the endmember. This makes this method particularly sensitive to the endmember, which is not easy to define in this region (spring at its outlet vs. porewater; Figures 5, 6). Using the karstic springs as an endmember (²²⁸Ra of 826 ± 45 dpm 100 L⁻¹; n = 3), this yields a volumetric SGD flux for T7 in December 2016 of 1.6 (±1.2)–3.2 (±2.3) *10⁴ m³ d⁻¹ km⁻¹ using an impacted layer of 5 and 10 m, respectively. The volumetric SGD fluxes for T5 in March 2017 are 0.7 (±0.5)–1.4 (±1.0) *10⁴ m³ d⁻¹ km⁻¹, respectively.

Using the mean DSi concentration in the spring (113 µM), the DSi fluxes thus estimated are 1.8 (±1.4)–3.6 (±2.9) *10³ mol d⁻¹ km⁻¹ for T7 in December and 0.8 (±0.6)–1.6 (±1.2) *10³ mol d⁻¹ km⁻¹ for T5 in March (for 5 and 10 m impacted depths, respectively). These estimations are on the same order of magnitude as the estimations reported in section Mixing Model: K<sub>h</sub> * Nutrient Offshore Gradient (Figure 10). Using the NO₃⁻ concentrations of the brackish springs (2.4 ± 2.6 µM; two springs considered out of the three), the NO₃⁻ fluxes estimated from the ²²⁸Ra gradient are 0.04 (±0.06)–0.08 (±0.10) *10³ and 0.02 (±0.03)–0.03 (±0.04) *10³ mol d⁻¹ km⁻¹ for T7 and T5, respectively. The estimations of NO₃⁻ fluxes are lower than the estimations made in section Mixing Model: K<sub>h</sub> * Nutrient Offshore Gradient. Note that the NO₃⁻ concentrations of the brackish springs are surprisingly low (2.4 ± 2.6 µM) compared to other Mediterranean karstic springs (e.g., Rodellas et al., 2015). This pattern could explain the NO₃⁻ fluxes being lower than the estimations made in section Mixing Model: K<sub>h</sub> * Nutrient
<figure><img src="image_9.png" /><figcaption>**FIGURE 7** | Plot of ln $^{224}$Ra$_{ex}$ (dpm 100 L$^{-1}$) as a function of distance offshore. The transects are arranged by sampling date (April 2016, October 2016, December 2016, and March 2017). Each transect is presented by a unique color. For each season, we report the correlation coefficient *r*, slope *m* (and associated uncertainty), *p*-value, and the number of values (*n*) used to build the linear regression.</figcaption></figure>

<table><caption>**TABLE 3** | Monthly averaged horizontal eddy diffusivity coefficients ($K_h$; m$^2$ s$^{-1}$), nutrient gradients (µmol L$^{-1}$ km$^{-1}$) and nutrient fluxes (mol d$^{-1}$ km$^{-1}$) along Côte Bleue.</caption><thead><tr><th rowspan="2">Sampling time</th><th>K<sub>h</sub> (<sup>224</sup>Ra)</th><th rowspan="2">DSi gradient</th><th>NO<sub>2</sub><sup>-</sup> + NO<sub>3</sub><sup>-</sup> gradient</th><th>DSi flux</th><th>(NO<sub>2</sub><sup>-</sup> + NO<sub>3</sub><sup>-</sup>) flux</th></tr><tr><th>m<sup>2</sup> s<sup>-1</sup></th><th>µmol L<sup>-1</sup> km<sup>-1</sup></th><th>mol d<sup>-1</sup> km<sup>-1</sup></th><th>mol d<sup>-1</sup> km<sup>-1</sup></th></tr></thead><tbody><tr><td>Apr-16</td><td>39 ± 22</td><td>-0.05 ± 0.04</td><td>-0.09 ± 0.04</td><td>1.7E+03 ± 1.8E+03</td><td>3.2E+03 ± 2.3E+03</td></tr><tr><td>Oct-16</td><td>96 ± 44</td><td>-0.08 ± 0.03</td><td>-0.04 ± 0.03</td><td>6.4E+03 ± 3.8E+03</td><td>3.7E+03 ± 2.7E+03</td></tr><tr><td>Dec-16</td><td>184 ± 112</td><td>-0.08 ± 0.03</td><td>-0.04 ± 0.03</td><td>1.3E+04 ± 9.0E+03</td><td>6.9E+03 ± 6.9E+03</td></tr><tr><td>Mar-17</td><td>54 ± 38</td><td>-0.07 ± 0.01</td><td>-0.05 ± 0.02</td><td>3.5E+03 ± 2.5E+03</td><td>2.2E+03 ± 1.8E+03</td></tr><tr><td>Average</td><td></td><td></td><td></td><td>6.2E+03 ± 5.0E+03</td><td>4.0E+03 ± 2.0E+03</td></tr></tbody></table>

The flux estimations reported here were made using 10 m as the water depth impacted by SGD.

Offshore Gradient and questions the validity of the NO$_3^-$ endmember used for this method. The spring endmember may be dominated by NH$_4^+$, which we did not measure in this study.

Alternatively, considering the maximum NO$_3^-$ concentrations of the porewater samples (83 µM corresponding to S = 21.1, with a $^{228}$Ra activity of 125 dpm 100 L$^{-1}$; Table S1) as the endmember leads to NO$_3^-$ fluxes of 8.7 (±6.4)–17 (±13) *10$^3$ and 3.8 (±2.8)–7.7 (±5.5) *10$^3$ mol d$^{-1}$ km$^{-1}$ for T7 and T5, respectively. These estimations are in better agreement, but greater than the estimations of section Mixing

Model: K$_h$ * Nutrient Offshore Gradient (Figure 10). Using the maximum DSi concentration of the porewater samples (133 µM corresponding to a salinity of 13.7; with a $^{228}$Ra activity of 18.2 dpm 100 L$^{-1}$) as the endmember leads to DSi fluxes of 95 (±69)–190 (±150) *10$^3$ mol d$^{-1}$ km$^{-1}$ of shoreline and 42 (±30)–84 (±65) *10$^3$ mol d$^{-1}$ km$^{-1}$ for T7 and T5, respectively (calculated for 5 and 10 m layers). Considering this porewater DSi endmember leads to estimates that are one to two orders of magnitude higher than the estimate reported in section Mixing Model: K$_h$ * Nutrient Offshore Gradient. This highlights
<figure><img src="image_10.png" /><figcaption>**FIGURE 8** | Plot of DSi (µM) as a function of distance offshore. The transects are arranged by sampling date ([A] April 2016, [B] October 2016, [C] December 2016, and [D] March 2017). Each transect is presented as a unique color. The dotted lines represent the linear regressions with significant *p*-values. The correlation coefficient *r*, slope *m* (and associated uncertainty), and *p*-values are listed for each transect (same color code as the transect symbols).</figcaption></figure>

the importance of properly selecting the SGD endmember in
determining nutrient fluxes from offshore ²²⁸Ra gradients.

We should be able to estimate an SGD-driven PO₄³⁻ flux using this method, whereas we could not from the method used in section Mixing Model: K<sub>h</sub> * Nutrient Offshore Gradient because there was no observable surface water PO₄³⁻ gradient. Considering the maximum PO₄³⁻ concentrations of the porewater samples (4.04 µM; **Table 2**), PO₄³⁻ fluxes are 0.8 (±0.7)–1.6 (±1.4) *10²* and 0.4 (±0.3)–0.7 (±0.6) *10²* mol d⁻¹ km⁻¹ for T7 and T5, respectively. PO₄³⁻ was not measured in the karstic springs sampled in the bay at the head of transect 1. In comparison, Tamborski et al. (2018) estimated an SGD-driven DIP flux along the sandy, alluvial shoreline of La Palme lagoon of 0.2 (±0.1) *10²* mol d⁻¹ km⁻¹ of shoreline.

In the following, we prefer to rely on the flux estimates
determined by combining horizontal eddy diffusivity coefficients
with surface water nutrient gradients (section Mixing Model: K<sub>h</sub> * Nutrient Offshore Gradient), particularly because knowledge of
the SGD endmember is not requisite for this method. In addition,
the ²²⁸Ra distribution during the different campaigns did not
show a decreasing trend with increasing distance offshore, which
prevented us from using the ²²⁸Ra mixing model for the majority
of the sampling campaigns. The Rhône River mouth is located

15 km west of Côte Bleue and represents the largest freshwater input into the Mediterranean Sea (mean flow of 1,700 m³ s⁻¹; Pont et al., 2002). The river plume is usually deflected westward (Estournel et al., 2001). However, Gatti et al. (2006) and Fraysse et al. (2014) showed that intrusion of Rhône river diluted waters could in some cases reach the Bay of Marseille, as a consequence of different physical processes (e.g., wind stress; presence of mesoscale eddy). In particular, the combination of an intrusion event with southeastern winds is expected to impact the Côte Bleue coastline (Fraysse et al., 2014). Such intrusion was shown to impact the biogeochemical functioning of the coastal waters by transporting large quantities of nutrients into the bay. Because ²²⁸Ra displays a relatively long half-life (5.75 y), its distribution offshore of Côte Bleue may thus be impacted by sporadic intrusions of Rhône River waters. In March 2017, the ²²⁸Ra activities are especially high at T1 near the coastline (~25 dpm 100 L⁻¹; **Figure 12**). During that same campaign, the nutrient concentrations were also high in the coastal waters. The transects located along the western half of the Côte Bleue coastline showed elevated NO₃⁻ concentrations during this period (T1–T4; **Figure 9**). As discussed above (section Mixing Model: K<sub>h</sub> * Nutrient Offshore Gradient), the surface salinity maps generated using the SYMPHONIE model indicate intrusion
<figure><img src="image_11.png" /></figure>

**FIGURE 9** | Plot of NO₂⁻ + NO₃⁻ (µM) as a function of distance offshore. The transects are arranged by sampling date [(A) April 2016, (B) October 2016, (C) December 2016, and (D) March 2017]. Each transect is presented as a unique color. The dotted lines represent the linear regressions with significant *p*-values. The regression coefficient *r*, slope *m* (and associated uncertainty), and *p*-values are listed for each transect (same color code as the transect symbols). Note that the plots have different y-axes scale.

of Rhône river waters that reached the entire Côte Bleue coastline on 15–16 March 2017 (Figure 11). One week later, on 21 and 22 March (March campaign), the western part of Côte Bleue is still impacted by these waters (Figure 11), which likely explain the observed nutrient and ²²⁸Ra patterns. The vertical profiles of salinity (Figure S2) indicate waters with lower salinity (<38) at T1, T2, and T4 during the March 2017 campaign that support this hypothesis. ²²⁴Ra, with its shorter half-life (3.66 d), may be less impacted by such intrusions that occur on average 7.6 times per year (Fraysse et al., 2014), unless sampling is performed just after an intrusion event. This is an additional argument to determine SGD from the method that combines K_h with the gradient of the chemical element of interest, rather than the method using the ²²⁸Ra gradient.

## Significance of the SGD Fluxes

Tamborski et al. (2018) estimated DSi and DIN fluxes along La Palme lagoon, a sandy alluvial stretch of the French Mediterranean coastline, of (2.4 ± 1.4) *10³ mol Si d⁻¹ km⁻¹ and (5.7 ± 3.2) *10³ mol N d⁻¹ km⁻¹ of shoreline into the Mediterranean Sea (Figure 13). These fluxes are similar in magnitude to the DSi and DIN fluxes reported here (Figure 10)

despite the different lithology of the two areas (unconsolidated sediment vs. karst systems). For the entire Mediterranean basin, Rodellas et al. (2015) reported DSi and DIN fluxes of (0.1–27) *10³ mol d⁻¹ km⁻¹ of shoreline and (0.1–64) *10³ mol d⁻¹ km⁻¹ of shoreline, respectively, in relative agreement with the fluxes estimated here (Figure 13). Our estimations are on the same order of magnitude as other karst aquifer SGD studies in the Mediterranean Sea. Weinstein et al. (2011) estimate a terrestrial SGD flux discharging in Dor Bay, southeastern Mediterranean Sea, of 1.5 *10³ and 1.4 *10³ mol d⁻¹ km⁻¹ of DSi and DIN, respectively. In different coves along Majorca Island, SGD supplies a DSi flux ranging from 22 to 13,000 mol d⁻¹ km⁻¹ and DIN fluxes of 21 to 6,500 mol d⁻¹ km⁻¹ (Tovar-Sánchez et al., 2014). These examples suggest that Côte Bleue is similar to other Mediterranean karstic areas.

We estimated nutrient fluxes supplied by local rivers using river gauging stations (hydro.eaufrance.fr) and monthly chemical element analyses (sierm.eaurmc.fr), where the nutrient flux is equivalent to the river discharge multiplied by the nutrient concentration. Here, we first compared the nutrient flux driven by SGD to the Huveaune River (flowing in Marseille; Figure 1) using the gauging station (Y4424040) in Aubagne and the
<figure><img src="image_12.png" /></figure>

**FIGURE 10** | Fluxes of DSi and NO₃⁻ (mol d⁻¹ km⁻¹) estimated along Côte Bleue for each transect using Ra isotopes (K_h * nutrient gradient). The mean fluxes (±standard deviation) are also reported for each season (April 2016, October 2016, December 2016, and March 2017) and are denoted AVG.

<figure><img src="image_13.png" /></figure>

**FIGURE 11** | Daily-averaged surface salinity and current simulated with the SYMPHONIE model for: (A) 27 October 2016, (B) 16 March 2017, (C) 22 March 2017. Note the change in the isocontours spacing at 37. These figures illustrate two different patterns for the Rhône river plume that enters into the Gulf of Lions: in most cases, the Rhône river waters are transported westward (October 2016 situation as an example) and in some rare cases, they are deflected eastward toward the Côte Bleue coastline (March 2017 situation). The arrow locates the Rhône river and the open circle indicates the Côte Bleue coastline.

chemical analyses station (06198100) in Marseille. It is important to note that the Huveaune River is over 15 km away from the eastern section of our study site; these waters do not impact our investigated shorelines. We then compared SGD fluxes in Côte Bleue to the largest river discharging in the French Mediterranean Sea: The Rhône River. We use the gauging station in Tarascon (Station: V7200015) and the chemical analyses made in Roquemaure, the most downstream sampling station (Station: 06121500). We compared these river fluxes to the SGD-driven NO₃⁻ and DSi fluxes estimated from the method described in section Mixing Model: K_h * Nutrient Offshore Gradient. The SGD-driven nutrient fluxes are multiplied by the coastal length of Côte Bleue impacted by terrestrial SGD; here we take a shoreline

length of 5 km, where the coastal ²²²Rn activities were the highest (Figure 3). Note that these fluxes will represent upper limits, as the mean SGD-driven nutrient fluxes are estimated only for transects in which we observed significant offshore gradients.

SGD represents between 3 and 22 times the DSi and NO₃⁻ fluxes driven by the Huveaune River, with significant temporal variability (Table 4). While the DSi flux driven by SGD is 4 times greater than the DSi flux driven by the Huveaune River in March 2017, SGD is 22 times greater than the river fluxes in December 2016. NO₃⁻ fluxes follow the same pattern, with SGD-driven fluxes slightly higher than the Huveaune river flux in March 2017 and reaching 5 times the river flux in October 2016. In comparison, the DSi fluxes driven by SGD is only between 0.1
<figure><img src="image_14.png" /><figcaption>**FIGURE 12** | Plot of $^{228}$Ra (dpm 100 L$^{-1}$) as a function of distance offshore. The transects are arranged by sampling date ([A] April 2016, [B] October 2016, [C] December 2016, and [D] March 2017). Each transect is presented as a unique color. Note that the plots have different y-axes scale. The dotted lines represent the linear regressions with significant *p*-values.</figcaption></figure>

and 1.4% of the DSi fluxes driven by the Rhône River (Table 4). The $NO_3^-$ flux driven by SGD is only between 0.1 and 0.3% of the Rhône River $NO_3^-$ inputs (Table 4). Nutrient inputs of SGD may still be significant compared to the Rhône River, as SGD is the only terrestrial input of nutrients to Côte Bleue when the Rhône River plume does not enter the Bay of Marseille and Côte Bleue.

CONCLUSION

Characterization of water and nutrient fluxes from coastal karst aquifers remains a significant challenge (Montiel et al., 2018). We used two different geochemical tracer methods to estimate SGD fluxes along a karst coastline (Côte Bleue, French Mediterranean Sea) where the presence of terrestrial groundwater discharge was suspected from the thermal infrared signature of nearshore waters and from a radon survey that was conducted along the coastline. Horizontal eddy diffusivity coefficients ($K_h$), derived from offshore $^{224}$Ra$_{ex}$ transects, combined with offshore nutrient gradients provided the most consistent flux results, as we could derive SGD nutrient fluxes for almost all the investigated transects during all sampling periods. Knowledge of the SGD endmember is not required for this method, which constitutes

a strong advantage in such an area where the endmember is difficult to constrain. This method is only applicable when the chemical element of interest behaves conservatively within the timeframe of coastal mixing processes, and when there is no significant external nutrient (or Ra) sources. In contrast, offshore $^{228}$Ra gradients combined with $K_h$ yielded significant flux estimations for only 2 of the investigated transects (out of 16, over 4 seasons). This approach is also disadvantageous for Côte Bleue in that it requires a well-constrained SGD endmember, which is difficult to determine in this specific region (karstic spring vs. porewaters).

The $^{228}$Ra and nutrient distributions observed during March 2017 suggest that the coastal waters may have been impacted by the intrusion of Rhône River diluted waters that were deflected eastward, a phenomenon that only occurs ~8 times a year (Fraysse et al., 2014), further complicating the use of long-lived $^{228}$Ra as a tracer for SGD along Côte Bleue. This is in contrast with the dominant situation when the plume is deflected westward and therefore does not impact the investigated region. The Rhône River constitutes the largest nutrient source to the Gulf of Lions and can support 23–69% of the primary production (Ludwig et al., 2009). We estimate that SGD along Côte Bleue
<figure><img src="image_15.png" /><figcaption>**FIGURE 13** | DSi and NO₃⁻ fluxes determined in this study along Côte Bleue. As a comparison, we report the (i) DSi fluxes determined by Tamborski et al. (2018) offshore La Franqui coastline (western Gulf of Lions), (ii) the DSi and N fluxes determined by Rodellas et al. (2015) using a Ra mass balance applied to the entire Mediterranean Sea (the fluxes reported on the figure have been normalized to the Mediterranean shore length -64,000 km—according to Rodellas et al., 2015), and (iii) the DSi and N fluxes associated with the Rhône river.</figcaption></figure>

<table><caption>**TABLE 4** | Estimation of nutrient fluxes (mol d⁻¹) driven by SGD, the Huveaune River and the Rhône River in April, October, and December 2016, and March 2017.</caption><thead><tr><th></th><th></th><th>SGD<br>(*10³ mol d⁻¹)</th><th>Huveaune river<br>(*10³ mol d⁻¹)</th><th>Rhône river<br>(*10⁶ mol d⁻¹)</th></tr></thead><tbody><tr><td rowspan="4">DSi</td><td>Apr-16</td><td>9 ± 9</td><td>1.7</td><td>13</td></tr><tr><td>Oct-16</td><td>32 ± 12</td><td>1.4</td><td>5.3</td></tr><tr><td>Dec-16</td><td>63 ± 45</td><td>2.9</td><td>4.4</td></tr><tr><td>Mar-17</td><td>17 ± 2</td><td>4.4</td><td>6.5</td></tr><tr><td rowspan="4">NO₃⁻ +NO₂⁻</td><td>Apr-16</td><td>16 ± 7</td><td>5.4</td><td>14</td></tr><tr><td>Oct-16</td><td>19 ± 11</td><td>3.9</td><td>9.4</td></tr><tr><td>Dec-16</td><td>34 ± 35</td><td>10.1</td><td>10</td></tr><tr><td>Mar-17</td><td>11 ± 9</td><td>13.3</td><td>14</td></tr></tbody></table>

DSi and NO₃⁻ fluxes reported here were estimated using Ra isotopes (Kₐ * nutrient gradients) considering 10 m as the water depth impacted by SGD. Note that the Rhône River fluxes are reported in 10⁶ mol d⁻¹.

supplies ≤1% of the DSi and NO₃⁻ inputs of the Rhône River. However, in the absence of eastward deflected Rhône River diluted waters, SGD is the only known nutrient source to the nearshore area of Côte Bleue. It is interesting to note that the nutrient fluxes reported here are similar in magnitude compared with the fluxes quantified along the sandy beach of La Franqui, in the western Gulf of Lions (Tamborski et al., 2018), despite the different lithology of the two areas (karst systems vs. unconsolidated sediment). We recommend that multiple

methods be used in future studies investigating SGD and nutrient cycling from coastal karst aquifers.

## DATA AVAILABILITY STATEMENT

All datasets generated for this study are included in the article/Supplementary Material.

## AUTHOR CONTRIBUTIONS

SB, JT, MS, and PB analyzed the samples for radium isotopes. OR, CC, and TS conducted the radon surveys (in situ analysis). SB, JT, PB, MS, TS, OR, CC, and EL contributed to the collection of samples at sea or along the coastline. OC, PC, and MP-P analyzed the nutrient concentrations. SB, JT, TS, and PB computed the TIR images. SB, JT, and PB calculated the SGD fluxes based on the radium isotopes. SB, JT, PB, TS, OR, CC, PC, and MP-P contributed to the interpretation of the data. SB, JT, PB, TS, OR, and PC contributed to the writing of the manuscript. CE generated the salinity and circulation maps using the SYMPHONIE model and contributed to the sections related to the interpretation of these maps.

## FUNDING

The Ph.D. thesis fellowship of SB and the postdoctoral fellowship of JT were supported by FEDER funded by Europe and Région
Occitanie Pyrénées-Méditerranée (SELECT project). This project was funded by (i) ANR-MED-SGD (ANR-15-CE01-0004; PB) and (ii) CNES for funding the airborne TIR images acquired in 2012 as part of the Geomether project (PI: Pascal Allemand, PB being responsible for the acquisition of TIR images in that project).

## ACKNOWLEDGMENTS

We thank Virginie Sanial for her participation to the survey that allowed us to acquire airborne TIR images. We thank captain and crew of RV Antédon II for help during sampling at sea. We thank Dorian Guillemain, Nagib Bhairy, Deny Malengros, and Christian Grenz at MIO for providing the CTD data. We thank M. Plantevin at Saumaty harbor. We thank European Union and Région Occitanie Pyrénées-Méditerranée for supporting the LAFARA underground laboratory through a FEDER funding (SELECT project). We are grateful to EDF (Electricité De France) for allowing us to run our germanium detectors in the tunnel of Ferrières, Ariège. We thank Caroline Ulses and Ivane Pairaud

for constructive discussions. We thank the two reviewers and the associate editor Henrietta Dulai for their constructive comments that allowed us to improve the quality of the paper.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2019.00205/full#supplementary-material

**Figure S1** | Daily precipitation recorded at Marseille-Marignane airport (data from meteociel.fr). Periods when field work was conducted are indicated by vertical gray lines.

**Figure S2** | Vertical profiles of salinity (CTD data) obtained at all transects investigated in March 2017.

**Table S1** | $^{223}$Ra, $^{224}$Raex, $^{228}$Ra activities, and $NO_3^-$ + $NO_2^-$, DSI concentrations determined in water samples collected offshore Côte Bleue, in nearshore waters, in porewaters and in springs (in April 2016, October 2016, December 2016, and March 2017). Shallow porewater and coastal water sample locations are presented in **Figure 2**; seawater sample locations are presented in **Figure 1**. The vertical profiles built in March 2017 are in bold. Porewater samples are underlined in gray.

## REFERENCES

Aminot, A., and Kérouel, R. (2004). *Hydrologie des Écosystèmes Marins: Paramètres et Analyses, Edn.* (Ifremer), 336.

Anderson, M. P. (2005). Heat as a ground water tracer. *Ground Water* 43, 951–968. doi: 10.1111/j.1745-6584.2005.00052.x

Arfib, B., Cavalera, T., and Gilli, E. (2006). Influence de l'hydrodynamique sur l'intrusion saline en aquifère karstique côtier. *Comptes Rendus Geosci.* 338, 757–767. doi: 10.1016/j.crte.2006.07.001

Arfib, B., and Charlier, J.-B. (2016). Insights into saline intrusion and freshwater resources in coastal karstic aquifers using a lumped Rainfall-Discharge-Salinity model (the Port-Miou brackish spring, SE France). *J. Hydrol.* 540, 148–161. doi: 10.1016/j.jhydrol.2016.06.010

Baudron, P., Cockenpot, S., Lopez-Castejon, F., Radakovitch, O., Gilabert, J., Mayer, A., et al. (2015). Combining radon, short-lived radium isotopes and hydrodynamic modeling to assess submarine groundwater discharge from an anthropized semiarid watershed to a Mediterranean lagoon (Mar Menor, SE Spain). *J. Hydrol.* 525, 55–71. doi: 10.1016/j.jhydrol.2015.03.015

Bejannin, S., van Beek, P., Stieglitz, T., Souhaut, M., and Tamborski, J. (2017). Combining airborne thermal infrared images and radium isotopes to study submarine groundwater discharge along the French Mediterranean coastline. *J. Hydrol. Reg. Stud.* 13, 72–90. doi: 10.1016/j.ejrh.2017.08.001

Beusen, A. H. W., Slomp, C. P., and Bouwman, A. F. (2013). Global land-ocean linkage: direct inputs of nitrogen to coastal waters via submarine groundwater discharge. *Environ. Res. Lett.* 8:034035. doi: 10.1088/1748-9326/8/3/034035

Burnett, W. C., Aggarwal, P. K., Aureli, A., Bokuniewicz, H., Cable, J. E., Charette, M. A., et al. (2006). Quantifying submarine groundwater discharge in the coastal zone via multiple methods. *Sci. Total Environ.* 367, 498–543. doi: 10.1016/j.scitotenv.2006.05.009

Burnett, W. C., Bokuniewicz, H., Huettel, M., Moore, W. S., and Taniguchi, M. (2003). Groundwater and pore water inputs to the coastal zone. *Biogeochemistry* 66, 3–33. doi: 10.1023/B: BIOG.0000006066.21240.53

Burnett, W. C., and Dulaiova, H. (2003). Estimating the dynamics of groundwater input into the coastal zone via continuous radon-222 measurements. *J. Environ. Radioactiv.* 69, 21–35. doi: 10.1016/S0265-931X(03)00084-5

Charette, M. A., and Buesseler, K. O. (2004). Submarine groundwater discharge of nutrients and copper to an urban subestuary of Chesapeake Bay (Elizabeth River). *Limnol. Oceanogr.* 49, 376–385. doi: 10.4319/lo.2004.49.2.0376

Charette, M. A., Buesseler, K. O., and Andrews, J. E. (2001). Utility of radium isotopes for evaluating the input and transport of groundwater-derived

nitrogen to a Cape Cod estuary. *Limnol. Oceanogr.* 46, 465–470. doi: 10.4319/lo.2001.46.2.0465

Charette, M. A., Gonneea, M. E., Morris, P. J., Statham, P., Fones, G., Planquette, H., et al. (2007). Radium isotopes as tracers of iron sources fueling a Southern Ocean phytoplankton bloom. *Deep Sea Res. Part II Top. Stud. Oceanogr.* 54, 1989–1998. doi: 10.1016/j.dsr2.2007.06.003

Cockenpot, S., Claude, C., and Radakovitch, O. (2015). Estimation of air water gas exchange coefficient in a shallow lagoon based on $^{222}$Rn mass balance. *J. Environ. Radioactiv.* 143, 58–69. doi: 10.1016/j.jenvrad.2015.02.007

Dulaiova, H., Peterson, R., Burnett, W. C., and Lane-Smith, D. (2005). A multi-detector continuous monitor for assessment of $^{222}$Rn in the coastal ocean. *J. Radioanal. Nucl. Chem.* 263, 361–363. doi: 10.1007/s10967-005-0063-8

Estournel, C., Broche, P., Marsaleix, P., Devenon, J. L., Auclair, F., and Vehil, R. (2001). The Rhone River plume in unsteady conditions: numerical and experimental results. *Estuarine Coast. Shelf. Sci.* 53, 25–38. doi: 10.1006/ecss.2000.0685

Estournel, C., Durrieu de Madron, X., Marsaleix, P., Auclair, F., Julliand, C., and Vehil, R. (2003). Observation and modelisation of the winter coastal oceanic circulation in the Gulf of Lions under wind conditions influenced by the continental orography (FETCH experiment). *J. Geophys. Res.* 108:8059. doi: 10.1029/2001JC000825

Fleury, P., Bakalowicz, M., and de Marsily, G. (2007). Submarine springs and coastal karst aquifers: a review. *J. Hydrol.* 339, 79–92. doi: 10.1016/j.jhydrol.2007.03.009

Frayse, M., Pairaud, I., Ross, O. N., Faure, V. M., and Pinazo, C. (2014). Intrusion of Rhone River diluted water into the Bay of Marseille: generation processes and impacts on ecosystem functioning. *J. Geophys. Res. Oceans* 119, 6535–6556. doi: 10.1002/2014JC010022

Garcia-Solsona, E., Garcia-Orellana, J., Masqué, P., and Dulaiova, H. (2008). Uncertainties associated with $^{223}$Ra and $^{224}$Ra measurements in water via a Delayed Coincidence Counter (RaDeCC). *Mar. Chem.* 109, 198–219. doi: 10.1016/j.marchem.2007.11.006

Garcia-Solsona, E., Garcia-Orellana, J., Masqué, P., Rodellas, V., Mejías, M., Ballesteros, B., et al. (2010). Groundwater and nutrient discharge through karstic coastal springs (Castello, Spain). *Biogeosciences* 7, 2625–2638. doi: 10.5194/bg-7-2625-2010

Gatti, J., Petrenko, A., Devenon, J.-L., Leredde, Y., and Ulses, C. (2006). The Rhone river dilution zone present in the northeastern shelf of the Gulf of Lion in December 2003. *Cont. Shelf Res.* 26, 1794–1805, doi: 10.1016/j.csr.2006.05.012

Gobler, C., and Sañudo-Wilhelmy, S. (2001). Effects of organic carbon, organic nitrogen, inorganic nutrients, and iron additions on the growth of
phytoplankton and bacteria during a brown tide bloom. *Mar. Ecol. Prog. Ser.* 209, 19–34. doi: 10.3354/meps209019

Hancock, G. J., Webster, I. T., and Stieglitz, T. C. (2006). Horizontal mixing of Great Barrier Reef waters: offshore diffusivity determined from radium isotope distribution. *J. Geophys. Res.* 111:C12019. doi: 10.1029/2006JC003608

Hwang, D.-W., Kim, G., Lee, Y.-W., and Yang, H.-S. (2005). Estimating submarine inputs of groundwater and nutrients to a coastal bay using radium isotopes. *Mar. Chem.* 96, 61–71. doi: 10.1016/j.marchem.2004.11.002

Jensen, A., Collins, K., and Lockwood, A. P. (2012). *Artificial Reefs in European Seas*. Springer Science and Business Media.

Jeong, J., Kim, G., and Han, S. (2012). Influence of trace element fluxes from submarine groundwater discharge (SGD) on their inventories in coastal waters off volcanic island, Jeju, Korea. *Appl. Geochem.* 27, 37–43. doi: 10.1016/j.apgeochem.2011.08.014

Jouvenel, J.-Y., Bachet, F., Harmenin, J.-G., and Bellan-Santini, D. (2004). "Suivi biologique d'une réserve marine de la Côte Bleue (Golfe de Marseille, Méditerranée, France)," in *Colloque «Importance de La Recherche Dans Les Aires Protégées: Des Fondements à La Gestion», Organisé En Guadeloupe Par La SNPN, 5-7 Juin 2002* (Paris: Société nationale de protection de la nature et d'acclimatation de France).

Knee, K. L., Garcia-Solsona, E., Garcia-Orellana, J., Boehm, A. B., and Paytan, A. (2011). Using radium isotopes to characterize water ages and coastal mixing rates: a sensitivity analysis. *Limnol. Oceanogr. Methods* 9, 380–395. doi: 10.4319/lom.2011.9.380

Knee, K. L., Street, J. H., Grossman, E. E., Boehm, A. B., and Paytan, A. (2010). Nutrient inputs to the coastal ocean from submarine groundwater discharge in a groundwater-dominated system: relation to land use (Kona coast, Hawaii, U.S.A.). *Limnol. Oceanogr.* 55, 1105–1122. doi: 10.4319/lo.2010.55.3.1105

Lecher, A. L., MacKay, K., Kudela, R., Ryan, J., Fisher, A., Murray, J., et al. (2015). Nutrient loading through submarine groundwater discharge and phytoplankton growth in Monterey Bay, CA. *Environ. Sci. Technol.* 49, 6665–6673. doi: 10.1021/acs.est.5b00909

Li, C., and Cai, W.-J. (2011). On the calculation of eddy diffusivity in the shelf water from radium isotopes: high sensitivity to advection. *J. Mar. Syst.* 86, 28–33. doi: 10.1016/j.jmarsys.2011.01.003

Ludwig, W., Dumont, E., Meybeck, M., and Heussner, S. (2009). River discharges of water and nutrients to the Mediterranean and Black Sea: major drivers for ecosystem changes during past and future decades? *Prog. Oceanogr.* 80, 199–217. doi: 10.1016/j.poccean.2009.02.001

Marsaleix, P., Auclair, F., Floor, J. W., Herrmann, M. J., Estournel, C., Pairaud, I., et al. (2008). Energy conservation issues in sigma-coordinate free-surface ocean models. *Ocean Model.* 20, 61–89. doi: 10.1016/j.ocemod.2007.07.005

Marsaleix, P., Michaud, H., and Estournel, C. (2019). 3D phase-resolved wave modelling with a non-hydrostatic ocean circulation model. *Ocean Model.* 136, 28–50. doi: 10.1016/j.ocemod.2019.02.002

Mejías, M., Ballesteros, B. J., Antón-Pacheco, C., Domínguez, J. A., Garcia-Orellana, J., and Garcia-Solsona, E., et al. (2012). Methodological study of submarine groundwater discharge from a karstic aquifer in the Western Mediterranean Sea. *J. Hydrol.* 464–465, 27–40. doi: 10.1016/j.jhydrol.2012.06.020

Millot, C. (1990). The Gulf of Lions' hydrodynamics. *Cont. Shelf Res.* 10, 885–894. doi: 10.1016/0278-4343(90)90065-T

Montiel, D., Dimova, N., Andreo, B., Prioeto, J., Garcia-Orellana, J., and Rodellas, V. (2018). Assessing submarine groundwater discharge (SGD) and nitrate fluxes in highly heterogeneous coastal karst aquifers: challenges and solutions. *J. Hydrol.* 557, 222–242. doi: 10.1016/j.jhydrol.2017.12.036

Moore, W. S. (1996). Using the radium quartet for evaluating groundwater input and water exchange in salt marshes. *Geochim. Cosmochim. Acta* 60, 4645–4652. doi: 10.1016/S0016-7037(96)00289-X

Moore, W. S. (1999). The subterranean estuary: a reaction zone of ground water and sea water. *Mar. Chem.* 65, 111–125. doi: 10.1016/S0304-4203(99)00014-6

Moore, W. S. (2000a). Determining coastal mixing rates using radium isotopes. *Cont. Shelf Res.* 20, 1993–2007. doi: 10.1016/S0278-4343(00)00054-6

Moore, W. S. (2000b). Ages of continental shelf waters determined from $^{223}$Ra and $^{224}$Ra. *J. Geophys. Res. Oceans* 105, 22117–22122. doi: 10.1029/1999JC000289

Moore, W. S. (2008). Fifteen years experience in measuring $^{224}$Ra and $^{223}$Ra by delayed-coincidence counting. *Mar. Chem.* 109, 188–197. doi: 10.1016/j.marchem.2007.06.015

Moore, W. S. (2010). The effect of submarine groundwater discharge on the ocean. *Annu. Rev. Mar. Sci.* 2, 59–88. doi: 10.1146/annurev-marine-120308-081019

Moore, W. S., and Arnold, R. (1996). Measurement of $^{223}$Ra and $^{224}$Ra on coastal waters using a delayed coincidence counter. *J. Geophys. Res. Oceans* 101, 1321–1329. doi: 10.1029/95JC03139

Moore, W. S., and Reid, D. F. (1973). Extraction of radium from natural waters using manganese-impregnated acrylic fibers. *J. Geophys. Res.* 78, 8880–8886. doi: 10.1029/JC078i036p08880

Niencheski, L. F. H., Windom, H. L., Moore, W. S., and Jahnke, R. A. (2007). Submarine groundwater discharge of nutrients to the ocean along a coastal lagoon barrier, Southern Brazil. *Mar. Chem.* 106, 546–561. doi: 10.1016/j.marchem.2007.06.004

Ollivier, P., Claude, C., Radakovitch, O., and Hamelin, B. (2008). TIMS measurements of 226Ra and 228Ra in the Gulf of Lion, an attempt to quantify submarine groundwater discharge. *Mar. Chem.* 109, 337–354. doi: 10.1016/j.marchem.2007.08.006

Pairaud, I. L., Gatti, J., Bensoussan, N., Verney, R., and Garreau, P. (2011). Hydrology and circulation in a coastal area off Marseille: validation of a nested 3D model with observations. *J. Mar. Syst.* 88, 20–33. doi: 10.1016/j.jmarsys.2011.02.010

Pavlidou, A., Papadopoulos, V. P., Hatzianestis, I., Simboura, N., Patiris, D., and Tsabaris, C. (2014). Chemical inputs from a karstic submarine groundwater discharge (SGD) into an oligotrophic Mediterranean coastal area. *Sci. Total Environ.* 488–489, 1–13. doi: 10.1016/j.scitotenv.2014.04.056

Paytan, A., Shellenbarger, G. G., Street, J. H., Gonneea, M. E., Davis, K., and Young, M. B., et al. (2006). Submarine groundwater discharge: an important source of new inorganic nitrogen to coral reef ecosystems. *Limnol. Oceanogr.* 51, 343–348. doi: 10.4319/lo.2006.51.1.0343

Petrenko, A. (2003). Variability of circulation features in the Gulf of Lion NW Mediterranean Sea. Importance of inertial currents. *Oceanol. Acta* 26, 323–338. doi: 10.1016/S0399-1784(03)00038-0

Petrenko, A., Dufau, C., and Estournel, C. (2008). Barotropic eastward currents in the western Gulf of Lion, north-western Mediterranean Sea, during stratified conditions. *J. Mar. Syst.* 74, 406–428. doi: 10.1016/j.jmarsys.2008.03.004

Pont, D., Simonnet, J.-P., and Walter, A. V. (2002), Medium-term changes in suspended sediment delivery to the ocean: consequences of catchment heterogeneity and river management (Rhône River, France). *Estuarine Coastal Shelf Sci.* 54, 1–18, doi: 10.1006/ecss.2001.0829

Reffray, G., Fraunie, P., and Marsaleix, P. (2004). Secondary flows induced by wind forcing in the Rhône region of freshwater influence. *Ocean Dyn.* 54, 179–196. doi: 10.1007/s10236-003-0079-y

Rodellas, V., Garcia-Orellana, J., Masqué, P., Feldman, M., and Weinstein, Y. (2015). Submarine groundwater discharge as a major source of nutrients to the Mediterranean Sea. *Proc. Natl. Acad. Sci. U.S.A.* 112, 3926–3930. doi: 10.1073/pnas.1419049112

Rodellas, V., Garcia-Orellana, J., Trezzi, G., Masqué, P., Stieglitz, T. C., Bokuniewicz, H., et al. (2017). Using the radium quartet to quantify submarine groundwater discharge and porewater exchange. *Geochim. Cosmochim. Acta* 196, 58–73. doi: 10.1016/j.gca.2016.09.016

Schmidt, S., and Reys, J.-L. (1996). Radium as internal tracer of Mediterranean Outflow Water. *J. Geophys. Res.* 101, 3589–3596.

Slomp, C. P., and Van Cappellen, P. (2004). Nutrient inputs to the coastal ocean through submarine groundwater discharge: controls and potential impact. *J. Hydrol.* 295, 64–86. doi: 10.1016/j.jhydrol.2004.02.018

Stieglitz, T., van Beek, P., Souhaut, M., and Cook, P. (2013). Groundwater discharge Karstic groundwater discharge and seawater recirculation through sediments in shallow coastal Mediterranean lagoons, determined from water, salt and radon budgets. *Mar. Chem.* 156, 73–84. doi: 10.1016/j.marchem.2013.05.005

Stieglitz, T. C., Cook, P. G., and Burnett, W. C. (2010). Inferring coastal processes from regional-scale mapping of $^{222}$Radon and salinity: examples from the Great Barrier Reef, Australia. *J. Environ. Radioact.* 101, 544–552. doi: 10.1016/j.jenvrad.2009.11.012

Sun, Y., and Torgersen, T. (1998). The effects of water content and Mn-fiber surface conditions on $^{224}$Ra measurement by $^{220}$Rn emanation. *Mar. Chem.* 62, 299–306. doi: 10.1016/S0304-4203(98)00019-X

Swarzenski, P. W. (2007). U/Th series radionuclides as coastal groundwater tracers. *Chem. Rev.* 107, 663–674. doi: 10.1021/cr0503761
Tamborski, J., Bejannin, S., Garcia-Orellana, J., Souhaut, M., Charbonnier, C., Anschutz, P., et al. (2018). A comparison between water circulation and terrestrially-driven dissolved silica fluxes to the Mediterranean Sea traced using radium isotopes. *Geochim. Cosmochim. Acta* 238, 496–515. doi: 10.1016/j.gca.2018.07.022

Tamborski, J. J., Cochran, J. K., and Bokuniewicz, H. J. (2017). Submarine groundwater discharge driven nitrogen fluxes to Long Island Sound, NY: terrestrial vs. marine sources. *Geochim. Cosmochim. Acta* 218, 40–57. doi: 10.1016/j.gca.2017.09.003

Tamborski, J. J., Rogers, A. D., Bokuniewicz, H. J., Cochran, J. K., and Young, C. R. (2015). Identification and quantification of diffuse fresh submarine groundwater discharge via airborne thermal infrared remote sensing. *Remote Sens. Environ.* 171, 202–217. doi: 10.1016/j.rse.2015.10.010

Tovar-Sánchez, A., Basterretxea, G., Rodellas, V., Sánchez-Quiles, D., Garcia-Orellana, J., Masqué, P., et al. (2014). Contribution of groundwater discharge to the coastal dissolved nutrients and trace metal concentrations in Majorca Island: Karstic vs Detrital systems. *Environ. Sci. Technol.* 48, 11819–11827. doi: 10.1021/es502958t

Trezzi, G., Garcia-Orellana, J., Rodellas, V., Masqué, P., Garcia-Solsona, E., and Andersson, P. S. (2017). Assessing the role of submarine groundwater discharge as a source of Sr to the Mediterranean Sea. *Geochim. Cosmochim. Acta* 200, 42–54. doi: 10.1016/j.gca.2016.12.005

van Beek, P., Souhaut, M., Lansard, B., Bourquin, M., Reyss, J.-L., von Ballmoos, P., et al. (2013). LAFARA: a new underground laboratory in the French Pyrénées for ultra low-level gamma-ray spectrometry. *J. Environ. Radioact.* 116, 152–158. doi: 10.1016/j.jenvrad.2012.10.002

van Beek, P., Souhaut, M., and Reyss, J.-L. (2010). Measuring the radium quartet (²²⁸Ra, ²²⁶Ra, ²²⁴Ra, ²²³Ra) in seawater samples using gamma spectrometry. *J. Environ. Radioact.* 101, 521–529. doi: 10.1016/j.jenvrad.2009.12.002

van Beek, P., Sternberg, E., Reyss, J.-L., Souhaut, M., Robin, E., and Jeandel, C. (2009). ²²⁸Ra/²²⁶Ra and ²²⁶Ra/Ba ratios in the Western Mediterranean Sea: barite formation and transport in the water column. *Geochim. Cosmochim. Acta* 73, 4720–4737 doi: 10.1016/j.gca.2009.05.063

Webster, I. T., Hancock, G. J., and Murray, A. S. (1995). Modelling the effect of salinity on radium desorption from sediments. *Geochim. Cosmochim. Acta* 59, 2469–2476. doi: 10.1016/0016-7037(95)00141-7

Weinstein, Y., Yechieli, Y., Shalem, Y., Burnett, W. C., Swarzenski, P. W., and Herut, B. (2011). What is the role of fresh groundwater and recirculated seawater in conveying nutrients to the Coastal Ocean? *Environ. Sci. Technol.* 45, 5195–5200. doi: 10.1021/es104394r

Windom, H. L., Moore, W. S., Niencheski, L. F. H., and Jahnke, R. A. (2006). Submarine groundwater discharge: a large, previously unrecognized source of dissolved iron to the South Atlantic Ocean. *Mar. Chem.* 102, 252–266. doi: 10.1016/j.marchem.2006.06.016

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Bejannin, Tamborski, van Beek, Souhaut, Stieglitz, Radakovitch, Claude, Conan, Pujo-Pay, Crispi, Le Roy and Estournel. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.