Assessing land-ocean connectivity via submarine groundwater discharge (SGD) in the Ria Formosa Lagoon (Portugal): combining radon measurements and stable isotope hydrology

Abstract. Natural radioactive tracer-based assessments of basin-scale submarine groundwater discharge (SGD) are well developed. However, SGD takes place in different modes and the flow and discharge mechanisms involved occur over a wide range of spatial and temporal scales. Quantifying SGD while discriminating its source functions therefore remains a major challenge. However, correctly identifying both the fluid source and composition is critical. When multiple sources of the tracer of interest are present, failure to adequately discriminate between them leads to inaccurate attribution and the resulting uncertainties will affect the reliability of SGD solute loading estimates. This lack of reliability then extends to the closure of local biogeochemical budgets, confusing measures aiming to mitigate pollution. Here, we report a multi-tracer study to identify the sources of SGD, distinguish its component parts and elucidate the mechanisms of their dispersion throughout the Ria Formosa – a seasonally hypersaline lagoon in Portugal. We combine radon budgets that determine the total SGD (meteoric + recirculated seawater) in the system with stable isotopes in water (δ2H, δ18O), to specifically identify SGD source functions and characterize active hydrological pathways in the catchment. Using this approach, SGD in the Ria Formosa could be separated into two modes, a net meteoric water input and another involving no net water transfer, i.e., originating in lagoon water re-circulated through permeable sediments. The former SGD mode is present occasionally on a multi-annual timescale, while the latter is a dominant feature of the system. In the absence of meteoric SGD inputs, seawater recirculation through beach sediments occurs at a rate of  ∼  1.4  ×  106 m3 day−1. This implies that the entire tidal-averaged volume of the lagoon is filtered through local sandy sediments within 100 days ( ∼  3.5 times a year), driving an estimated nitrogen (N) load of  ∼  350 Ton N yr−1 into the system as NO3−. Land-borne SGD could add a further  ∼  61 Ton N yr−1 to the lagoon. The former source is autochthonous, continuous and responsible for a large fraction (59 %) of the estimated total N inputs into the system via non-point sources, while the latter is an occasional allochthonous source capable of driving new production in the system.


