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

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 (223Ra, 224Ra, 228Ra) 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 223Ra and 224Ra were used to derive horizontal eddy diffusivity coefficients that were subsequently combined with surface water nutrient gradients (NO2− + NO3−, DSi) to determine the net nutrient fluxes from SGD. The estimated DSi and (NO2− + NO3−) fluxes are 6.2 ± 5.0 *103 and 4.0 ± 2.0 *103 mol d−1 per km of coastline, respectively. We attempted to further constrain these SGD fluxes by combining horizontal eddy diffusivity and 228Ra 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 228Ra (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 228Ra 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 scale, DSi fluxes driven by SGD vary between 0.1 and 1.4% of the DSi inputs of the Rhône River, while the (NO2− + NO3−) 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).

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 3 and 4.0 ± 2.0 * 10 3 mol d −1 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 halflife 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 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 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., 2015Rodellas et al., , 2017Tamborski et al., 2018). In particular, little information is available on the nutrient fluxes associated with SGD along the French Mediterranean coastline 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 SGDdriven nutrient estimates based on a basin-wide 228 Ra mass balance (Rodellas et al., 2015). The 228 Ra inventory, after considering known 228 Ra sources and 228 Ra sinks, is divided by the 228 Ra submarine groundwater endmember activity to determine a SGD flux. Importantly, Rodellas et al. (2015) note that 228 Ra SGD endmembers reported throughout the Mediterranean Basin span a wide-range, from 10 2 to 10 5 dpm m −3 . 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 h ), that were subsequently combined with surface water nutrient gradients (DSi, NO 2 − + NO 3 − ) in order to determine the net nutrient flux from SGD. Offshore transects of 228 Ra were also investigated with the aim to derive SGDdriven nutrient fluxes by combining K h and 228 Ra gradients with the SGD endmember. Fluxes were investigated over a 1-year period (April 2016, October 2016, December 2016, and March 2017

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.
Frontiers in Environmental Science | www.frontiersin.org 3 February 2020 | Volume 7 | Article 205 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.
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 −1 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 222 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 WTW TM 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 2coated acrylic fibers ("Mn-fiber") at a flow rate <1 L min −1 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 223 Ra and 224 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 Ra ex ). 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. 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 semiplanar 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 HgCl 2 (10 µg L −1 ) and stored at 4 • C in the dark. In the laboratory, nitrate (NO 3 − ), nitrite (NO 2 − ), phosphate (PO 3− 4 ) and dissolved silica (DSi) concentrations were measured on a continuous flow autoanalyzer Technicon R AutoAnalyzer II at LOMIC, Banyuls-sur-Mer (Aminot and Kérouel, 2004). The analytical precision of NO 3 − , NO 2 − , PO 3− 4 , 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 1/2 = 11.4 d; 224 Ra, t 1/2 = 3.66 d; 226 Ra, t 1/2 = 1,600 y; 228 Ra, t 1/2 = 5.75 y), as well as 222 Rn (t 1/2 = 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 halflives (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): where A is the radium isotope activity, K h is the horizontal eddy diffusivity coefficient (m 2 s −1 ), d is the distance offshore and λ is the decay constant of the Ra isotope used. Assuming steady state, the solution of Equation (1) is: 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 Ra ex ) 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 = λ/m² where m is the slope of the linear regression). Horizontal eddy diffusivity coefficients (m 2 s −1 ) 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 Frontiers in Environmental Science | www.frontiersin.org 5 February 2020 | Volume 7 | Article 205 measurement (Garcia-Solsona et al., 2008). The uncertainty of K h is thus: 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 −1 km −1 ) is defined as the slope of the plot of a nutrient concentration as a function of offshore distance. A nutrient flux (µmol s −1 km −1 ) 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 −1 km −1 ), 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(Marsaleix et al., , 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.

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 Carryle-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 −1 ) than other in situ measurements taken along the coast (<370 dpm 100L −1 ), 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 Ra ex 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 Frontiers in Environmental Science | www.frontiersin.org 6 February 2020 | Volume 7 | Article 205  in activity, with a maximum value in October (4.8 dpm 100 L −1 ) and a minimum value in March (1.0 dpm 100 L −1 ). 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 222 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 224 Ra ex activities ranging from below detection limit (i.e., no excess 224 Ra) to 3.6 dpm 100 L −1 , with an average (±STD) value of 1.3 (±0.9) dpm 100 L −1 (n = 16; Table 1). The 228 Ra activities in the coastal seas (<8 km) range between 1.7 and 24 dpm 100 L −1 , with a mean value of 4.2 ± 2.3 dpm 100 L −1 (n = 98;  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)  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 228 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). 228 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 −1 , 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). 223 Ra and 224 Ra ex activities were between 2-39 and 21-324 dpm 100 L −1 in porewaters while 228 Ra activities were 2.1-203 dpm 100 L −1 (Figure 5). Finally, three karstic springs that discharge in Laurons bay at transect T1 showed average 223 Ra, 224 Ra ex , and 228 Ra activities of 37 (±25), 1,485 (±658), and 826 (±45) dpm 100 L −1 , respectively (n = 3; Table 1). Taken together, these lines of evidence suggest that the chemical enrichments observed along this 22 km coastline ( 223 Ra, 224 Ra ex , 228 Ra, and 222 Rn; Table 1) are driven by SGD. This includes the discharge of terrestrial groundwater (e.g., via 222 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).
Porewater concentrations were approximately one order of magnitude higher than coastal surface waters, with values between 0.01 and 4.0 µM for PO 3− 4 , 1.9 and 133 µM for DSi, and 22 and 194 µM for NO 3 − 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 3 − 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 3− 4 was not analyzed for these springs.