Introduction
Freshwater inputs into the coastal zone are important pathways for the transfer of land-borne solutes and particulates into the sea.Even if channeled freshwater flows such as rivers are relatively well gauged worldwide, sub-surface sources are more difficult to quantify in coastal settings.This dif-Published by Copernicus Publications on behalf of the European Geosciences Union.
ficulty has hindered the understanding of current drivers of coastal ecosystem decline (Carpenter et al., 1998;Finkl and Krupa, 2003).Indeed, on a global scale, an estimated 6 % of the freshwater input into the sea, carrying an anticipated 52 % of the total dissolved salts crossing the land-ocean interface, was estimated to occur via SGD -submarine groundwater discharge -by Zektser and Loaiciga (1993).This early estimate has since been updated by Kwon et al. (2014), who show that global SGD is 3-4 times greater than the freshwater flow into the oceans by rivers.This revision means that SGD is by far the largest contributor of terrestrial solutes to the global ocean, hence implying that some global biogeochemical budgets of major elements need revision.However, mass flows defining the contribution of SGD to coastal biogeochemical budgets are difficult to quantify in a systematic way (Burnett et al., 2001a).
To understand the contribution of groundwater-seawater interactions to marine biogeochemistry (Moore, 1996(Moore, , 2006;;Moore and Church, 1996;Church, 1996), the definition of SGD encompasses any flow of water across the sea floor, regardless of fluid composition or driving force (Burnett et al., 2003).This is because reactivity of solutes when meteoric water and seawater mix and travel through porous media significantly alters the composition of the discharging water with respect to both original contributions (Moore, 1999(Moore, , 2010)).Submarine groundwater discharge is therefore not limited to fresh groundwater discharge, but includes seawater recirculation through coastal sediments (Li et al., 1999) and seasonal repositioning of the saltwater-freshwater interface (Michael et al., 2005;Edmunds, 2003;Santos et al., 2009).All of these promote changes to the rates of transfer, mixing and chemical reaction at the subterranean estuary (Moore, 1999;Charette et al., 2005;Charette and Sholkovitz, 2006;Robinson et al., 2007), altering the original chemical signatures in a non-uniform way at system scale (Slomp and van Cappellen, 2004;Spiteri et al., 2008).
Tracer-based assessments of basin-scale SGD are well developed (Burnett et al., 2001a(Burnett et al., , b, 2003(Burnett et al., , 2008)), but because the flow and discharge mechanisms involved cover a wide range of spatial and temporal scales (Bratton, 2010;Santos et al., 2012), quantifying SGD while discriminating its source functions is still a challenge (e.g., Mulligan and Charette, 2006).Indeed, the most common approaches to estimating SGD are (a) radioactive tracer studies specifically looking at radon ( 222 Rn, T 1/2 = 3.8 days; Burnett et al., 2001a, b) and radium isotopes (Moore and Arnold, 1996); (b) direct measurement of discharge fluxes over small areas (Lee, 1977;Michael et al., 2003;Taniguchi et al., 2003); and (c) modeling.Direct measurements offer limited spatial coverage and are labor intensive (e.g., Leote et al., 2008), making reliable flux estimates at the system scale difficult.Modeling approaches depend on the water and/or salt budgets, hydrograph separation techniques, or descriptions of interfacial flow dynamics based on Darcy's law.Frequently, however, they incorporate assumptions of a steady state inventory and homogeneity of hydraulic conductivity over large-scale lengths and fail to include seawater recirculation.In addition, there is often a mismatch between the spatial and/or temporal scale of the model outputs and those necessary to close coastal biogeochemical budgets (Prieto and Destouni, 2010).
Radioactive tracer studies produce spatially integrated estimates of flux (Cable et al., 1996;Moore, 1996) while simultaneously dampening the effects of short-term variability (Burnett et al., 2001a).However, while radon budgets produce an estimate of "total" SGD, i.e., freshwater inputs + re-circulated seawater (Mulligan and Charette, 2006), radium budgets primarily assess the salty component of SGD given that radium is normally absent in fresh groundwater but might be mobilized from sediment particles in case of saline water influence (Webster et al., 1995).Even so, the variety of ubiquitous temporally and spatially variable sediment-water exchange mechanisms that also act as sources of radon (Cable et al., 2004;Martin et al., 2004;Colbert et al., 2008a, b) and short-lived radium isotopes to surface waters (Webster et al., 1994;Hancock and Murray, 1996;Hancock et al., 2000;Colbert andHammond, 2007, 2008;Gonneea et al., 2008) cannot be ignored.Correctly identifying both the fluid source and composition is thus an important task (Mulligan and Charette, 2006;Burnett et al., 2006).When multiple tracer sources of interest are present, failure to adequately discriminate between them will lead to inaccurate attribution and the resulting uncertainties will affect the reliability of SGD solute loading estimates.
Indeed, as noted by Beck et al. (2007), SGD-borne chemical load into coastal systems is usually predicted by combining measurements of source composition with SGD estimates.Linking these two data sets requires care and is underpinned by our ability to correctly identify and quantify the different SGD pathways into any one system.This is because the final SGD solute-load estimate not only depends on how accurate our recognition of the SGD source functions is, but also on the ability to track their path within the system, since this is required to evaluate the biogeochemical history of the source components prior to their mixture into receiving waters.Not fulfilling this requisite therefore constitutes the major obstacle to prognosticating upper boundary or "potential" SGD-related impact, and, more importantly, confidently attributing causality.Indeed, the end-member is usually the greatest source of uncertainty in any tracer or solute mass balance.It follows that determining the end-member concentration in the area(s) most likely to be the source(s) of groundwater would decrease uncertainty in SGD estimates, on the one hand, and in biogeochemical budgets derived from those estimates on the other.The current panorama of SGD research at the system scale therefore begs the question of which end-member to use when selecting a source solute concentration in attempts to quantify pollutant fluxes associated with SGD.
We contribute an answer to this conundrum with a study conducted in a seasonally hypersaline lagoon in southern Hydrol.Earth Syst. Sci., 20, 3077-3098, 2016 www.hydrol-earth-syst-sci.net/20/3077/2016/ Portugal where we combine two data sets: radon surveys are used to determine total SGD in the system, while stable isotopes in water ( 2 H, 18 O) are used to specifically identify SGD sources and characterize active hydrological pathways.We show that, in combination with radon budgeting, stable isotope hydrology is a reliable tool for identifying different SGD sources in a very complex coastal system, even though it has not been used to this end before.This under-use of the methodology has two main reasons.The first is a disciplinary divide: the technique has been the domain of freshwater hydrologists; correlations between δ 18 O and δ 2 H are central to research into the effect of evaporation and mixing on surface waters (Gat et al., 1994;Gibson and Edwards, 2002) and contribute to the disentanglement of different water sources affecting catchments (Rodgers et al., 2005).The other is the paucity of paired δ 18 O-δ 2 H data on coastal sea-water (e.g., Rohling, 2007), even if stable isotope data sets might help constrain the origins of freshwater inputs into the ocean when coupled with salinity data (Munksgaard et al., 2012;Schubert et al., 2015) or as part of a methodological arsenal in SGD studies combining physical and chemical measurements with radioactive and stable isotope tracers (e.g., Povinec et al., 2008).Hence we also bridge the disciplinary gap between marine chemists and hydrogeologists currently extant in SGD studies by using a combined approach merging techniques from both disciplines.The occurrence of SGD comprising significant freshwater contributions was first detected in the Ria Formosa in 2006-2007 and subsequently described as a prominent source of nutrients, in particular nitrogen derived from fertilizers, to the lagoon (Leote et al., 2008;Rocha et al., 2009;Ibánhez et al., 2011Ibánhez et al., , 2013)).However, the unpredictable nature of fresh-water availability in the region, coupled with a mixed-source (i.e., a variable mix of groundwater abstraction and surface water collected in reservoirs) management of public water supply to meet demand (Monteiro and Costa Manuel, 2004;Stigter and Monteiro, 2008), made it unclear whether meteoric groundwater would be a persistent feature of SGD in the system.This made it difficult to clarify the contribution of SGD to the nitrogen budget of the Ria Formosa, with obvious consequences for environmental management strategies.The overarching aims of the study were therefore to identify the sources of SGD, distinguish its component parts and elucidate the mechanisms of their dispersion throughout the Ria Formosa.The outcomes are then employed to distinguish and quantify nitrogen loads carried into the lagoon by different SGD modes.
2 Study site

Geomorphology and hydrodynamics
Located in southern Portugal (36 • 58 N, 8 • 02 W-7 • 03 N, 7 • 32 W), the Ria Formosa (Fig. 1) is a leaky (Kjerfve, 1986) lagoon system separated from the Atlantic by a multi-inlet barrier island cordon.The system covers a surface area of ∼ 111 km 2 and has an average depth of 2 m.The tide is semi-diurnal, with average ranges of 2.8 m for spring tides and 1.3 m for neap tides (Vila-Concejo et al., 2004;Pacheco et al., 2010).The maximum average tidal volume as estimated by the Navy Hydrographical Institute (IH, 1986) is ∼ 140 × 10 6 m 3 .Lagoon water is exchanged with the Atlantic Ocean through six tidal inlets with an average tidal flux of ∼ 8 × 10 6 m 3 (Balouin et al., 2001).Estimates for the submerged area amount to ∼ 55km 2 at high spring tide and between 14 and 22 km 2 at low spring tide (IH, 1986).From west to east (Fig. 1), inlets (Barra, in Portuguese) are identified as Ancão, Faro-Olhão (Barra Nova), Armona (Barra Velha), and Fuzeta, Tavira and Lacem.Barra Nova, Barra Velha and Ancão jointly capture ∼ 90 % of the total tidal prism: 61, 23 and 8 % of the total flow during spring tides and 45, 40 and ∼ 5 % during neap tides, respectively (Pacheco et al., 2010).With the exception of the Barra Nova, all inlets are ebb dominated with residual circulation directed seaward (Dias and Sousa, 2009).

Hydrogeological setting
The regional climate is semi-arid, with an average annual temperature of 17 • C and averages of 11 and 24 • C during winter and summer.The surrounding watershed covers 740 km 2 and receives effective precipitation of 152 mm yr −1 (Salles, 2001), corresponding to an annual rainfall amount of ∼ 1.2 × 10 6 m 3 .There are five minor rivers and 14 streams discharging into the lagoon.Most are ephemeral and dry out during the summer, the exception being the River Gilão, which intermittently discharges almost directly into the At-lantic through the Tavira inlet at the eastern limits of the system.
Three aquifer systems (Fig. 1) border the Ria Formosa (Almeida et al., 2000).These are the Campina de Faro (M12), Chão de Cevada-Quinta João de Ourém (M11) and São João da Venda-Quelfes (M10).The main lithologies supporting these units are Plio-Quaternary, Miocene and Cretaceous formations, comprising, respectively, Pliocene sands and gravels, Quaternary dunes and alluvial deposits; sandy limestones of marine facies; and limestones and detritic limestones.The oldest formation dips to the south, and is found at depths in excess of 200 m near the city of Faro.It is overlain by the Miocene formation extending below the Ria Formosa into the Atlantic Ocean.Sand dunes, sands and gravels of the Plio-Quaternary cover the Miocene and Cretaceous formations within the coastal area.The Campina de Faro (M12, Fig. 1, 86.4 km 2 ) comprises a superficial unconfined aquifer (Pleistocene deposits) with a maximum thickness of 30 m and an underlying Miocene confined multi-layered aquifer, which Engelen and van Beers (1986) suggest discharges directly into the Atlantic Ocean, bypassing the lagoon.The unconfined Pleistocene aquifer is hydraulically connected to the underlying Miocene aquifer.The São João da Venda-Quelfes aquifer (M10, Fig. 1, 113 km 2 ) includes a surface 75 m thick layer of Wealdian facies and an underlying Cretaceous layer of loamy limestone.It contacts with the M12 (Campina de Faro) aquifer and the M11 (Chão de Cevada-Quinta João de Ourém) to the south, and the main flow direction on the eastern side is towards the southeast.Groundwater flow is divergent toward the southeast and the southwest from a central point (Almeida et al., 2000).

Lagoon radon inventory during ebb and flood
Water radon ( 222 Rn) content was measured continuously in situ using two electronic Durridge RAD-7 radon-in-air monitors deployed in tandem on a moving rubber boat during winter (December 2009) and spring (May 2010).Each monitor was coupled to an air-water equilibrator (Durridge RAD-Aqua Accessory) via its own air loop.Non-cavitating centrifugal pumps were used to flush water from ∼ 50 cm below the water surface directly into the equilibrators, at a flow rate of 1.8-2.5 L min −1 .HOBO ™ temperature sensors and a CTD diver (Schlumberger ™ ) continuously recorded the temperature in the mixing chambers and the salinity and temperature of the water being pumped.The counting interval was set at 20 min on each RAD-7 monitor, with the two machines staggered by a 10 min period, allowing for simultaneous replication of 20 min integration periods over the route and increased temporal resolution.Full equilibration between the air within the air loop and the pumped seawater was achieved before surveys started.Sampling began near low tide and continued without an interruption for 24 h.The survey path, recorded with an on-board GPS unit, and the timing were designed to cover the main navigable sectors of the whole lagoon at different tidal stages (ebb and flood) within the course of two complete tidal cycles.In-water radon activity was calculated from the temperature and salinity dependant gas-water equilibrium (Schubert et al., 2012).Radon activities obtained this way were then corrected by the local 226 Ra supported activity, to obtain excess (i.e., unsupported) radon activities.For mass balance purposes, the excess radon inventories were calculated by multiplying the unsupported radon activity from the continuous measurements by the local bathymetric depth, and then normalized to mean tidal height (Burnett and Dulaiova, 2003).

Tidal variability of radon activity at fixed locations
Time series of radon activity were obtained synchronously at two fixed locations within the Faro channel (Fig. 1) during June 2010.The locations were chosen in order to gain insight into the exchange of radon between the lagoon and the adjacent coastal zone through the Barra Nova (Fig. 1) and between the inner reaches of the lagoon and the latter via the Faro channel (Quatro Águas, Fig. 1).Radon activity was measured as described previously, with the added deployment of a CTD diver (Schlumberger ™ ) recording depth, salinity and temperature at the bottom of the channel.The Barra Nova tidal cycle data was then used to calculate the net exchange of radon with the adjacent coastal zone through the main inlet, assuming a vertically well-mixed water column.Exchange of radon through the inlet cross section driven by oscillating tidal flow was determined by first calculating the instantaneous directional flux, F Rn ( t), where t is the counting interval, A Rn ( t) the activity of radon integrated across the counting interval and dh/dt the change in tidal height (r.m.s.l., relative mean sea level) occurring over that interval: The total radon flux was obtained for both the flood and ebb periods by integrating the instantaneous directional fluxes calculated for each counting period (Eq. 1) over time.Radon outflow (when fluxes were negative) and inflow (when positive) are hence obtained for each complete semi-tidal period.The difference between successive outflow and inflow periods gives us the net transfer across the channel during a complete tidal cycle.Data for a minimum of three successive complete tidal cycles, giving three different values for net transfer, were used, and the exchange values determined for each cycle were then averaged to obtain the net exchange flux along the channel at each sampling site.

Complementary radon measurements
Measurements of air temperature, wind speed and atmospheric radon activities were taken on land while the lagoon radon survey progressed.Atmospheric evasion losses (radon degassing flux) were calculated as described in Burnett and Dulaiova (2003), using the equations given in Macintyre et al. (1995) and Turner et al. (1996).Sediment-water diffusive fluxes of radon were measured as described in Corbett et al. (1998) in samples (n = 16) collected throughout the lagoon and directly analyzed in the laboratory upon collection.To obtain these samples, undisturbed sediment cores (35 cm length) were collected using polycarbonate core liners (∅ 5.5 cm) in both sub-tidal (n = 8) and intertidal environments (n = 8), with each environment sub-sampled for sandy and muddy sediments in equal proportions.Resulting fluxes from all analyzed cores where then averaged, and the latter value, with its associated uncertainty, used in subsequent mass balance calculations.

SGD flux estimates based on Rn mass balances
Lagoon radon budget under steady state assumptions The advective flux of radon associated with SGD is determined by the closure of a radon budget incorporating all known sources and sinks of radon in the system (Burnett and Dulaiova, 2003).Mass conservation accounting for the change in the inventory of radon was expressed as where I Rn is the radon inventory measured within the Ria Formosa, t the time, Rn diff the radon flux across the sediment-water interface by diffusion, Rn dg the radon degassing flux, i.e., atmospheric evasion, Rn dy the radon decay flux in the lagoon (i.e., the internal sink), Rn exp and Rn imp the exchange fluxes across inlets, seaward (export) and landward (import), respectively, and Rn adv the advective radon flux putatively associated with SGD.Usually, an additional term accounting for the radon influx via river flow is added if the water and particulate flux associated with river discharge is significant.However, the only perennial river in the Ria Formosa is the Gilão, located at the eastern limit of the lagoon.Salinity measured at the estuary mouth was 29.6 (Table S1 in the Supplement), which in combination with its location implied very low if any inputs of freshwater carrying radon into the system, so we neglected the term.
Assuming steady state of all sinks and sources over the lifetime of radon in the system, then where Rn net is the residual radon exchange flux with the ocean.

Mass balance of radon during ebb and flood
Inventories of radon in the lagoon were determined during ebb and flood.Taking the tide as a travelling wave, the change in the inventory of radon as the tide floods and ebbs has to be balanced by all known radon fluxes occurring within the traversed system during the travel period.If we then take the mean tide level (MTL) as a reference, it follows that the Rn adv term may be calculated for different periods: the period (T ) at which the tidal height in the lagoon is below MTL (Rn adv (T < MTL), i.e., the trough of the tidal wave or low tide, and the one when it is above MTL (Rn adv (T > MTL), corresponding to the peak of the wave, or high tide.Assuming constant mean amplitude for the tidal wave, the corresponding mass conservation equations may be written as follows: Rn adv(T <MTL) = Rn adv(T >MTL) = where I f and I e are the flood and ebb inventories of radon in the lagoon, t the period of the wave (∼ 0.5 day) and Rn adv (T < MTL) and Rn adv (T > MTL) the radon advective fluxes associated with each semi-period (trough and peak stages, respectively).The corresponding continuity equation, describing the net advective flux of radon on a daily basis (note that for semi-diurnal tidal periodicity we assume 1 day ∼ 2 tidal periods), is then 3.2 Stable isotope hydrology

Sampling location and timing
Water samples for stable isotope analysis were collected in triplicate from all possible water sources to the lagoon (endmembers) during winter on various occasions between 2007 and 2011 (Tables 2 and S1).These include the marine endmember, sampled in 2009; groundwater from local aquifer units (M10, M12, unconfined aquifer lenses in the barrier  For the latter, quasi-synoptic distributions of δ 18 O and δ 2 H in water at different tidal stages were obtained.For this purpose, we followed the division of the lagoon into two sectors, comprising western and eastern areas (see Fig. 1), with the separation line lying between the city of Faro and the Barra Nova.This division was based on the known divergent flow of groundwater in the M12 and M10 aquifers from a central point (Rio Seco-Chelote line, Fig. 1) as described (see Sect. 2.2) in Almeida (2000).High-powered boats were deployed, one from the city of Faro, on 2 December 2009, and the other from the city of Olhão, on 5 December 2009 (Fig. 1).The boats followed the tide outflow (or inflow) while covering all the pre-defined sampling points (western sector stations 1-5 and 1B to 5B, eastern stations A to I, Fig. 1).Each region of the lagoon was covered at each tidal stage in Table 2. Precipitation records over the region during the sampling campaigns described by this study, as measured at the São Brás de Alportel meteorological station (www.snirh.pt,Ref 31J/C).Monthly precipitation is contrasted with rainfall during the sampling campaigns and compared with historical monthly averages in order to evaluate the relative wetness of the periods in the wider temporal context.Accumulated precipitation during the 3 months prior to the month fieldwork took place is also shown and similarly compared to the historical record average.For a more detailed contextual assessment, the chronological record of daily precipitation for the period 2006-2013 is shown in Fig. 4d, with the sampling periods overlain for easy reference when evaluating the stable isotope hydrology of the catchment defined by this study and previous research.Under "Sampling", and "Type", the type of end-member collected for stable isotope analysis is shown, except when radon survey campaigns were executed in parallel -in this case "Radon survey" is added to the column.More details on the individual samples are shown in Table S1 no more than 2 h around slack tide.Coastal seawater adjacent to the Ria Formosa was sampled 2 nautical miles (∼ 3.8 km) offshore from the town of Quarteira to the west and from the Barra Velha (Armona inlet, Fig. 1, reference J).

Sampling and analytic methodology
Water was directly filtered through Rhizon SMS ™ membranes into sterile glass Vaccutainer ™ vials in the field.Subsequently, the cap area including the rubber septum was sealed with a layer of hot glue encased in Parafilm ™ .The vials were kept preserved at 4 • C until analysis could occur (typically within 6 months from the date of collection).
Samples were sent for standard analysis of δ 18 O and δ 2 H to GEOTOP Canada (Micromass Isoprime ™ dual inlet coupled to an Aquaprep ™ system), Durham University (LGRliquid water isotope analyzer, DT100) and at UFZ's stable isotope laboratory facilities in Halle, Germany (Laser cavity ring-down spectroscopy (Laser CRDS) Picarro water isotope analyzer L-1120i).Following standard reporting procedures (Craig, 1961a), delta values (δ) are reported as deviations in permil (‰) from the Vienna Standard Mean Ocean Water (V-SMOW), such that δ sample = 1000((R sample /R V-SMOW )-1), where R is the relevant isotopic ratio (i.e., either 2 H / 1 H or 18 O / 16 O).The mean analytical uncertainty is reported for each data point as ± 1 standard deviation (SD) of the mean of n analysis results obtained for n replicate samples in ‰ for δ 18 O and for δ 2 H (see Table 2).Each laboratory uses stringent protocols and reporting of stable isotope values using internationally calibrated standards; hence, reported stable isotopes values of water between the different labs used in this study are directly comparable.

Inter-annual comparability of isotopic data
Sampling campaigns where carried out strategically following a field-adaptive protocol.Of primary concern was to capture the extent of temporal end-member variability in isotopic signatures under maximum freshwater flow (high-flow) conditions, in order to (a) guarantee coherence of source compositions to feed into mixing models when necessary while assessing the hydrology of the lagoon over wider temporal scales and to (b) minimize logistics and costs while guaranteeing inter-comparability.For this purpose, the winter season was chosen given that ∼ 61 % of the mean annual precipitation falls on the region between November and February (34 % in the months of December and January).
Stable isotope sampling in winter had the added advantage of minimizing kinetic effects over stable isotope signatures given the lower evaporation potential.Sampling in winter 2007 was exploratory, with two main objectives: firstly, to characterize the isotopic signature of M12 groundwater and surface lagoon waters in the western sector, particularly in the area that could be potentially influenced by both SGD and the WWTP outflow under maximum dilution potential (hence high tide), and secondly, to conduct an exploratory survey of potential seepage areas along the Ancão peninsula, keeping in mind that the location of at least one of the important SGD seepage sites was known (Leote et al., 2008).Detection of the isotopic signature of groundwater in porewaters at the seepage face at stations Pw_e and Pw_f (Table S1) led to the installation at their location of a nested piezometer transect array in January 2010.This was subsequently used to obtain porewater samples in the 2010/2011 winter season (December 2010 and January 2011).
To capture inter-annual variability, the M12 aquifer was sampled twice (winters of 2007 and 2009), with the provision of one common location (Ramalhete) for cross-referencing.Following the same reasoning, the M10 aquifer was sampled in December 2010 while simultaneously sampling Rio Seco (belonging to M12, Table S1).This ensured intercomparability between groundwater isotopic signatures in 2009 and 2010.Campaigns were planned in advance considering the precipitation over the region to ensure similarity in the hydrological regime and ultimately guaranteeing intercomparability of results.The sampling itself took place in dry conditions as much as possible, and never after intensive rain that could have promoted flooding (Table 2, Fig. 4d).For example, while January 2007 was a dry month (8.8 mm) compared to the historical average (138 mm), the accumulated precipitation during the previous 3 months was 369.7 mm, consistent with the historical average (Table 2).By contrast, both December 2009 and 2010 were relatively wet months (392.2 and 269.6 mm), but followed relatively dry 3-month periods (Table 2).So porewater samples were also taken in January 2011, hence complementing winter 2010/2011.January 2011 followed a wet 3-month period (414.7 mm) and was hence comparable with January 2007, also relatively dry but on the back of 3 wet months (369.7 mm cumulative).The combined data set therefore contains results from repeated measurements for end-member isotopic composition under high-flow conditions across different years.These are in addition compared to historical data (Table S1, Fig. 4), leading to a temporally coherent quantitative overview of stable isotopic hydrology over the catchment.2 for summarized data).EMMWL: Eastern Mediterranean Meteoric Water Line (Gat and Carmi, 1970); WMMWL: Western Mediterranean Meteoric Water Line (Celle-Jeanton et al., 2001); GMWL: Global Meteoric Water Line (Clark and Fritz, 1997); LMWL: Local Meteoric Water Line (Carreira et al., 2005).

Spatial and temporal distribution
The activity ranges and spatial distribution of 222 Rn were similar in winter and spring.Because the weather was stormy during winter sampling, the uncertainties associated with determination of the radon evasion fluxes affecting the overall lagoon radon inventory were much higher than in spring (see Table 1).Indeed, using a mass balance, used estimate fluxes have been shown to be sensitive to parameterization of gas exchange (k) with the atmosphere, with potential uncertainties reaching 58 % (Gilfedder et al., 2015).Hence only the spring survey data are presented and discussed.Excess ).The highest local water column inventories (318 and 267 Bq m −2 during flood and ebb, respectively) were found in the Faro channel, covering stations 3 to A during ebb and 4 and 5 during flood.The eastern sector water column inventories were much higher during flood than during ebb.Given the non-random spatial distribution of radon, the median of each data set was used to calculate whole-lagoon inventories.The MAD (median absolute deviation, Hampel, 1974) was then used to propagate uncertainty in the radon budget calculations (Table 1).Radon inventories (median ± MAD) were 54.2 ± 17.8 and 74.0 ± 17.6 Bq m −2 , respectively, during ebb and flood (Table 1).

Along-channel tidal radon fluxes
Radon activity at Quatro Águas and Barra Nova was strongly anti-correlated with water level.At Quatro Águas, radon activities varied between 0 and 40 Bq m −3 , while at Barra Nova they varied between 1 and 31 Bq m −3 .Tidal variability at these two points was therefore consistent with the ranges in radon activities found during the lagoon survey.Time series of instantaneous Rn fluxes obtained as described by Eq. ( 1) are depicted for both locations in Fig. 3.The plots show consistency in the magnitude of upstream and downstream radon fluxes (grey area under the curves) through successive tidal cycles.The net daily tidal exchanges of radon through the Barra Nova and the Quatro Águas site (8.0 ± 0.5 × 10 4 and 9.9 ± 2.0 × 10 3 Bq d −1 , respectively) were both directed landward.This finding is consistent with the Barra Nova being a flood-dominated inlet (channeling ∼ 64 % of the flood and ∼ 59 % of the ebb prism of the Ria Formosa during spring tides: Dias and Sousa, 2009;Pacheco et al., 2010b).
To calculate the total residual exchange of radon between the Ria Formosa and the adjacent coastal area, we assumed the radon flux occurring at the other inlets to be proportional in equal measure to the individual residual tidal prisms.After adjustment to the lagoon surface area at MTL, the net exchange was just −9.3 (±1.6) × 10 −4 Bq m −2 d −1 (Table 1), so small as to be well within the uncertainty of all other quantities in the mass balance, implying that the radon inventory within the lagoon is controlled by internal fluxes.

SGD estimates based on radon mass balance
Solving Eq. (3) for a radon inventory of 65.9 ± 19.6 Bq m −2 (Table 1) gave a result for Rn adv of 7.14 ± 5.18 Bq m −2 day −1 , which adjusted to the submerged area at mean tide level (Tett et al., 2003) gives an SGD-derived radon flux of 4.14 (±3.00) × 10 8 Bq day −1 for the entire lagoon.Alternatively, the advective radon fluxes calculated as per Eqs.(4a) and (4b) for low and high tide periods were, respectively, 46.8 ± 38.8 and −32.5 ± 27 Bq m −2 day −1 .The positive and negative signs imply an advective flux of radon (Rn adv ) into the lagoon water column at low tide, while a net loss occurs during high tide.The resultant net Rn adv (Eq.4c) occurring during a full tidal period is 7.15 ± 8.4 Bq m −2 day −1 , statistically equivalent to the flux calculated via the assumption of steady state of the system over the lifetime of radon on a daily timescale (Eq.3), and yielding an equivalent SGD-derived radon flux of 4.14 (±4.87) × 10 8 Bq day −1 for the entire lagoon.
4  Jeanton et al., 2001).In conjunction with the average isotopic composition of groundwater in the catchment, that of seawater (Carreira, 1991) and adjacent coastal water, a precipitation-seawater mixing line (PP-SW Mix, Fig. 4) may be defined (δ 2 H = 5.37 × δ 18 O-1.7,r 2 = 0.99).The slope of this mixing line is similar to that found by Munksgaard et al. (2012) for the Great Barrier Reef (i.e., 5.66).Additional relationships framing the isotopic composition of the waters in the catchment in δ-space include the Local Meteoric Water Line (LMWL), defined by Carreira et al. (2005) as δ 2 H = (6.44 ± 0.24) × δ 18 O + (3.41 ± 1.13), and the Eastern Mediterranean Meteoric Water Line (EMMWL, Gat and Carmi, 1970).This is introduced as an extreme boundary to the isotopic composition of precipitation in southern Portugal.Indeed, rain with high d-excess originating either from the eastern Mediterranean or aligned with extreme precipitation events might fall in the region (see Fig. 4c), particularly during summer and/or fall (e.g., Frot et al., 2007) In winter 2007, the stable isotope composition of groundwater in M12 reveals slight evaporative enrichment by comparison to the GMWL and LMWL, plotting along the precipitation seawater mixing line (Fig. 4b).The isotopic compositions of surface waters (WWTP settling lagoons and lagoon surface waters) and porewaters plotted between the LMWL and the PP-SW mixing line (Fig. 4b), suggesting their composition was controlled by the interplay between the mixture of seawater and groundwater and evaporation-condensation cycles occurring along the hydrological travel path.In winter 2009, however, the range of isotopic compositions of surface water samples (∼ 2.87 ‰ for δ 18 O and ∼ 3.96 ‰ for δ 2 H) was significantly different (see inset, Fig. 4c).Their composition then fell between the WMMWL and the PP-SW mixing line.Even though the number of samples taken in winter 2007 was lower than those taken later and tide-specific sampling was absent, comparison of samples taken in both winters at high tide slack (Table S1; stations 2, 3, 4, A and 3B) shows the isotopic composition of water in the Ramalhete and Faro channels was distinct -the observed difference in range cannot therefore be attributed to the sampling strategy.Groundwaters across the catchment could be divided into three distinct groups: samples from Pechão Gimno, Pechão Serra and Pechão Zona industrial (Table 2), all from unit M10, plot above the GMWL and the LMWL, while samples taken from the unconfined aquifer wells in the outer barrier islands belonging to the unconfined M12 aquifer (i.e., Deserta, Table S1) plot distinctly below the PP-SW mixing line.In between, M12 samples plot along (Ramalhete) and below (Costa, Chelote, Rio Seco) the PP-SW mixing line.Samples from unit M10 plot along a local evaporation line (LEL) with slope ∼ 4.5, while samples from unit M12, excluding the ones located within the Ria Formosa, plot along a LEL with slope ∼ 4.1.

Isotopic composition of beach porewater
The porewater isotope compositions differed significantly between the winter of 2007 and that of 2010/2011.Beach groundwater was sampled both during spring and neap tides from sediment depths ranging from 50 cm to 3.5 m below MTL across a beach profile from the upper to lower intertidal during the latter period.δ 18 O ranged from 0.96 to −0.20 ‰ and δ 2 H from 2.5 to 8.5 ‰ and plotted close to the LMWL (Fig. 5a) along an evaporation line defined by δ 2 H = (4.02± 0.56) × δ 18 O + (4.51 ± 0.31), n = 24, r 2 = 0.702 (not shown).The slope of this LEL is slightly lower than those of the groundwater LELs (4.1 for the M12 and 4.5 for the M10).The data fell into three distinct groups (Fig. 5a, b) according to the relative position of the sampling point within the beach section.The first group of samples (average δ 18 O of 0.0 ± 0.13 ‰ and δ 2 H of 3.5 ± 0.93 ‰, n = 5) corresponded to the un-saturated and intermediate zones (upper intertidal), while the second (average δ 18 O of 0.4 ± 0.31 ‰ and δ 2 H of 6.1 ± 0.47 ‰, n = 10) and third groups (average δ 18 O of 0.7 ± 0.18 ‰ and δ 2 H of 8.0 ± 0.37 ‰, n = 9) were isotopically heavier and included in that order porewater from the deeper (> 2 m below the surface) and shallower (< 1 m below the surface) areas of the beach section.The respective average porewater stable isotope compositions plotted close to the LMWL (Fig. 5a), showing enrichment in opposition to distance from the surface in the saturated zone and depletion in the unsaturated recharge zone, probably due to capillarity effects (Barnes and Allison, 1988).The dependence of d-excess (Dansgaard, 1964) on δ 18 O (Fig. 5b) illustrates the deviation of porewater composition from Craig's (1961b) GMWL (δ 2 H = 8 × δ 18 O + 10) along significantly linear slopes dependent on local evaporation conditions.Indeed, porewater d-excess from deeper within the beach plots along the line defined by d = −6.7 ( ± 0.27) × δ 18 O + 5.57 (±0.13) (n = 10, r 2 = 0.987, P < 0.0001), while that from shallower areas plots along the line defined by d = −7.1 (±0.69) × δ 18 O + 7.28 (±0.52) (n = 9, r 2 = 0.937, P < 0.0001).These define slopes in δspace close to 1 and are consistent with the flow paths taken by beach groundwater between the seawater infiltration point at the higher beach face (higher d-excess) and the exfiltration point at the seepage face (lower d-excess).For the intermediate group of samples, longer flow paths (larger d-excess range) and less evaporative enrichment (lower average δ 18 O) are consistent with tidal-forced circulation at larger depths within the beach face.Conversely, shorter flow paths (relatively narrow d-excess range) and more evaporative enrichment (higher average δ 18 O) characterize shallower circulation pathways.
Inter-annual variability was also significant.The range of ∼ 1.16 ‰ for δ 18 O and ∼ 6.03 ‰ for δ 2 H found in 2009-2011 was 50 and 36 %, respectively, of the 2007 range, in spite of a common sampling location.Furthermore, isotopic compositions for porewater collected in 2007 plotted in δspace clearly in between the LMWL and the PP-SW mixing line (Fig. 4b), while the 2010-2011 samples overlap the LMWL (Figs. 4c and 5a).This occurs in spite of fewer samples being taken in 2007 and their depth of 50 cm below the surface, in contrast to the wide range of sediment depths sampled during 2009-2011.Paired ranges of porewater salinity also differ, varying between 21 and 36 in 2007 and between 36 and 43 in 2010 and 2011.These results suggest different water source functions were present during each sampling period.

Tidal variability of surface water δ 18 O and δ 2 H
Tides have a significant effect on the range of isotopic composition of surface water within the lagoon (see Fig. 6).In both lagoon sectors, the isotopic compositional range of water was much wider at low tide (Fig. 6a) than at high tide  S1) at different depth levels below the surface at the saturated zone and the dynamics of the beach groundwater table.(a) frames the compositional range and the subdivision of the isotopic characteristics through three groups, corresponding to different circulation paths within the beach (for an explanation, see Sect.4.2.3).(b) frames the same samples in a deuterium-excess (d) vs. δ 18 O plot, illustrating the progression of evaporative enrichment throughout the three zones and its relationship with the LMWL (Local Meteoric Water Line, Carreira et al 2005).Crosses and attached error bars represent average compositions for each group.Error bars represent the ± 1 SD PP-SW mix: precipitation-seawater mixing line (Sect.4.2.1);EMMWL: Eastern Mediterranean Meteoric Water Line (Gat and Carmi, 1970); WMMWL: Western Mediterranean Meteoric Water Line (Celle-Jeanton et al., 2001); GMWL: Global Meteoric Water Line (Clark and Fritz, 1997).(Fig. 6b), but this variability was more apparent in the western sector.During low tide there δ 2 H ranged from 5.3 ‰ (Station 2B) to 7.9 ‰ (Station 2) and δ 18 O from −0.82 ‰ (Station 2B) to 2.05 ‰ (Station 3).By contrast, 2 H ranged from 5.1 ‰ at Station 3B to 7.3 ‰ at 4B, while δ 18 O varied from −0.16 ‰ (Station 4) to 0.86 ‰ (stations 1B and 2B).The water mass at Station 2B was the most depleted in 18 O during low tide (Fig. 6a) and the most enriched in 18 O during high tide (Fig. 6b), but remains at the lower end of the δ 2 H range covered by all collected samples dur-ing both tidal stages.Aspects of tide-induced circulation are also revealed when the western and eastern sectors are compared for identical tidal stages (Fig. 6a, b).During low tide (Fig. 6a), the isotope compositions of water collected at the Ramalhete channel and the associated Ancão basin (stations 1B to 5B, Fig. 1) plot to the left of the LMWL, with the most isotopically depleted water found at Station 2B and the most enriched found at Station 1B.Conversely, water samples collected in the Faro channel (stations 1 to 5) plot to the right of the LMWL.The situation is reversed during high tide (Fig. 6b), with isotopic compositions of water from stations 1B to 4B plotting to the right of the LMWL as a result of mixing with seawater and coastal water and all others plotting to the left (mixing with internal lagoon water, including porewater).

Radon source attribution
In order to derive an SGD rate for the Ria Formosa, we divide the end-member source activity by the advective radon flux (4.14 ± 3.00 × 10 8 Bq day −1 ) calculated from the mass balance.However, because radon budgets include 222 Rn sourced in seawater recirculation, the discharging fluid composition is important for discriminating between potential sources of SGD.In fact, the two modes of SGD may be separated according to whether they drive a net influx of freshwater to the system (Santos et al., 2012).Indeed, there are three identified potential sources for advective radon input to the lagoon, i.e. (Table 1), water in freshwater lenses under the outer barrier islands (outer reaches of the M12 aquifer) represented by the Deserta well (mean 0.95 salinity), porewater in sandy beaches (mean 40.6 salinity) mobilized by tidal pumping (seawater recirculation), and, finally, meteoric water travelling through the subterranean pathway (M12 aquifer), represented by samples taken from the Ramalhete borehole (mean 5.06 salinity).The corresponding volumetric discharges, if each of these potential sources is considered in turn to be the only source of SGD into the lagoon, are 4.42 (±4.25) × 10 6 , 1.36 (±1.28) × 10 6 and 6.26 (±4.63) × 10 4 m 3 day −1 , corresponding, respectively, to ∼ 3.16, ∼ 0.97 and ∼ 0.04 % of the mean daily flood prism (1.40 × 10 8 m 3 ).When defining the radon source function, salinity is occasionally used as the discriminating parameter because of its conservative nature (Crusius et al., 2005;Swarzenski et al., 2006;Stieglitz et al., 2010).However, the low estimated SGD to tidal prism ratio combined with saline intrusion into the local aquifers (Silva et al., 1986; Table S1) advises against this option, as the estimated discharge volumes would not have a discernable impact on the overall salinity of the Ria Formosa, leaving us without a way with which to verify the reliability of the choice.Furthermore, porewater salinity at the site where the piezometer transect is located (Fig. 1) was always very high (> 35; Table S1), but could be as low as 21 in 2007, suggesting different SGD modes might be active in different years.So how do we confidently identify the source of radon?
Our mass balances (see Sect. 4.1.3)for each tidal stage suggest that radon is removed from the water column during the flood period.In the absence of any other realistic explanation we might accept that it had to be advected into the unsaturated intertidal zone during beach recharge.The daily flux of radon into unsaturated sandy sediments would then amount to 16.25 ± 13.5 Bq m −2 day −1 .Conversely, the input of radon into the water column during ebb was 23.4 ± 19.4 Bq m −2 day −1 .Because the mean radon inventory during high tide was 19.3 ± 4.74 Bq m −3 , a flux of 16.25 ± 13.5 Bq m −2 day −1 into unsaturated sediments would equate to a beach recharge rate of ∼ 1.2 m day −1 .This figure is consistent with the discharge rates measured during 2006 by Leote et al. (2008) at the lower intertidal, which reached 1.9 m day −1 .If we therefore assume that beach discharge balances recharge on a volumetric basis at daily timescales, then the area of water infiltration would be ∼ 1.13 × 10 6 m 2 .Given the porosity of sandy beach sediments on site of ∼ 0.3-0.4(Rocha et al., 2009), recharge would need only occur through ∼ 7.5-10 % of the maximum surface intertidal area of the lagoon (see Sect. 2.1).Hence tidal pumping is a realistic explanation for the radon advected into the water column on a daily basis.Still, the radon data alone do not provide irrefutable proof that SGD estimated through the radon mass balance for June 2010 originates from seawater recirculation through beaches and porewater exchange mechanisms.
This proof is important: an example of how an unsupported choice of radon end-member might significantly affect quantification of nitrate loading to the lagoon through SGD could be given at this stage to illustrate the effects of the lack of irrefutable source attribution.The mean nitrate concentration (in mg L −1 , spring tides, 2009 to 2011) was 0.1 for the lagoon water column, 0.81 for beach porewaters, 2.22 in the Deserta well, and 130 for the Campina de Faro aquifer (M12).Our discharge estimates based on the Hydrol.Earth Syst.Sci., 20, 3077-3098, 2016 www.hydrol-earth-syst-sci.net/20/3077/2016/ radon balance would then result in potential average SGDborne nitrate loading to the Ria Formosa of 0.96, 9.8 and 8.14 Ton N day −1 , if the source of excess radon was, respectively, seawater recirculation through beach sands or fresh groundwater originating from either the lens under the dune cordon or the landward section of M12 aquifer.Two cautionary notes on these numbers should be obvious: (a) the latter would drive net N additions to the lagoon water budget, while the former would not, implying that (b) the loadings based on directly multiplying fresh SGD by the average nutrient concentrations found in the end-member samples ignore any transformations occurring within the interface before the mixture arrives at the lagoon, and therefore are likely to be overestimated.Ferreira et al. (2003) estimated total N fluxes to the lagoon at 1028 Ton N yr −1 (2.82 Ton N day −1 ), with 58 % (1.64 Ton N day −1 ) originating from diffuse sources.Simple extrapolation from our data would suggest that ∼ 34 % of the total N fluxes to the lagoon, and ∼ 59 % of the nonpoint source loading, would arise from seawater recirculation through beaches, while the meteoric SGD sources would multiply the total N loading into the system by a factor of 6 or 5 on a daily basis, depending on the composition of fresh groundwater.These two latest figures compound our cautionary notes above.Furthermore, during winter 2010 and 2011, when porewater salinities were very high, nitrate available in porewaters at the littoral fringe was likely sourced from benthic mineralization of local organic matter (autochthonous source) and not in fresh groundwater input.Conversely, because nitrate contamination of the Campina de Faro aquifer is anthropogenic, freshwater inflow via SGD into the lagoon would also define the associated nitrate inputs as allochthonous, or "new" contributions to the system's nutrient budget.Depending on SGD mode, therefore, there would be an order of magnitude difference between allochthonous and autochthonous sources of nitrate into the lagoon, even if the former might be overestimated as discussed.Accurately identifying the SGD source function would therefore be absolutely necessary to understand the biogeochemical workings of the lagoon, but this is not possible with the radon data alone, even in combination with the salinity data.
However, the stable isotope signatures of surface water bring clarity to the problem.The Local Evaporation Line (LEL-1, Fig. 6a) fitted by linear regression of the samples taken within the eastern sector at low tide intersects the LMWL close to the average isotopic signature of beach porewater in the unsaturated zone (Figs.5a and 6a).This indicates the original composition of the surface water before evaporation and mixing take place within the lagoon.The origin of the surface water is the recharge into the unsaturated beach area, which then reveals isotopic enrichment in proportion to its permanence within the system and the consequent extent of evaporative loss.Indeed, water in the upper intertidal at low tide will see its isotopic signature depleted within the sedimentary matrix -in the unsaturated zone, the isotopic concentration decreases quickly from a maximum at the zone of evaporation (phreatic surface) within the sediment matrix to a minimum close to the surface because of the movement of water vapor through the pores toward the surface (Barnes andAllison, 1983, 1988).While this is clear for the eastern sector, within the western sector there is another surface source of water (WWTP) that further complicates the picture.This water joins the lagoon close to Station 2B (Fig. 6a).So, the porewater in the unsaturated sediments mixes over time with the lagoon recharge at high tide and water already present within the tidal wedge (cf.Robinson et al., 2007), whereupon it leaves during beach discharge at low tide, either through shallow or deeper flow paths (Fig. 5b), and mixes with other meteoric sources and seawater (MX-1, MX-2, Fig. 6a).
For the period between the winter of 2009 and that of 2010/2011, therefore, the combined stable isotope and radon tracer approach allows definite attribution of the SGD source into the Ria Formosa.SGD arises from seawater recirculation through the permeable beach sediments of the lagoon driven by the tide.In the absence of meteoric SGD inputs, a significant amount of the tidal prism (∼ 1 %) circulates through local sandy sediments driven by tidal pumping, at a rate of ∼ 1.4 × 10 6 m 3 day −1 .This implies that the entire tidal-averaged volume of the lagoon (140 × 10 6 m 3 ) is filtered through its sandy beaches within 100 days, or about 3.5 times a year.Based on our nutrient data, the average nitrate loading driven by this SGD mode to the Ria Formosa can now be confidently put at an average of 0.96 Ton N day −1 , ∼ 59 % of the non-point source nitrogen loading estimated by Ferreira et al. (2003).
Salinity (see Table S1) does not correlate well with either δ 18 O and δ 2 H, though, particularly for samples with δ 18 O > 1 ‰ and/or δ 2 H > 1 ‰ and S > 37 ‰.With reference to surface water δ 18 O values, these comprise, respectively, the most isotopically enriched waters found during low tide and the innermost stations in the eastern sector (stations G,H and F;Figs. 1 and 6a) and at locations within the Faro channel (stations 1-4; Figs. 1 and 6a), as discussed earlier.It is also the case for most porewater samples.Indeed, even if the mean composition of porewater from different sections of the beach plots along well-defined mixing and evaporation lines (Fig. 5a, b), the average salinities of each group do not change significantly with δ 18 O  enrichment (40.2 ± 1.78, 40.6 ± 2.57 and 40.6 ± 2.07, respectively).While this observation is consistent with theory (Craig and Gordon, 1965) and previous analysis of the covariance of δ 18 O, δ 2 H and salinity in seawater (Rohling, 2007), it also implies that the joint use of these tracers to infer the relative contribution of different source functions has to be done with care in semi-confined coastal water bodies subject to significant evaporation.As further support to this observation, we note that the mixing lines (MX-1 and MX-2, Fig. 6a) between the porewater within the beach tidal wedge and the most enriched waters found in the www.hydrol-earth-syst-sci.net/20/3077/2016/ Hydrol.Earth Syst.Sci., 20, 3077-3098, 2016 western sector (δ 2 H = (0.97 ± 0.08) × δ 18 O + (5.70 ± 0.09), r 2 = 0.87, n = 21) and between the Ramalhete channel and Ancão basin (stations 3B, 4B, and 5B) and the water mass near Olhão at stations G and H (δ 2 H = (1.02± 0.18) × δ 18 O + (7.13 ± 1.01), r 2 =0.84, n = 16) are virtually the same as that characteristic of the modern surface ocean (δ 2 H = 1.05 × δ 18 O + 6.24, r 2 = 0.21, n = 62) within a comparable salinity range (Rohling 2007).This observation suggests, in coastal ocean regions and areas of restricted exchange like lagoons, that the stable isotope signature of seawater reflects important contributions arising from porewater exchange driven by tidal pumping, amongst other mechanisms.Identifying and discriminating these contributions brings insights also into the hydrological paths active within these systems and therefore provides an invaluable tool for supporting reliable biogeochemical budgets.