Determination of the Horizontal Eddy Diffusivity Coefficients K h
All of the surface water transects were used to estimate horizontal eddy diffusivity coefficients, K h (Figure 7). We do not observe any significant difference between the 223 Ra and 224 Ra ex derived Frontiers in Environmental Science | www.frontiersin.org 8 February 2020 | Volume 7 | Article 205 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 2 s −1 ( Table 3). In October 2016, 5 transects were investigated and the mean K h is 96 (±44) m 2 s −1 . The highest K h was estimated for December 2016, with a K h of 184 (±112) m 2 s −1 , 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 2 s −1 , 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.
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 3 − fluxes are similar between the different sampling periods (Figure 10; Table 3). NO 3 − fluxes estimated for transect 4 were similar in April and October. The average (±STD) DSi and NO 3 − fluxes, for all the transects in which a nutrient gradient was estimated, equal to 6.2 (±5.0) * 10 3 and 4.0 (±2.0) * 10 3 mol d −1 km −1 , 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 224 Ra/ 228 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; 228 Ra offshore gradients), in addition to the nutrient gradient method, to constrain SGD along Côte Bleue.
228 Ra Gradient Method The regional view of the entire dataset is that the 228 Ra activities do not decrease with increasing distance offshore (Figure 12). Only two transects among the 16 transects investigated here display statistically significant offshore 228 Ra gradients (p < 0.01) and we thus derived SGD fluxes from these data only. The 228 Ra gradient for T7 in December 2016 is −0.164 (±0.065) dpm 100 L −1 km −1 while the 228 Ra gradient for T5 in March 2017 is −0.225 (±0.058) dpm 100 L −1 km −1 . The 228 Ra fluxes are estimated by combining these gradients with the 224 Ra-derived Frontiers in Environmental Science | www.frontiersin.org 10 February 2020 | Volume 7 | Article 205 FIGURE 6 | (A) DSi (µM), (B) NO 2 − +NO 3 − (µM), and (C) PO 3− 4 (µ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.
horizontal eddy diffusivity coefficient corresponding to the appropriate month (section Determination of the Horizontal Eddy Diffusivity Coefficients K h ). The 228 Ra fluxes are thus 2.6 (±1.9) * 10 8 and 1.8 (±0.8) * 10 8 dpm d −1 km −1 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 228 Ra fluxes by the 228 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 ( 228 Ra of 826 ± 45 dpm 100 L −1 ; n = 3), this yields a volumetric SGD flux for T7 in December 2016 of 1.6 (±1.2)−3.2 (±2.3) * 10 4 m 3 d −1 km −1 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 4 m 3 d −1 km −1 , 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 3 mol d −1 km −1 for T7 in December and 0.8 (±0.6)−1.6 (±1.2) * 10 3 mol d −1 km −1 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 h * Nutrient Offshore Gradient (Figure 10). Using the NO 3 − concentrations of the brackish springs (2.4 ± 2.6 µM; two springs considered out of the three), the NO 3 − fluxes estimated from the 228 Ra gradient are 0.04 (±0.06)−0.08 (±0.10) * 10 3 and 0.02 (±0.03)−0.03 (±0.04) * 10 3 mol d −1 km −1 for T7 and T5, respectively. The estimations of NO 3 − fluxes are lower than the estimations made in section Mixing Model: K h * Nutrient Offshore Gradient. Note that the NO 3 − 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 3 − fluxes being lower than the estimations made in section Mixing Model: K h * Nutrient Frontiers in Environmental Science | www.frontiersin.org 11 February 2020 | Volume 7 | Article 205  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 the importance of properly selecting the SGD endmember in determining nutrient fluxes from offshore 228 Ra gradients. We should be able to estimate an SGD-driven PO 3− 4 flux using this method, whereas we could not from the method used in section Mixing Model: K h * Nutrient Offshore Gradient because there was no observable surface water PO 3− 4 gradient. Considering the maximum PO 3− 4 concentrations of the porewater samples (4.04 µM; Table 2), PO 3− 4 fluxes are 0.8 (±0.7)−1.6 (±1.4) * 10 2 and 0.4 (±0.3)−0.7 (±0.6) * 10 2 mol d −1 km −1 for T7 and T5, respectively. PO 3− 4 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 2 mol d −1 km −1 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 h * Nutrient Offshore Gradient), particularly because knowledge of the SGD endmember is not requisite for this method. In addition, the 228 Ra distribution during the different campaigns did not show a decreasing trend with increasing distance offshore, which prevented us from using the 228 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 3 s −1 ; 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 228 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 228 Ra activities are especially high at T1 near the coastline (∼25 dpm 100 L −1 ; 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 3 − concentrations during this period (T1-T4; Figure 9). As discussed above (section Mixing Model: K h * Nutrient Offshore Gradient), the surface salinity maps generated using the SYMPHONIE model indicate intrusion 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 228 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. 224 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 228 Ra gradient. We estimated nutrient fluxes supplied by local rivers using river gauging stations (hydro.eaufrance.fr) and monthly chemical elements 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  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 3 − 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 222 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 3 − fluxes driven by the Huveaune River, with significant temporal variability ( Table 4) 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 longlived 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  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. DSi and NO 3 − fluxes reported here were estimated using Ra isotopes (K h * nutrient gradients) considering 10 m as the water depth impacted by SGD. Note than the Rhône River fluxes are reported in 10 6 mol d −1 .
supplies ≤1% of the DSi and NO 3 − 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.