Hydrological pathways and dispersion of SGD in the Ria Formosa Lagoon
The amount-weighed isotopic composition of precipitation over Faro (GNIP: IAEA/WMO, 2013) plots (Fig. 4a) at the intercept point of the GMWL, the LMWL (slope ∼ 6.4) and the precipitation-seawater mixing line (slope ∼ 5.4).The isotopic signature of precipitation hence plots close to that of groundwater, indicating that local aquifers are directly recharged by precipitation, in agreement with prior reports (e.g., Engelen and van Beers, 1986).The isotopic composition of surface waters also reveals that the lagoon and the adjacent coastal water may be classified as a coastal boundary zone similar to that described elsewhere (Blanton et al., 1989(Blanton et al., , 1994;;Moore, 2000), in which the isotopic signatures result from the mixing between offshore seawater and continental meteoric sources affected by surface evaporation.Accordingly (Fig. 6), the stable isotope composition of water within the lagoon varies with tidal stage and will be affected on the one hand by the magnitude, origin and pathways taken by the meteoric inputs and on the other by internal mixing, driven by lagoon hydrodynamics and by the local evaporation regime.Nevertheless, the porewater endmember is part of the surface water mixture in both sampled periods, although in different ways: some porewaters (Pw_e and Pw_f; see Table 2) collected at the same site were significantly more depleted in both 18 O and 2 H during 2007 (Fig. 4b) when compared to 2009-2011 (Fig. 4c), and these are characterized by comparatively low salinities (21 and 23, Table 2).Station 2B is the closest to the Faro WWTP outlet; during low tide the water mass joining the lagoon mixture there has an isotopic signature close to the Western Mediterranean Water Line (Fig. 6a), suggesting that a meteoric source of water joins the lagoon there, presumably as part of the WWTP discharge.On the other hand, the exchange in position of the isotopic signature of water at stations 1-5 and 1B-3B with reference to the LMWL in δ 18 Oδ 2 H-space during flood (Fig. 6b) suggests a hydrodynamic connection between the Ramalhete channel, the Ancão inlet and the water masses in the eastern sector.This connection would occur via the Faro-Olhão inlet and associated channels as ebb progresses onto flood, linking both the stations closest to the city of Olhão (stations E, F, and G) and the ones closer to the coastal ocean (stations A, B, and C), to the water masses originally present in the western sector.Indeed, stations 1 to 4 in the Faro channel display depletion of 18 O during high tide (Fig. 6b) by comparison to low tide (Fig. 6a).This provides evidence that the meteoric source present within the Ramalhete channel also influences the water in the Faro channel during high tide.Furthermore, the isotopic data suggest that part of the water mass out flowing through the Ramalhete channel during ebb tide (stations 2B-5B) eventually ends up being present at stations F, G and H close to the city of Olhão via the inner portion of the system (Station 1B), having mixed with shallow beach groundwater (MX-2 in Fig. 6a), while water from the same region might also be led to stations A, B and C in the eastern sector via Station 5 after mixing through the beach water table (MX-1 in Fig. 6a, b).The dominant alongshore drift in the area is eastward, and, in fact, Pacheco et al. (2010) show that a strong hydraulic connection exists between the the Ancão, Barra Nova (Faro-Olhão) and Armona (Barra Velha) inlets, whereby the excess flood prism at Barra Nova is directed toward both the Ancão and Armona ebb-dominated inlets.The combination of data indicates that the body of water ebbing in the first instance through the Ramalhete channel is partially retained within the system and ends up in the Faro channel before the subsequent flood moves it eastward, either via an internal pathway eastward from the Ancão inlet basin and/or externally, looping back into the lagoon via the Faro-Olhão inlet after exiting through the Ancão inlet (Fig. 6a, b).
The combination of flood lag time between the Ancão and Barra Nova inlets, the eastward alongshore drift and the meteoric source of water at the WWTP plant outlet (closest to Station 2B) creates the characteristic inversion observed in δ 18 O-δ 2 H relationships and highlighted in Fig. 6a,  b.This circulation path inferred from the isotopic composition of water is also consistent with the radon data, since the radon-enriched water masses found in the Ramalhete and Faro channels (Fig. 2a) during low tide would eventually be transported toward the eastern sector via the distribution of the excess flood prism at Faro-Olhão (Pacheco et al., 2010).This would help explain why the radon inventory in the eastern sector is higher during flood tide (Fig. 2b) and why the net exchange of radon is directed into the lagoon at both Quatro Águas and Barra Nova (Table 1), as part of the radon associated with beach seepage would be retained in the lagoon and/or transported back into the system via the Barra Nova after exiting through the Ancão inlet.

Inter-annual comparison of lagoon hydrology using deuterium-excess
Because of the relatively higher enrichment in 18 O compared to 2 H in the residual water (Gat, 1996) = 0.992, n = 14, P ∼ 0).These relations reveal the two different pathways into the ria followed by groundwater from the M12 aquifer in 2007 (Fig. 7a).The surface water circulation pathway (P1) originates when water from the public supply (sourced in local aquifers) is treated at the WWTP and subsequently discharged into the lagoon, whereupon it circulates into the Ancão basin, mixing with coastal water and seawater.This pathway is consistent with the internal circulation path discussed earlier.In contrast, the groundwater pathway (P2) followed by water originating in the same aquifer crosses the subterranean estuary and emerges later (the d − δ 18 O correlation slope magnitude is higher than P1) within the lagoon where it mixes with surface waters, including seawater and the WWTP outlet emissions (Fig. 7a).Hence the isotope data conclusively show two aspects of the local water balance in 2007: on the one hand, water for public consumption was essentially extracted from groundwater sources, while on the other, SGD into the lagoon comprising a net water input into the system was present.
The situation later (2009-2011) was substantially different (Fig. 7b = 0.979, n = 37, and P ∼ 0. These highlight other aspects of the local water balance.Firstly, P3 suggests that groundwater from the M10 aquifer mixes with water in M12, and that the local groundwater flow follows a northeasterly to southeasterly general direction (cf. the locations of M10 and M12 in Fig. 1), eventually communicating under the Ria Formosa with freshwater lenses present in the barrier islands, where the d-excess signature of groundwater is lowest.Secondly, P4 shows that water used for public consumption in the catchment was mainly withdrawn from a direct meteoric source (position of rainwater signature, Fig. 7b).This water, upon leaving the WWTPs, then mixes with surface and re-circulated seawater, establishing the mixing line for the lagoon (Figs.6a and 7b).It is also evident that the surface water samples collected in the lagoon in 2007 plot close to the P4 line, suggesting that the magnitudes of the factors driving evaporation and internal circulation in the lagoon are generally stable on a multi-annual basis.This comparative approach confirms, additionally, that the subterranean pathway was not present in 2009-2011, and hence SGD at this time was comprised entirely of saline water re-circulated through the sandy beaches by tidal pumping.
The difference observed in water sources for public water supply and their isotopic signature in the catchment and subsequently released through the WWTPs into the lagoon is consistent with the changes occurring in the regional water management strategy: while water to meet irrigation and public consumption demand relied almost entirely on groundwater abstraction until the 2000s (Stigter et al., 2006), from this period onwards it was to be drawn almost exclusively from surface reservoirs north of the littoral zone.However, a substantial number of the local groundwater captions remained active in support of irrigation, while some of the major municipal captions had to be re-activated after the 2005 drought (EM-DAT, 2013) to support consumption demand when surface reservoirs became depleted.In fact, because of the unpredictability of scarcity periods, the current operational thinking tends toward mixing both water sources to face demand, with the primary source being surface water reservoirs (Monteiro and Costa Manuel, 2004;Stigter and Monteiro, 2008).Our approach clearly indicates that this is the case for 2009-2011 as the WWTP plant water signal shows the water being discharged as meteoric in origin (Figs.6a and 7b).Following the implementation of a mixed source water supply chain, the activity of the SGD subterranean pathway into the ria becomes dependent on whether groundwater levels in M12 are sufficient to establish a hydraulic gradient driving the flow, as was apparently the case in 2007 (Fig. 7a).Increased water mining and reduced aquifer recharge would provide the counterbalance by reducing groundwater levels and consequently the hydraulic gradient driving SGD of meteoric origin into the system via the subterranean estuary.

Concluding remarks
We compared hydrological scenarios in a semi-arid coastal lagoon across two different periods, aiming to distinguish SGD modes and correctly identify end-member contributions to the water mixture within the system.While it has been established that radon mass conservation allows for the determination of total SGD, i.e., meteoric plus re-circulated water flow, we show that combining this information with stable isotope hydrology contributes to define and distinguish origins and pathways followed by SGD into the system.While δ 18 O and d-excess paired data helped define the active hydro-logical pathways in the Ria Formosa, δ 2 H vs. δ 18 O plots provided insights into water source functions and their dispersion through the lagoon.Using our combined approach, SGD occurring in the Ria Formosa could be separated into a discharge incorporating net meteoric water input into a receiving ecosystem (2007) and an input with no net water transfer (2009)(2010)(2011).We conclude that whilst the Ria Formosa receives SGD through tidal pumping (as in [2009][2010][2011], it is also occasionally subject to SGD inputs of meteoric origin (as in 2007) directly associated with the contaminated M12 aquifer.
In the absence of meteoric SGD inputs, part of the tidal prism circulates through local sandy sediments driven by tidal pumping, at a rate of ∼ 1.4 × 10 6 m 3 day −1 .This implies that the entire tidal-averaged volume of the lagoon (140 × 10 6 m 3 ) is filtered through its sandy beaches within 100 days, or about 3.5 times a year, driving an estimated load of ∼ 350 Ton N yr −1 into the lagoon.Conversely, using the estimates for the upper bound of N concentration found in the freshwater component of SGD during 2006 (0.4 mmol L −1 ) and the associated SGD-borne freshwater discharge of ∼ 1.1 × 10 7 m 3 yr −1 estimated by Leote et al. (2008) based on seepage meter measurements, meteoric SGD inputs could add a further ∼ 61 Ton N yr −1 to the lagoon.If for the former the source is autochthonous and responsible for a rather large fraction (59 %) of the estimated nitrogen inputs into the system via non-point sources (Ferreira et al., 2003), leaving no direct mitigation options in the context of environmental management, it is not so for the latter, as specific measures could be implemented in support of mitigation (e.g., Almasri and Kaluarachchi, 2004).Nevertheless, the potential loadings delivered from two distinct vectors differ in magnitude, frequency and origin, and could therefore cause different ecosystem-level impacts.Hence while simple or weighted averages of end-member radon activities might be useful under well-defined circumstances (Crusius et al., 2005;Swarzenski et al., 2006;Kroeger et al., 2007;Blanco et al., 2011) in radon budgets to evaluate SGD as a potential pollutant source in comparison to other vectors (local surface drainage, riverine input, etc.), these are of little value for effectively providing environmental managers with the causal chain alluded to in the introduction: without actual source identification and attribution, there is little that can be done to manage potential pollutant loading of coastal ecosystems via SGD.
The Supplement related to this article is available online at doi:10.5194/hess-20-3077-2016-supplement.

Figure 1 .
Figure1.Map showing the locations of the sampling sites within the Ria Formosa and its geographical context.The top panel shows the full geographical extent of the system, with the operational separation of the region of interest into western and eastern lagoons and the names of all the inlets.The lower panel shows an amplified map of the region of interest, including major channels, locations of sampling and tidal stations, as well as boundaries of the aquifers bordering the lagoon (M10, M11, M12).

Figure 2 .
Figure 2. Map showing the distribution of radon inventories (Bq m −2 ) within the main channels, during ebb (a) and flood (b), for the radon survey conducted in 2010.For more details regarding the radon budget in both December 2009 and June 2010, see Table1.
island) taken from boreholes and wells (Fig.1), in January 2007 and December 2009 and 2010; precipitation, taken at the city of Faro in December 2009; and beach porewater collected in January 2007, December 2010 and January 2011.In 2007, samples were extracted from 50 cm below the sediment-water interface at various locations along the Ancão peninsula's inner dune cordon (Fig.1), while in 2010 and 2011 they originated from various depths in the sediment (2 to 7 m below r.m.s.l.) and where collected using a cross-shore array of nested, multi-level sampling piezometers (Fig.1) installed in the inner margin of the outer dune cordon in January 2010 at the point of maximal freshwater seepage rates found in 2007.Surface water reservoirs near Quinta do Lago used for irrigation and settling lagoons in the wastewater treatment plant near the city of Faro (WWTP) were sampled in July 2007, the River Gilão (Fig.1) in December 2010, and surface water from the lagoon was sam-pled during flood tide (western sector, Fig.1) in January 2007 and during both high and low tide in December 2009.

Figure 3 .
Figure 3. Tidal variability of instantaneous radon fluxes, respectively, at the inner Barra Nova inlet (a) and Quatro Águas station (b), for the radon survey conducted in 2010.For more details on calculation methods, please see Sect.3.1.2.

Figure 4 .
Figure 4. Catchment isotope hydrology.Counterclockwise, from top left: (a) shows the main meteoric water lines framing the isotopic composition of precipitation within the catchment, including the precipitation-seawater mixing line (PP-SW Mix, Sect.4.2.1).(b) plots the isotopic compositional range of water samples taken during 2007, while (c) plots the isotopic compositional range of water samples taken during the period 2009-2011; the lagoon surface water samples (inset) are shown in more detail in Fig. 6.(d) provides the complete record of daily precipitation over the region for the period 2006-2013 for contextual support (see also Table2for summarized data).EMMWL: Eastern Mediterranean Meteoric Water Line(Gat and Carmi, 1970); WMMWL: Western Mediterranean Meteoric Water Line(Celle-Jeanton et al., 2001); GMWL: Global Meteoric Water Line(Clark and Fritz, 1997); LMWL: Local Meteoric Water Line(Carreira et al., 2005).
radon activities measured in water varied between 3.5 and 37 Bq m −3 , with a narrower range (5-25 Bq m −3 ) measured during ebb.The highest activities within the western sector during this stage (> 25 Bq m −3 ) were measured close to the city of Faro and in the Ramalhete channel, and close to the city of Olhão (∼ 20 Bq m −3 ) in the eastern sector.Radon activities generally declined from the northwest to the southeast during ebb tide, with the lowest values (∼ 5 Bq m −3 ) found in the Olhão channel northeast of the Barra Nova.Conversely, the lowest activities during flood (∼ 5 Bq m −3 ) were measured close to the Ancão inlet and at the outer end of the Faro channel, suggesting radon-poor coastal water intrusion during flood tide.The mean radon activities throughout the lagoon were 19.3 ± 4.74 and 15.59 ± 4.54 Bq m −3 , respectively, during flood and ebb.Relative accumulation of radon occurred at specific locations in the lagoon (Fig.2a,Hydrol.Earth Syst.Sci., 20, 3077-3098, 2016   www.hydrol-earth-syst-sci.net/20/3077/2016/ b

Figure 5 .
Figure 5. Isotopic composition of porewater extracted in winter 2010/2011 (TableS1) at different depth levels below the surface at the saturated zone and the dynamics of the beach groundwater table.(a) frames the compositional range and the subdivision of the isotopic characteristics through three groups, corresponding to different circulation paths within the beach (for an explanation, see Sect.4.2.3).(b) frames the same samples in a deuterium-excess (d) vs. δ 18 O plot, illustrating the progression of evaporative enrichment throughout the three zones and its relationship with the LMWL (Local Meteoric Water Line,Carreira et al 2005).Crosses and attached error bars represent average compositions for each group.Error bars represent the ± 1 SD PP-SW mix: precipitation-seawater mixing line (Sect.4.2.1);EMMWL: Eastern Mediterranean Meteoric Water Line(Gat and Carmi, 1970); WMMWL: Western Mediterranean Meteoric Water Line(Celle-Jeanton et al., 2001); GMWL: Global Meteoric Water Line(Clark and Fritz, 1997).

Figure 6 .
Figure 6.Tidal variability of the isotopic composition of surface waters in the lagoon, framed by significant local evaporation (LEL), mixing (MX), and meteoric lines, as well as the average composition of adjacent coastal water and seawater (historic data).(a) Low tide and (b) high tide.For more details, see Sects.4.2.4 and 5.2.

Figure 7 .
Figure 7. Hydrological pathways within the Ria Formosa, as defined by stable isotope data.(a) 2007 situation -SGD with net input of meteoric water present; (b) 2009-2011 -SGD essentially derived from tidal pumping.Detailed explanations are available in Sect.5.3.
, deuterium-excess (dexcess = d = δ 2 H -8 × δ 18 O) decreases in water as evaporation progresses (i.e., as δ 18 O increases).It follows therefore that a plot of d-excess vs. δ 18 O (in a similar fashion to Fig. 5b for porewater) might reveal the path taken by a particular water mass within a catchment area, because (a) the magnitude of the fractionation imposed by evaporation along the travel path affects the d-excess of residual water (setting the slope of paired d − δ 18 O relationships), and (b) water of different origins would have different d-excess values.The slope of the d − δ 18 O covariance line shows the deviation of isotopic compositions from Craig's meteoric water line(Craig, 1961b).Therefore its magnitude in absolute terms is proportional to the extent of evaporative enrichment, a function of the exposure time of the water to evaporation.Conversely, following the line along decreasing δ 18 O values would lead us to the original isotopic composition of the water, set before the evaporative regime changed.These characteristics allow us to disentangle and identify the main hydraulic pathways active in the Ria Formosa and compare the two periods under scrutiny to reveal the distinct nature of SGD within the system (Figs.5b and 7a, b).Accordingly, four significant d − δ 18 O correlation lines are identified in the basin (Fig.7).In 2007, two pathways (P1 and P2) connecting the composition of M12 groundwater with water sampled in the lagoon are revealed: P1, with d = (−1.10± 0.02) × δ 18 O + (4.41 ± 0.1), r 2 = 0.997, n = 6, P ∼ 0; and P2, with d = (−1.85± 0.05) × δ 18 O + (0.72 ± 0.11), r 2 ).Two major hydraulic pathways are shown in the isotopic data (P3, P4): P3, with d = (−7.8± 1.2) × δ 18 O-(22.76± 5.04), r 2 = 0.813, n = 10, and P = 0.0002; and P4, with d = (−7.43± 0.18) × δ 18 O + (6.45 ± 0.18), r 2