The effects of subterranean estuary dynamics on nutrient resource ratio availability to microphytobenthos in a coastal lagoon

be regulated by the annual oscillation of the local STE, itself driven by groundwater recharge variability.WestudySTEout ﬂ owsamplesgatheredmonthlyforayearin theRiaFormosalagoon,examiningthetem- poraldynamicsofsalinity,E H ,pH(Totalscale),dissolvedoxygenandnutrient(PO 43 − ,NO 2 − ,NO 3 − ,NH 4+ ,andSi(OH) 4 ) concentrations under the local hydrological regime. The objectives were threefold: (1) to determine the annual vari- ability ofnutrient content andN:P:Sistoichiometryin SGD intothelagoon;(2) toidentifythe main driversofvariabil-ity in SGD composition and stoichiometry and their interactive effects; (3) to discuss links to, and implications for, ecosystem function that could help de ﬁ ne expectations of cause-effect relationships and be useful for environmental management of the lagoon and similar systems elsewhere. We ﬁ nd that the terrestrial groundwater recharge cycle drives the expansion and contraction of the subterranean es-tuaryon annualtimescales,causingthepHofSGDto ﬂ uctuateinoppositiontocontinentalgroundwaterlevel.Thean- nual dynamics of the STE and the resulting pH oscillation determine the annual variability of nutrient composition ratio in SGD and shape benthic primary production dynamics. When saltwater intrusion occurs, the pH within SGD increases, enhancing nitri ﬁ cation and desorption of exchangeable phosphorus, while silicate ﬂ uxes increase with sea- waterretreat.Theresultisthatnutrientresourceratioavailabilityforcoastalprimaryproductiondependsonthefresh groundwater level. This implies that ecosystem function in such systems is more tightly related to the dynamics of linked groundwater reservoirs than previously thought.


H I G H L I G H T S G R A P H I C A L A B S T R A C T
• The pH of subterranean estuary (STE) outflows covaries with annual variability in continental recharge. • Landward movement of the saltwater wedge toe drives pH higher in the STE, promoting nitrification and phosphate release. • Groundwater recharge dynamics drive annual variability in nutrient composition ratio of submarine groundwater discharge (SGD). • Coastal lagoon primary production is interconnected with the oscillation of the saltwater/freshwater interface in STE. • The continental groundwater recharge cycle regulates coastal lagoon benthic primary production dynamics.

A B S T R A C T A R T I C L E I N F O Editor: Jurgen Mahlknecht
Keywords: Groundwater recharge dynamics Subterranean estuary Submarine groundwater discharge Semi-arid climate Coastal lagoon Benthic primary production Nutrients Causal links between subterranean estuary (STE) dynamics, their climatological drivers, and the ecology of coastal ecosystems have remained elusive. Yet, establishing these connections is essential for fully integrated management of coastal ecosystems. We test, in a semi-arid climate, whether the composition of submarine groundwater discharge (SGD) to a lagoon can be regulated by the annual oscillation of the local STE, itself driven by groundwater recharge variability. We study STE outflow samples gathered monthly for a year in the Ria Formosa lagoon, examining the temporal dynamics of salinity, E H , pH (Total scale), dissolved oxygen and nutrient (PO 4 3− , NO 2 − , NO 3 − , NH 4 + , and Si(OH) 4 ) concentrations under the local hydrological regime. The objectives were threefold: (1) to determine the annual variability of nutrient content and N:P:Si stoichiometry in SGD into the lagoon; (2) to identify the main drivers of variability in SGD composition and stoichiometry and their interactive effects; (3) to discuss links to, and implications for, ecosystem function that could help define expectations of cause-effect relationships and be useful for environmental management of the lagoon and similar systems elsewhere. We find that the terrestrial groundwater recharge cycle drives the expansion and contraction of the subterranean estuary on annual timescales, causing the pH of SGD to fluctuate in opposition to continental groundwater level. The annual dynamics of the STE and the resulting pH oscillation determine the annual variability of nutrient composition ratio in SGD and shape benthic primary production dynamics. When saltwater intrusion occurs, the pH within SGD Science of the Total Environment 851 (2022) 157522

Introduction
Submarine Groundwater Discharge (SGD) is 'all flow of water on continental margins from the seabed to the coastal ocean, regardless of the fluid composition or driving force' (Burnett et al., 2003). It usually consists of a mixture of fresh groundwater and seawater recycled through permeable sediments (Moore, 1999(Moore, , 2010, and often within subterranean estuaries (STEs) at the coast . In contrast to channeled surface inputs (e.g., rivers), outside of karst-dominated seaboards SGD occurs mainly as diffusive flow wherever a coastal permeable aquifer with a positive head relative to sea level connects to the ocean (Johannes, 1980). This makes SGD difficult to quantify as a nutrient transport vector, and early research largely ignored its contribution to marine biogeochemical budgets (Burnett et al., 2003). More recently, however, radium tracer studies have shown that SGD is the dominant water flux within coastal regions on a global scale (Kwon et al., 2014) and has therefore a considerable impact on coastal nutrient biogeochemistry worldwide Zhou et al., 2019). This is because SGD acts as a large nutrient source at local to regional scale, particularly when the freshwater portion originates in contaminated aquifers (e.g., Cho et al., 2018;Rocha et al., 2015;Leote et al., 2008;Ullman et al., 2003).
Nutrient fluxes associated with SGD have been evaluated to date at >50 study sites spanning 4 continents ). Yet, our understanding of the nutrient stoichiometry of SGD, its temporal variability and its controls is absent, mostly because STEs are still poorly understood . More than the magnitude of SGD-borne nutrient inputs to the coastal ocean, it is important to understand whether these are sourced internally within the land-ocean continuum (mineralization of organic matter within the coastal aquifer or coastal system), or externally to it (e.g., contamination of aquifers by nitrogen fertilizers). Indeed, autochthonous (internal) nutrient loads barely influence multiannual trends of whole-system net primary production, but allochthonous (external) nutrient inputs cause net nutrient enrichment over time and drive 'new' primary production, causing eutrophication and environmental degradation (e.g., Howarth and Marino, 2006). Nutrient loading may also shift species composition (Sugimoto et al., 2017): if net ecosystem production, or its potential rate, is limited by access to nutrients, the temporal dynamics of nutrient ratio resource availability in coastal waters may be followed by an increase in biomass of opportunistic species that take advantage of shifts in relative nutrient concentrations (Lecher and Mackey, 2018;Rocha et al., 2002). These shifts might trigger harmful algae blooms (HABs) and benthic hypoxia caused by organic matter accumulation , which threaten ecosystem services and cause great economic loss (e.g., Alorda-Kleinglass et al., 2021).
Estimates of nutrient loading driven by SGD have often been obtained by extrapolating the composition of the terrestrial groundwater endmember to the outflowing SGD mixture (e.g., Jiang et al., 2021;Lee and Kim, 2007;Taniguchi et al., 2008;Wu et al., 2013). For nonconservative species, such as nitrate, this approach is difficult to justify. It is generally taken because nitrogen concentrations may be very significant in fresh groundwater due to anthropogenic pollution, hence posing a significant threat to coastal ecosystems (Bowen et al., 2007;Valiela et al., 1990;Zhang et al., 2020). However, in the 'real world', the composition of SGD is modulated by the biogeochemical history of the water flows associated on the one hand with the transit path followed from land to the STE, and on the other to the biogeochemical reactions occurring within the STE, where fresh groundwater mixes with saltwater recycled through the coastal aquifer (Moore and Joye, 2021;Robinson et al., 2018;Santos et al., 2009a;Charette andSholkovitz, 2006, 2002). How the combined effects of these factors regulate SGD composition and impact ecosystem metabolism has not been shown.
The Ria Formosa in Southern Portugal is a leaky lagoon system with persistent signals of functional eutrophication, the drivers of which are still debated after decades of study (Domingues et al., 2017;Barbosa, 2010;Newton et al., 2003). SGD within the system was first detected and described by Leote et al. (2008). Almost a decade later, Rocha et al. (2016) showed how SGD inputs to the lagoon could be quantitatively separated into fresh groundwater inputs and saline recirculation flows using a combination of environmental tracers, including stable isotopes in water (δ 2 H, δ 18 O), radon ( 222 Rn) and salinity. The study provided evidence for the episodic nature of fresh groundwater inputs into the lagoon, but showed that seawater recirculation occurs permanently and drives an estimated gross load of~350 Ton N y −1 into the system. Using direct discharge measurements, combined with 222 Rn budgets where endmember attribution was achieved by water stable isotope signatures, Rocha et al., 2016 also showed that fresh, land-borne SGD could add a further~61 Ton N y −1 to the lagoon and was therefore capable of driving new production and/or shifting species composition. Hugman et al. (2017) estimated similar nitrogen loading (350 Ton N y −1 ) into the lagoon by SGD using numerical flow and transport modeling.
Despite the large number of studies carried out all over the world , coastal primary production dynamics have yet to be mechanistically linked to the variability in annual SGD composition and its hydrological and biogeochemical drivers. This causal chain has been elusive partially because changeable mixing of marine and continental waters in the coastal aquifer generates significant biogeochemical dynamics in the STE, resulting in SGD with highly variable composition throughout any one year (Ibánhez and Rocha, 2017;Jiang et al., 2018;Ibánhez et al., 2021). However, SGD surveys generally follow 'seasonal' schedules (i.e., summer, winter), are biased toward the best weather periods, or focus on either the hydrology, the biogeochemistry, or the ecology, paired combinations of these, and very seldom on all. This results in a lack of sufficiently detailed datasets to tackle the problem at the necessary level of granularity. The extensive data collected at high temporal resolution covering both annual and interannual (2005)(2006)(2007)(2009)(2010)(2011) timescales therefore make the Ria Formosa an ideal case study to determine the functional drivers of nutrient stoichiometry in SGD.
Here, we take advantage of this globally unique dataset to test a hypothesis for the first time in the real world: Can the composition of submarine groundwater discharge (SGD) to a lagoon be regulated by the annual oscillation of the local STE, itself driven by groundwater recharge variability and, in turn, are the resultant compositional dynamics reflected in ecosystem function?
To test this hypothesis, we revisit our Ria Formosa dataset and first examine the temporal dynamics of salinity, E H , pH, dissolved oxygen, nutrient concentrations (N, P, and Si) and stoichiometry in SGD flows; we then assess the annual variability of SGD composition and SGDdriven nutrient fluxes in the context of the local hydrological regime; finally, we compare these to synoptic benthic primary production measurements to ultimately understand the impact of SGD on primary production over the course of a year. Our objectives are: (1) to determine the annual variability of nutrient content and N:P:Si stoichiometry in SGD into the lagoon; (2) to identify the main drivers of variability in SGD composition and stoichiometry and their interactive effects; (3) to discuss links to, and implications for, ecosystem function that could help define expectations of cause-effect relationships between SGD and primary production dynamics for environmental management of the lagoon and similar systems elsewhere.

General characterization
The Ria Formosa (Fig. 1) extends from 36°58′N 8°20′W to 37°30′N 7°32′W in the south of Portugal. It is a shallow (average depth~2 m) coastal lagoon separated from the Atlantic Ocean by a chain of five barrier islands and two sandy spits. The system was formed by rapid infilling of palaeovalleys during the early Holocene marine transgression (10000-8000 cal yrs. BP), followed by back barrier infilling from 7250 BP (Sousa et al., 2019). The regional climate is semi-arid, with an average annual temperature of 17°C, cooler winters (mean 11°C) and warmer summers (mean 24°C). Mean annual precipitation is~480 mm, but the watershed (740 km 2 ) only receives effective precipitation of 152 mm year −1 due to high potential evapotranspiration (Salles, 2001). Except for the River Gilão ( Fig. 1), which intermittently discharges almost directly into the Atlantic at the eastern limits of the system, the other four minor rivers and fourteen streams discharging into the lagoon are ephemeral and dry out in spring/summer.

Hydrogeology
The lagoon is connected to the Campina de Faro coastal aquifer system (M12, c.f., Fig. 1,~86.4 km 2 ), composed by three distinct units (Almeida et al., 2000). The oldest and deepest is a cretaceous aquifer, formed by limestone, marls and dolomites resting on Jurassic sedimentary formations. A sandy limestone Miocene aquifer overlays this cretaceous unit, thickening from north to south and reaching a depth of >200 m close to the littoral zone (Geirnaert et al., 1982). On the surface lies an unconfined (phreatic) aquifer composed of detritic sediments (fine sands), sandy limestones and sandstones from the Miocene, Pliocene sands and gravels, alluvial deposits, and sand dunes from the Quaternary. These deposits reach an average depth of 20-30 m across the system (Almeida and Silva, 1987) with a Fig. 1. Location of the sampling points within the wider geographical setting of the Ria Formosa lagoon (top, inset), including meteorological stations from which data was gathered for this study, SNIRH stations São Brás de Alportel (31 J/01C), Estói (31 J/04UG) and Quelfes (31 K/02UG); the lower inset shows the seepage meter array site in the Ancão peninsula, SNIRH monitoring wells and boreholes Ludo (610/167), Vale das Almas (610/6), Areal Gordo (611/217) and Horta do Barrote/Pechão (611/91) screened in the unconfined Campina de Faro aquifer system (M12), as well as the site of monthly microphytobenthos (MPB) Chl a measurement by Brito et al. (2009). maximum thickness of~60 m in some areas north of Faro and connect with modern sand dunes and alluvial deposits within the lagoon (Silva, 1988). The northern limit of the aquifer system lies against outcrops of permeable and dolomitic limestones of the Jurassic and cretaceous formations. The maximum length (measured N-S) is~6.6 km, found north of our seepage meter deployment station (Fig. 1), while the maximum width (E-W) is 24 km. Recharge occurs directly by precipitation (Almeida et al., 2000;Silva et al., 1986). Infiltration may occur over the land surface or through the Jurassic and cretaceous outcrops at the northern limit of the aquifer unit (Almeida and Silva, 1987). Engelen and van Beers (1986) propose that the intermediate Miocene unit is confined and discharges into the Atlantic Ocean, bypassing the lagoon. There is substantial evidence, however, that it is hydraulically connected to the surface unconfined aquifer. This connection is supported by co-contamination of the deeper aquifer by nitrate (NO 3 − ) pollution derived from ammonium sulphate fertilizers employed in agriculture on the coastal plain (Almeida and Silva, 1987), and by stable isotope signatures of groundwater taken from both units .

SGD to the Ria Formosa
SGD comprising significant freshwater contributions was first detected in the Ria Formosa in 2006-2007 and described as a potentially important source of nutrients, in particular nitrogen, to the lagoon (Leote et al., 2008). Lobo-Ferreira et al. (2007) Brzezinski, 1985;Brzezinski et al., 2003).

Sampling and chemical analysis
Monthly sampling campaigns were conducted on a beach at the inner side of the Ancão peninsula throughout 2006 (Fig. 1). The tide is semidiurnal with average amplitudes of 2.8 m for spring tides and 1.3 m for neap tides. The local sediment is classified as medium sand, with a median grain size (d 50 ) of 0.44 mm, a clay+silt content ranging from 0.8 % to 1.6 %, mean porosity of 0.35, and an hydraulic conductivity of 5.5 × 10 −3 cm s −1 , or~4.8 m d −1 (Rocha et al., 2009). Details on the sampling strategy, methods, and nature of the data gathered during these field surveys may be found in Leote et al. (2008), Rocha et al. (2009) and Ibánhez et al. (2013Ibánhez et al. ( , 2011. Briefly, four to seven seepage meters (Lee, 1977) with attached water collection bags (two per meter) were deployed along the beach profile each month at spring low tide and allowed to settle in for 24 h prior to sampling. Water samples were then collected at subhourly intervals, following the precautions listed in Cable et al. (1997). Sampling covered a minimum of two and a maximum of four successive tidal cycles every month. Discharge rates measured in-situ (annual range: 1.7-190.6 cm d −1 ) were intrinsically related to tidal height, with peak discharge observed at low tide (Leote et al., 2008;Rocha et al., 2009).
Seepage water was analysed in-situ for salinity, oxidation-reduction potential (4 M Ag/AgCl reference electrode), oxygen concentration, and pH with an YSI 600 multi-parameter probe (Yellowspring Instruments®). Probe sensors where individually calibrated following manufacturer instructions prior and again after each survey to check for potential drift. For nutrient analysis, water samples were filtered through polyethersulphone (PES) membrane filters (Rhizon SMS, Eijkelkamp Agriresearch Equipment®; pore diameter of 0.1 μm) directly into nonadditivated sterile vacutainers. The combination of ultra-filtration through a hydrophilic membrane directly into a sterile environment in isolation from the atmosphere and environmental contamination vectors preserves nutrient concentrations with a high degree of accuracy, while substantially mitigating the risk of cross contamination or bacterial adulteration of the analytes (Seeberg-Elverfeldt et al., 2005;Jiang et al., 2017). The PES samplers were intensively flushed with sample water prior to the collection of the filtrate to avoid NH 4 + sorption artifacts as described in Ibánhez and Rocha (2014). Concentrations of nitrite (NO 2 − ), ammonium (NH 4 + ), SRP and DRSi were determined following the spectrophotometric methods described by Grasshoff et al. (1999) within 48 h of collection. The NO 3 − concentration was quantified by the cadmium reduction approach of Jones (1984). The method detection limits (MDL, in μM) for the determination of each analyte, calculated following IUPAC recommendations (Analytical Methods Committee, 1987;Long and Winefordner, 1983) were: 0.08 ( Measurements of ORP and pH (NBS) were temperature compensated. For each sample, the pH measured in-situ (NBS) was corrected analytically for ionic strength, using the Davies extension of Debye-Hückel theory on the behaviour of electrolyte solutions (Davies, 1938). Results are expressed in the total hydrogen ion concentration scale (Dickson, 1993). The median of each monthly set of values measured on separate SGD samples (26 < n < 55) was used to characterize the central tendency of the distribution of pH, O 2 concentration and ORP measured in SGD each month, and the median absolute deviation (MAD) employed as a robust estimator of associated uncertainty (Hampel, 1974).

Calculation of SGD-driven nutrient fluxes
SGD-driven nutrient fluxes were quantified by multiplying the nutrient content of water collected from the seepage meters by concurrent measurements of volumetric discharge. For each seepage meter, the nutrient flux (F, mmol m −2 d −1 ) was calculated as: where Nu is the total number of samples collected from each seepage meter deployed in any monthly field campaign, C i is the nutrient concentration in a specific seepage water sample, V i is the volume of the seepage water sample, T is the sample collection timespan for each volume and A the surface area covered by the seepage meter. Subsequently, nutrient fluxes derived from different seepage meters located at identical distances from the lowwater mark (LWM) along the cross-shore transect were averaged. The standard deviation of the average was employed to perform end-to-end uncertainty analysis, carried out following general rules for the propagation of unknown uncertainties, taken here as the standard deviation (± 1σ) of an estimated mean value (Taylor, 1997). The total daily SGD nutrient flux (mmol m −1 d −1 ) was then obtained by spatial integration of the mean nutrient fluxes derived as above for different positions with respect to the LWM along the cross-shore transect. These were then normalized to the length scale covered by the seepage meter array to obtain areal fluxes (mmol m −2 d −1 ). The fresh and saline volumetric fractions of SGD (fSGD and sSGD, respectively) were estimated using a linear salinity-mixing model. In the breakdown of seepage into fSGD and sSGD, the seawater salinity end member was taken as the average of the salinities above the 97.5th percentile measured in outflow, rather than the highest salinity observed within the chambers during deployment to avoid an artificial 'freshening' of the fractional composition caused by the effect of evaporation on porewater salinity (Geng et al., 2016;Santos et al., 2009b).

Groundwater recharge and inland water table fluctuations
Meteorological data for the region and water table levels for the superficial unconfined aquifer of the M12 system were taken from the online Assessing the accuracy of any recharge quantification method is difficult. Because we sought an indicator of actual recharge throughout the year, and not a precise determination, we prudently estimated 'potential' aquifer recharge rates (R P ) by following two different approaches. In the first, we used monthly water table heights for each borehole to estimate potential recharge by the water table fluctuation (WTF) method (Healy and Cook, 2002), as: Here, Sy is the specific yield for the aquifer, and dh/dt the temporal variation between water table height at the high points during the year and the antecedent height, obtained by extrapolation of the recession curve. For this study, and since the Campina de Faro aquifer is mainly composed of detritic sediments, including fine sands, we took the specific yield Sy to be constant during the year and equal to 0.25 (Johnson, 1967). The resultant potential recharge rates are plotted in Fig. 2 for the period 2005-2008. In parallel, potential monthly recharge for the year of study (2006) was estimated from the difference between the mean monthly precipitation estimated from the three meteorological stations covering the coastal plain ( Fig. 1) and local potential evapotranspiration estimated according to Thornthwaite's formula (Thornthwaite and Mather, 1957). The annual variation of potential recharge obtained this way compared favourably with the long term annual monthly means (Fig. 3a), determined from concurrent data of precipitation and pan evaporation measured from 1980 to 2000 in São Brás de Alportel (31 J/01C).

Transgressive oscillation of the STE
The annual variability of SGD rates (fresh and saline) at the STE mouth depends directly on the change in hydraulic head driving land-sea water flow and indirectly on the dynamics of the STE (Michael et al., 2005;Robinson et al., 2006). Here, we describe the amplitude of oscillations of the STE throughout the year with reference to the toe of the saltwater wedge, employing the Glover solution to delimit its position (Glover, 1959). This assumes a sharp saltwater-freshwater interface, dynamic equilibrium between fresh and salty groundwater following the Ghyben-Herzberg relationship, a stagnant saltwater mass, isotropic media and an aquifer body that is confined at the top and of infinite depth. These assumptions do not describe our aquifer closely, since it is unconfined, has average thickness of~30 m and sits on an impermeable layer (see Section 2.2). Hence, we follow the re-interpretation of the Glover solution for the case of an unconfined surface with a positive freshwater head gradient toward the sea, a horizontal seepage face and an impermeable layer at depth that intersects the saltwater-freshwater interface to form a saltwater wedge that extends some distance inland (Cheng and Ouazar, 1999). This approach finds only small discrepancies between the saltwater toe positioning given by the Glover formulae for the two conceptual models of aquifer.
Because the flow is unconfined, and the slope of the water table is small (in our case, L, the cross-section length of the aquifer, is measured in kilometres for an average aquifer depth, D, of just 30 m, implying L >> D), we may take the Dupuit-Forchheimer approximation of Darcy's law (Dupuit, 1863;Forchheimer, 1930) to simplify by one dimension and treat SGD as a 2D flow problem. This approximation relies on the following assumptions: a) Flow is horizontal over any vertical cross-section of the aquifer b) The velocity is constant over the depth of the aquifer. c) To calculate the velocity, we employ the slope of the free surface (dh/dL) as the hydraulic gradient d) The slope of the water table is relatively small, i.e., Δh << L. This is reasonable to accept for our case, since the amplitude of water table fluctuation ranges from 0.29 m to 1.87 msee Table 1 for a cross-sectional length of the aquifer ranging between 2300 and 12000 m, respectively, see Fig. 1).
The formulation of Darcy's law with the Dupuit-Forchheimer approximations to calculate the mean discharge to the sea per unit length of aquifer shoreline Q fw (m 3 m −1 d −1 ), for a cross section length of the aquifer L and a mean saturated aquifer depth D is (e.g., Cushman and Tartakovsky, 2016): where K is the hydraulic conductivity (4.8 m d −1 ), the mean depth D of the aquifer is 30 m and dh/dL is the hydraulic gradient (m m −1 ). For the latter, we are interested in the component of the flow vector directed southward ( Fig. 1), so we use the north-south distance between each of the four piezometric stations and our seepage meter array location for L (note that L is measured landward from the point of seepage (x = 0) to the piezometric stations inland (x = −2300 to x = −12,000 m), so that L < 0 and there is no negative sign appearing in Eq. (3)), while dh is taken monthly from the piezometric time-series at each station. The average fresh groundwater discharge from the M12 unconfined aquifer to the lagoon calculated this way was 0.126 ± 0.017 m 3 m −1 d −1 for 2006, comparing very favourably to direct measurement by our seepage meter array (0.115 ± 0.079 m 3 m −1 d −1 ). We are therefore confident that the saltwater-freshwater interface position in this aquifer is quite well described by the variation on the Glover approach developed by Cheng and Ouazar (1999) and that our seepage meter array produces flow magnitude data that are consistent with theoretical predictions based on continental water table fluctuations. In order to delineate the position of the saltwater wedge toe, and hence the STE, we first calculate the horizontal gap through which freshwater flows through the seafloor to the coastal zone, S g,fw , as: where γ is the excess of the specific gravity of seawater over that of freshwater, given by (ρ sw -ρ fw )/ρ sw , where ρ fw and ρ sw are the densities of freshwater and saltwater, respectively taken as 1000 and 1025 kg m −3 . The position of the saltwater wedge toe, Sw L , can then be estimated to be at a distance from the coastal boundary (positive inland) where the freshwater seepage gap develops (L = 0) as: We then calculate the monthly position of the saltwater wedge toe in relation to all 4 piezometers and express it with reference to its location when the water table height measured at each piezometer is equal to the annual mean. This method yields a curve that illustrates the annual oscillation of the saltwater wedge around a reference point for every piezometer site, with positive values indicating saltwater intrusion and negative values indicating saltwater retreat. We then normalize the extension of that oscillation to the largest departure in length, either inland or seaward, as measured from the central reference point for each annual curve. The result yields normalized intrusion and recession curves oscillating between a maximum of +1 (inland) and a minimum of −1 (seaward), which can be compared to water table height oscillations measured at the respective piezometer sites treated in the same way ( Fig. 3c).

Spatial heterogeneity in STE outflow composition
High spatial and temporal variability of seepage outflow rate and composition determinations are to be expected when SGD is monitored by seepage meters (Taniguchi et al., 2019). Out of an abundance of caution, we therefore approach the treatment of nutrient composition in the STE outflow for the purpose of studying relationships with primary production in two different ways. Apart from the monthly nutrient fluxes, which provide an integrative measure of the spatial and temporal variability of individual measurements, we also evaluate the nutrient stoichiometry of monthly sample sets using frequency distribution analysis. This approach allows us to assess whether the uniformization of data resulting from the spatial and temporal integration process necessary to obtain monthly fluxes does not obscure natural variability patterns.
To this end, we investigated the frequency distribution of samples in which any of the nutrient species was impoverished in relative terms ('nutrient-limited', as it were), i.e., with a concentration below the R-B ratio, and how the frequency of this scenario varied with the fraction of freshwater contained in SGD along the year. To do so, the salinity range covered by the water sample pool (n = 391) was first divided into classes (value intervals), with a minimum of 20 samples per class. This ensured uniformity of the distribution, balancing the requirements of both interval quantity and sample size within each range for subsequent frequency analysis. For each class, the relative frequency with which an individual element (N, P, Si) fell below the idealized R-B ratio (N:P = 15, Si:N = 1, and P:Si (or N) = 0.07) was calculated by: where F N is the relative frequency of a putative nutrient-poor sample occurring within a salinity range-defined class, and N N , N Si and N P represent the frequency with which N-poor, Si-poor and P-poor samples were found within the same salinity range. The same process was followed in turn for Si and P. Subsequently, the discrete frequencies of relative nutrient depletion found for each salinity interval were fitted by Dirichlet regression (Aitchison, 1982) to model the frequency of nutrient-depleted samples per chemical element over entire salinity range covered by the sample pool. To fit the data, a system of quadratic regression equations, defined as follows, was applied: where F N S , F DSi S , and F P S are the model-fitted frequencies for N, Si and P being the limiting nutrient, S is salinity, and A, B, and C are parameters generated by the best fit of the quadratic regression system to the data. The system of equations was solved by the approach minimizing the standard error (Nelder and Mead, 1965) in MATLAB®.

Microphytopbenthos production (as Chlorophyl a)
Benthic chlorophyll (Chl) a was measured during 2006 close to our seepage meter array on the same sand spit (Fig. 1) by Brito et al. (2009). To evaluate the possible link between monthly SGD-driven nutrient fluxes and microphytobenthos (MPB) production, we used logistic regression (Verhulst, 1938, in Cramer, 2002. We model the measured MPB (in μg Chl a g −1 sediment) as a function of monthly SGD nutrient flux (F), using the logistic growth expression where MPB(L) is the carrying capacity, k the rate of growth, and b a constant. Boundary conditions were MPB = MPB(L) when F → ∞ and The model was fit to MPB data as a function of monthly SGD-driven nutrient fluxes (DIN, SRP and DRSi) determined in our study by minimizing the sum of squares of differences between model output and the data of Brito et al. (2009) while varying k, using the Microsoft™ Solver routine in Excel. Initial and carrying capacity MPB values (MPB(L 0 ) and MPB(L)) where first fixed to the annual minimum and maximum MPB found on site by Brito et al. (2009) of~2 and~20 μg Chl a g −1 , until k was found, and then allowed to fluctuate within the bounds established by the standard deviation of the mean of determinations throughout the year (±10 %). Goodness of fit was assessed by regression analysis comparing model outputs and MPB data at the 95 % confidence level.

Precipitation, recharge, and STE oscillation
The twenty year-long datasets  collected at S. Brás de Alportel (31 J/01C), suggest monthly precipitation over the Ria Formosa peaks in December (20-year mean of 211 ± 194 mm) or November (idem, 149 ± 132 mm) and is lowest in August (2.4 ± 3.3 mm). The (a) Annual variability of potential recharge estimated as the monthly effective precipitation, taken as the difference between precipitation (P) and pan evaporation (E) or evapotranspiration (ET). The dashed curve shows the 20-year historical means (± S.D.), obtained from the mean monthly differences between precipitation (P) and pan evaporation (E) measured at SNIRH station São Brás de Alportel (31 J/01C) for the period 1980-2000. The solid curve shows the regional P-ET estimate for 2006, with precipitation taken as the average of monthly rainfall (± S.D.) measured at SNIRH stations São Brás de Alportel (31 J/01C), Estói (31 J/04UG) and Quelfes (31 K/02UG) and evapotranspiration estimated monthly using Thornthwaite's expressions (Thornthwaite and Mather, 1957). (b) Time series analysis of the water level in the four wells screened in the Campina de Faro surface aquifer for the period January 2005 to December 2008, revealing the degree of correlation between lagged monthly water table height for the four-year time series of piezometric head in all stations, using an autocorrelation function (ACF) model. Highly significant (p < 0.001) correlation between water table heights within the aquifer are found at seasonal lags of 6 and 12 months, with negative coefficient for the former and positive for the latter. (c) Annual (2006)  long-term average pan evaporation rates peak in July (257 ± 25 mm) and are lowest in January (45.7 ± 11.9 mm). During the four-year period that frames our SGD dataset (2005)(2006)(2007)(2008), the average precipitation recorded by the three meteorological stations in the watershed (Fig. 1) (Fig. 2). Recharge estimated by the water table fluctuation method followed the annual precipitation pattern (Fig. 2), but the more inland stations (610/167-Ludo and 611/91 Horta do Barrote/Pechão) registered higher amplitude fluctuations on water level (Table 1) than the southernmost stations (610/6 and 611/217). This is reflected for example on the recharge rates calculated for November 2006, which reached 38 and 53 cm month −1 respectively at 610/167 and 611/91, while only 7.8 and 2.8 cm month −1 were found respectively for 610/6 and 611/217. The closer proximity of 610/167 and 611/91 to the boundaries of adjacent regional aquifers (M11, Quinta João de Ourém, and M10, São João da Venda/Quelfes) suggest that recharge of the aquifer may also occur by inflows from these, as well as through more rapid infiltration through the Jurassic and Cretaceous outcrops at its northern limit, as suggested by Almeida and Silva (1987). Overall, the data support an overall North-South direction of groundwater flow and fast recharge of the M12 aquifer during the wetter winter months (November-January) followed by a slower fall of the water table throughout the Spring-Summer period, up to October (Fig. 2).
The comparison between historical P-ET data obtained in S. Brás de Alportel with that obtained for 2006 using Thornthwaite's expression for evapotranspiration ( Fig. 3a) suggests that, normally, direct recharge by infiltration through the surface soil would only reliably occur during the first and last two months of the year (January and February, November, and December). However, in 2006 this period of effective precipitation could include March and October while December was quite dry (~30 mm) in historical terms. Apart from establishing historical relevance, though, we needed to investigate whether the behaviour of the water table during 2006 was representative of the normal aquifer dynamic, and this way speak of seasonality of STE oscillation, rather than random variability. To establish this point, we measure the correlation between lagged monthly water table height for the four-year time series (2005)(2006)(2007)(2008) of piezometric head in all stations, using an autocorrelation function (ACF) model (Fig. 3b). Highly significant (p < 0.001) correlation between water table heights within the aquifer are found at seasonal lags of 6 and 12 months, with negative coefficient for the former and positive for the latter (Fig. 3b). This shows that groundwater table falls to its lowest levels in June/July and rises to their highest in November/December every year. This dynamic is consistent with the potential recharge data (Figs. 2; 3a) and indicates a strong water balance cyclicity dictated by the seasonality of precipitation over the basin. This strong seasonal signal on water table fluctuation implies that the saltwater-freshwater interface also shifts, expanding and contracting within the coastal aquifer, establishing a cyclic dynamic for the STE (Fig. 3c). During 2006, this oscillation resulted in a phase of saltwater wedge retreat (SWR) seaward until March/April, coinciding with freshwater flushing of the coastal aquifer borne from high recharge. This phase was followed by landward intrusion (SWI) until October driven by the continuous drop in water table height before the SWR phase restarted from November (Fig. 3c). The oscillating behaviour of the saltwater wedge driven by the annual continental recharge pattern can perhaps be better visualized in Fig. 3d. The plot shows how the transgressive migration of the STE occurs on an annual cycle (an oscillation) that is consistent both with historical data on effective precipitation and with P-ET estimates for the year 2006.

Saltwater wedge oscillation and pH of STE outflows
The transgressive movement of the STE through the two phases (retreat and intrusion, expansion, and contraction) is reflected on the dynamics of pH measured in STE outflow at our seepage meter array (Fig. 4). Monthly effective precipitation measures correlate well with monthly pH measured in STE outflow, with higher pH observed during the SWI phase and lower during the SWR stage, if either historical (Fig. 4a) or contemporary (Fig. 4b) data for effective precipitation is used. During the period in which the saltwater toe is located seaward of the fulcrum point (January-June, Fig. 3c, d) the STE outflow pH follows a different dependence on effective precipitation than that observed during the SWI period (July-December) but for the historically anomalous December of 2006 (Fig. 4a). However, this dichotomous dependence of pH on effective precipitation is re-established if we compare the 2006 pH data with P-ET estimates for the same year (Fig. 4b). In addition, the antithetical behaviour of STE outflow pH and continental recharge is quite clear when both are plotted monthly for 2006 (Fig. 4c, inset). Combined, our observations suggest a dynamic, but strong, connection between coastal aquifer recharge (inferred from effective precipitation) on the coastal plain and the chemistry of the STE outflows, over yearly time scales.
The cyclical nature of the annual variability of median pH in STE outflows is illustrated as a function of its transgressive oscillation in Fig. 5a. Note the sudden increase of pH in December is linked to a fall in effective precipitation and therefore recharge for that month by comparison to historical levels (Fig. 4b), implying that pH increases result from contraction of the STE inland and the opposite occurs during saltwater wedge retreat. However, the same dynamic is not reflected on redox potential (E H ). This has important geochemical consequences. Fig. 5b shows the Pourbaix diagram for iron in our system, with the pH-pE pairs (median ± MAD of all our monthly data) superimposed. The results highlight a linear dependence of redox potential on pH, with higher pE/pH ratios measured during the freshwater flushing stage of the STE excursion (SWR stage, Fig. 3c, d,  Fig. 4), but indicate that pE of STE outflows changes little throughout the year. Assuming rapid equilibration of the iron system, particularly with regards to the stability of amorphous iron hydroxides, the diagram also indicates that within the STE, it is very unlikely that P would be in free-solute form because of the absence of iron adsorbent phases, with possible exception of January, where the median pH of the solution drops below the conditions for thermodynamic stability of Fe(OH) 3 and favour predominance of Fe in its reduced form (Fe 2+ ). However, this is also when the mean redox potential of the outflow is the highest. This is consistent with oxygen data, which show oxygenation levels in STE outflows range between 81.7 and 36.4 (±3) percent saturation, measured in January and December respectively. Regardless, the pH, but not the E H , of the STE outflow is affected by continental recharge cycles, via the effect they have on the dynamics of the saltwater-freshwater interface within the coastal aquifer (Fig. 5a).

Magnitude and variability of SGD rates in 2006
Total SGD (tSGD = saline + fresh groundwater discharge) rates had a bimodal dependence on precipitation throughout the year (Fig. 6a). Perhaps counterintuitively, the highest tSGD rates (505 ± 157 L m −2 d −1 ) during 2006 were measured in July, for a monthly precipitation of 0.9 mm. Generally, tSGD rates decreased from the driest months to a minimum in February (68 ± 9.5 L m −2 d −1 ) for a monthly precipitation of 58 ± 2 mm, and then increased again into the wetter months to November, when tSGD = 377 ± 87 L m −2 d −1 for a precipitation rate of 194 ± 21 mm. This data confirms that factors other than precipitation alone drive SGD rates through the coastal aquifer into the lagoon throughout the year. In June and July (Fig. 6b, c), SGD into the Ria Formosa contained mostly recirculated seawater (freshwater contribution <2 % of tSGD). The spring-summer period split the year in two: earlier, rates of fresh groundwater discharge increase from January to March (6.1 to 29.4 L m −2 d −1 and from 8.2 to 19.6 % of tSGD), and later progressively from June (3.9 L m −2 d −1 and 2 % of tSGD) to November (36.3 L m −2 d −1 and 10 % of tSGD). Both fresh groundwater seepage rates (fSGD, Fig. 6c) and relative freshwater content in SGD (Fig. 6b) reach annual peaks in March and November, while annual lows are observed during most spring-summer months (April to August). Fresh groundwater discharge in December dropped significantly from November (fSGD = 8.2 L m −2 d −1 , 6 % of tSGD) as did rainfall (Fig. 2).
The change in monthly fSGD rates throughout the year was consistent with mean recharge dynamics estimated for the aquifer by the water table fluctuation method (Fig. 2). Higher fSGD rates (March, November) also coincide with months of positive effective precipitation (Fig. 3a). It is interesting to note also that the highest SGD rates, either for fresh groundwater (March, November) or salty groundwater (July) happen at the approximate time the STE reaches its maximum excursion seaward or incursion landward (fresh SGD) or the fulcrum point in between, when peak saline SGD rates are measured (Figs. 3c, d; 6b, c).

Nutrient composition in SGD
The salinity and NO 3 − concentration of the seeping water varied significantly during the year (Fig. 7a). Nitrate concentrations in SGD ranged between 180.3 μM and the minimum detectable concentration, correlating significantly with salinity. The mixing plots show the data fall into two different subsets corresponding respectively to a 'fall' term, with higher slope of the mixing line (R1, August to December), and a 'spring' term, with lower slope of the mixing line (R2, January to July). On the other hand, the concentrations of both NO 2 − and NH 4 + in SGD were much lower than those of NO 3 − and showed no obvious conservative behaviour along the sample salinity range (Fig. 7b, c) (Fig. 7d). The intercepts (S = 0) for DIN (354 ± 15 and 99 ± 3.5 μM, respectively during the "fall" and "spring" terms) are much lower than the estimated fresh groundwater endmember (2.1 mM, see Section 2.3), indicating that N is not conserved along the entire extent of the groundwater-seawater mixing line.
The DRSi concentrations in SGD varied between 0.3 μM and 55.3 μM (Fig. 7f) and were not consistently conservative within the salinity range throughout the year, with the most obvious outliers being the colder months of January and December. When these two months are excluded, the remaining data fall into a significant linear correlation (df = 336, Fig. 5. Annual variation of the monthly pH of STE outflows and its geochemical impact (a) cyclic nature of the annual change of STE outflow pH as a function of salt wedge toe oscillation: note December 2006 position results from a contraction of the STE resulting from an unusually dry month in historical terms; (b) Pourbaix diagram for the Fe-CO 2 -H 2 O system in the STE with monthly pE-pH pairs overlain (C T (total inorganic carbon) = 1.6 × 10 −3 M, soluble Fe species = 1 × 10 −9 M. Solid phases (s) are Fe(OH) 3 (amorphous), Fe(OH) 2 (s)). The pH of STE outflow again oscillates within tight physicochemical bounds, with the highest pH corresponding to August and the lowest to January. r 2 = 0.462, p = 0), with intercept at zero salinity of 69 ± 3 μM, statistically indistinguishable from the mean groundwater concentration of 65 ± 7 μM measured from 2010 to 2013 (Section 2.3). This suggests a continental (fresh) groundwater origin for the silicate arriving at the beach aquifer, and since Si diagenesis is not affected by redox reactions, some additional silica mobilization in transit to explain the January and December data. Indeed, the highest DRSi concentrations could be measured in both the more saline (S > 36, during January and December) and the more brackish SGD samples (S < 20). Combined with the convex shapes of the mixing curves in January and December (Fig. 7f), this suggests the beach aquifer could be a net contributor of DRSi to SGD during the colder months. Dissolved P was low throughout the year, with maximum measured SRP concentrations of 8 μM and >40 water samples (> 10 % of the sample pool) below minimum detectable concentrations (Fig. 7e). This is consistent with the pE-pH stability data (Fig. 5b).

Seasonal variability in SGD nutrient loading and stoichiometry
Monthly DRSi fluxes, driven by SGD (Fig. 8a) varied between a minimum of 6.0 ± 0.4 mmol m −1 d −1 in June and a maximum of 40.6 ± 8.9 mmol m −1 d −1 in March, correlating well with monthly fSGD rates (Fig. 6c, Students t-test, p < 0.02). This is consistent with earlier data indicating that Si loading is controlled mostly by the fresh groundwater endmember (Fig. 7f). NO 3 − fluxes were significantly higher during the period of August to November (Fig. 8b)  SRP fluxes (Fig. 8c) measured from July to December were on average 4 times higher than those determined between January and June (means of 0.4 and 1.65 mmol P m −1 d −1 , respectively). Overall, the salient features of the annual variability of SGD-driven nutrient fluxes were the difference between the first part of the year (until June) and the second (July onwards) for SRP and DIN, the relatively low SGD-driven nutrient fluxes in July, when the saline component of SGD is highest (with SRP being the exception), and the monthly NO 2 − SGD flux and the Si:P SGD flux quotient being significantly correlated with the fSGD/sSGD ratio (Students t-test, p < 0.05 and p < 0.01 respectively). As a result, SGD contributed to enrichment of N over P from February to May and from August to November (Fig. 8d) and over Si in February, May, and August to November (Fig. 8e) while in July a deficit of N and Si relative to P was observed.

SGD-driven nutrient fluxes and primary production
In 2006, the mean water column Chl a was lowest in July (0.37 mg m −3 ) and peaked in April (4.1 mg m −3 ) and again in September (4.2 mg m −3 ) (Brito et al., 2012). The MPB Chl a, measured over the period 2006-2008 at the sand spit station (Station 'Ponte' in the original paper) every fortnight from April to November ranged from 2 to 25 μg Chl a g −1 (Brito et al., 2009). For 2006, MPB Chl a rose steadily from 8.7 μg g −1 (dry sediment) in April to 18.2 μg g −1 in August, dropped to 17.1 μg g −1 in September and increased again after that to a peak in November (22.3 μg g −1 ). For the period 2007-2008, MPB Chl a peaked above 15 μg g −1 in April July and December 2007, and January 2008.
The microphytobenthic production could be modelled as a logistic function of monthly SGD-driven fluxes for DRSi, DIN and SRP (Fig. 9a, b, c) in a statistically significant way (p < 0.0005, p < 0.02, p < 0.005, respectively with df = 6). Starting MPB levels (at SGD nutrient flux = 0) were found to fit the growth curve well at 2.65 μg Chl a g −1 for all the nutrients. Carrying capacity (μg Chl a g −1 MPB) was 20.93 for DRSi, 19.46 for DIN, and 18.71 for SRP, fluctuating around the mean peak MPB of~20 μg Chl a g −1 found in-situ and well within the maximum within-survey sample variability of ±5 μg Chl a g −1 reported by Brito et al. (2009). The best-fit growth rate constant k was 1.86 d −1 (±0.24 S.E.) for DRSi, 1.59 d −1 (±0.30 S.E.) for DIN, and 40.6 d −1 (±8 S.E.) for SRP. The proportions between nutrient-specific growth rate constants, taken as the reciprocals (i.e., for N:P, k P /k N ) reflected the ideal R-B ratios relatively well: 0.9 ± 0.2 for k N :k Si (N:Si~1), 22 ± 5 for k P :k Si (Si:P~15), and 26 ± 7 for k P : k N (N:P~16), consistent with MPB production being regulated by nutrient resource ratios in SGD and P being most often the limiting nutrient (Fig. 10c, d).

Nutrient stoichiometry in SGD
The frequency with which each individual nutrient in SGD water deviates the most toward scarcity from the idealized R-B ratio (i.e., is the most depleted with regards to N:P:Si = 16:1:15) is plotted as a function of the freshwater portion in the seeping water in Fig. 10, separating the 'spring' from the 'fall' term on the basis of the behaviour of NO 3 − with salinity ( Fig. 7a), as well as the putative location of the salt wedge toe with regards to its equilibrium position in the year (Fig. 3c, d). The 'fall' period would correspond to saline intrusion and expansion of the STE inland and toward the surface, while the 'spring' term corresponds to saltwater retreat during the wetter months of the year and the ensuing contraction of the STE. The best fit for the Dirichlet regression model (Eqs. (7)- (10)) is then overlain on the various discrete histograms. All possible scenarios, i.e., where either N, P, or Si were the most depleted in individual SGD samples occur. However, the frequency with which each nutrient fell into the 'limited' category with both the freshwater portion of the SGD outflow mixture and with regards to whether samples were collected during the first (spring) or the latter (fall) part of the year. Because brackish STE water contains a significant concentration of NO 3 − (Fig. 7a), samples taken from January to July with >20 % of fresh groundwater always carry an excess of N over the other nutrients (Fig. 10a). This is reinforced during the latter part of the year (August-December), given that the amount of NO 3 − contained in the fresher groundwater (S~15) seeping out of the beach face is larger (Fig. 7a). During the latter part of the year a ratio of >18 % fresh groundwater to total SGD is enough for N to be always above the R-B ratio (Fig. 10b). The frequency of N-depletion in samples decreases rapidly with increasing freshwater fraction in SGD and this in turn determines the relative impoverishment of the complementary nutrients, P (Fig. 10c, d) and DRSi (Fig. 10e, f) at higher freshwater contents. From January to July, when some SGD samples contain the highest portion of fresh groundwater measured in-situ (>30 %), P is always scarce in comparison to N and DRSi (Fig. 10). Conversely, when SGD is composed exclusively of recirculated seawater (fresh groundwater portion = 0 %), N-depletion occurs in~60 % of the samples, followed by P-depletion (~30 %) and finally by DRSi depletion (Fig. 10a, c, e). In between these extremes, the frequency with which P-depletion occurs increases with mounting fresh groundwater contribution, rapidly until the makeup of the outflow is~18 % fresh groundwater, in parallel with the decrease in the frequency of N depletion and a slight increase in the frequency of DRSi-scarcity. Beyond the~18 % fresh groundwater content breakpoint, N-depletion never occurs, and DRSi-depletion occurs only in 4 % of the samples. From August to December, the frequency with which the compositional range of the samples falls below the R-B ratio is quite like that shown for the previous term with regards to N-scarcity, but very different for P-and Sidepletion (Fig. 10d, f). Water samples with high fresh groundwater fractions frequently had a lower Si:N ratio: Si-depletion (N:Si > 1) occurred iñ 60 % of the samples when SGD was composed of >30 % of fresh groundwater, while P-depletion (N:P > 16) occurred in~40 % of the samples, compared to respectively 0 % and 100 % during the earlier part of the year. However, the SGD fractional composition breakpoint of~18 % fresh groundwater was still observable with respect to the nutrient stoichiometry in the outflows. From 0 % freshwater content to this breakpoint, the occurrence of the different scarcity scenarios was quite like that observed during the earlier part of the year, but with fresh groundwater content increasing beyond this point, the frequency with which P-depletion occurred decreased while Si-depletion increased, in stark contrast to observations during the spring term. This is also observable in the monthly-integrated flux ratio of Si to P (Fig. 8f). This implies that the STE outflow composition in terms of NO 3 − (Fig. 7a) is an important driver of nutrient stoichiometry of SGD, but that the annual variation of the Si:P compositional ratio stems from a distinction of biogeochemical histories of SGD happening during the earlier and later parts of the year.

Discussion
Cause-effect relationships that link changes to benthic primary production with nutrient composition of SGD and subterranean estuary dynamics are not well developed. To effectively discuss our results, we must first answer two questions. The foremost is whether our understanding and description of the spatiotemporal dynamics of the STE at the Ria Formosa site is realistic and physically solid; the next is whether scaling up individual measurements of flow and SGD composition taken from seepage meters covering a limited area to monthly quantities obscures the underlying trends.

Describing the oscillation of a subterranean estuary
Groundwater system dynamics in barrier islands are very complex. Consequently, the morphology of freshwater lenses and dynamics of associated transition zones (the subterranean estuary, or STE) are difficult to describe mathematically. The physical problem does not have explicit analytical solutions without considerable simplification (Fetter, 2001(Fetter, , 1972. Our objective is not to provide one solution, but instead to offer a mechanistic framework to understand the link between subterranean estuary dynamics, the composition of STE outflows at the seepage face and their impact on nutrient ratios and benthic primary production. Our SGD measurements clearly showed that there was saltwater as well as freshwater outflow and a well-developed seepage face, implying the saline groundwater body is not static and ruling out the Ghyben and Herzberg principle (Ghyben, 1888;Herzberg, 1901) to define the shape of freshwater lens (Hubbert, 1940). To describe oscillations of the STE at the annual scale, we first accept that these occur around a point of equilibrium (a fulcrum for the oscillation), in response to seasonal changes in recharge throughout the year. We then accept that the oscillation of the STE throughout the year can be reasonably described by taking the monthly steady state solution to the problem of the position of the saltwater wedge toe based on the water table height in the main aquifer. In support of this approach, we first show that the water table in the Campina de Faro aquifer reacts quickly (<1 month) to precipitation (Fig. 2), that the annual fluctuation of piezometric levels follows the effective precipitation over the system (Fig. 3a) and these fluctuations are cyclic on an annual basis over multi-annual timescales (Fig. 3b). This allows us to accept that the transition zone between fresh and salt groundwater also fluctuates on a similar time scale, and hence take the steady-state monthly position of the saltwater wedge calculated as function of the hydraulic gradient on land as an indicator of the location of the subterranean estuary throughout the year. Detailed studies yielding the position of the freshwater-saltwater interface for the case of flow being essentially horizontal (Dupuit assumptions) have been carried out amongst others by Glover (1959), Bear and Dagan (1964), Strack (1972) and Van Der Veer (1977), with the latter accommodating the simultaneous flow of the two fluidshowever, its explicit solution requires an unknown quantity, dependant on the discharge at the shoreline and on the accumulated recharge at the seepage face. Comparisons between all suggested solutions yield very little difference on the position of the tip of the saltwater wedge itself however (Vacher, 1988a). We therefore adopt the formulation of Glover (1959) following the assessment of Cheng and Oazar, (1999). (a, c, e): spring term (January-July) and (b, d, f): fall term (August-December) as discussed in the main text. Statistical information for the six fitted polynomial curves (Eqs. (7)-(10)): R1: R 2 = 0.93, P < 0.001, R2: R 2 = 0.92, P < 0.001, R3: R 2 = 0.85, P < 0.001, R4: R 2 = 0.91, P < 0.001, R5: R 2 = 0.82, P < 0.001, R6: R 2 = 0.90, P < 0.001.
We extend the precaution in our approach by expressing the position of the saltwater wedge as a relative quantity, referenced to the 'equilibrium' position as calculated by its location derived from the mean annual piezometric head, and compare results for four different piezometric series, which yield very similar results (Fig. 3c). Plotting the mean relative position of the saltwater wedge toe against the potential recharge (the excess precipitation) in normalized form (Fig. 3d) yields an annual cycle when 20-year mean excess precipitation records are used, and the notion that relative positioning depends essentially on recharge (see deviation in December 2006, because of the anomalously low precipitation) with the STE extending inland most of the year. We take these results to mean that our approach, though necessarily approximative, is realistic for the purpose of discussing our data.

Scaling of seepage measurements to monthly SGD fluxes
The large spatial and temporal heterogeneity in discharge magnitude and composition of SGD across the beach raises the question of whether scaling up results from the seepage meter array provides representative data for the month without loss or distortion of information. This is an important question -this on-site variability is caused by the combination of tidal and inland pressure gradients acting in opposition (Li et al., 1999), making both the rate of discharge and the relative proportions of fSGD to sSGD vary in the outflowing mixture with beach slope and distance from the tide line during a tidal period (Li et al., 2008). This ratio is also spatially variable in parallel to the coastline due to density effects and local variations in matrix permeability. Therefore, seepage salinity and nutrient stoichiometry become spatially and temporally heterogeneous at the surface length scale of the beach. As shown by Welti et al. (2015), this variability may induce the patchiness of microphytobenthos observed in the Ria Formosa site by Brito et al. (2009), as indeed it regulates subsurface biogeochemical zonation (Ibánhez and Rocha, 2016;Waska et al., 2019). Three separate pieces of evidence ensure our monthly data is realistic and representative of on-site SGD magnitude and composition: Firstly, Michael et al. (2003) studied the differences between deployment of small and large numbers of seepage meters at the 50 m 2 scale and found that arranging them in transects drastically reduced uncertainty in scaling up of fluxes. Ten seepage meters (i.e., 1 per 5 m 2 ) arranged in transects provided the most accurate results. This is the coverage we have for our seepage meter array, and this ensures that the most accurate results possible were obtained at the beach length scale to account for spatial heterogeneity. Temporal heterogeneity is accounted for by the long-term deployment of the seepage meters and the number of samples taken at sub-hourly intervals for never less than two full tidal cycles. The discharge rates are also consistent over time with data obtained from piezometers installed within the beachdata from both methods are compared and discussed in Rocha et al., 2009; Secondly, our seepage meter array produces SGD measurements per unit length of coastline that are consistent with whole-basin groundwater tracer ( 222 Rn) budgets with source attribution by stable isotopes in water, at both seasonal and tidal scales, in 2009 and 2010 . The results for total Rn-derived SGD per metre of coastline into the Ria Formosa obtained then (9.6 ± 9.1 m 3 m −1 d −1 ) are directly comparable to the monthly SGD measurements we obtain from seepage meters in 2006, which vary between 0.67 ± 0.1 and 3.54 ± 1.1 m 3 m −1 d −1 (see Fig. 6b).
Thirdly, as mentioned before, the average fresh groundwater discharge from the M12 unconfined aquifer to the lagoon, per day and per linear metre of shoreline, estimated (see Section 3.4) during 2006 from monthly piezometric data was 0.126 ± 0.017 m 3 m −1 d −1 , statistically identical to the 0.115 ± 0.079 m 3 m −1 d −1 obtained by direct measurement by our seepage meter array.
In addition, we can also assess whether information on nutrient structure of SGD is obscured when scaling individual concentrations (Fig. 7) up to monthly transport rates (Fig. 8a-c). We compare the information given by the annual change in nutrient ratios in monthly fluxes (Fig. 8d-f) with that provided by a frequency analysis of relative limitation/excess of N, P and Si relative Redfield ratios throughout the year as a function of freshwater content (Fig. 10) in individual samples. We find (Fig. 10a, b) that i) the frequency with which N becomes limited with respect to P is drastically reduced with increasing freshwater content of sampleswhich is essentially the message provided by the integrated fluxes (see Fig. 8d for the N:P monthly flux ratio, and Fig. 6b for the fraction of fresh groundwater in total SGD throughout the year); ii) the differences between the spring and fall periods found in nutrient concentration dependence on salinity and particularly in DIN (Fig. 7) is also apparent when looking at potential P limitation (N (or Si):P > 15) and Si limitation (Si:N < 1), at both the individual and the annual scale -the availability of P relative to the other nutrients increases as SGD becomes more saline (Fig. 10c, d), which corresponds to the scaled up monthly flux annual pattern ( Fig. 8d and f); finally, iii), Si:P ratios change from spring to fall particularly when SGD contains more freshwater (>~18 %), which is coherent with the finding that Si:P ratios in SGD (Fig. 8f) depend on freshwater fraction composition (Fig. 6b), which varies during the year as a function of excess monthly precipitation, as illustrated by the comparison between 2006 and the historical means (Fig. 3a). In combination, the data supports the hypothesis whereby nutrient resource ratio availability in SGD is controlled by recharge, and therefore by the annual dynamics of the transitional zone, or subterranean estuary. This conclusion, arising from the analysis of all our data at two different scales, also reinforces the coherence of our scaling up method and adds confidence to ensuing interpretations.

Groundwater discharge composition and dynamics
Our results show that the highest rates of saline groundwater discharge into the Ria Formosa lagoon occur during the driest months of the year. Fresh groundwater discharge magnitudes are higher during the wetter months of the year (January-March, October-November), but freshwater constitutes at most 20 % of the total amount of SGD seeping out into the lagoon (Fig. 6). In the broad sense, the fresh groundwater content of SGD co-varies with the annual continental recharge cycle and the total SGD flux follows in counter cycle to it. The same pattern has been observed by Kelly and Moran (2002), who found SGD into the Petaquamscutt estuary (Rhode Island, USA) reached seasonal highs in summer and lows in winter, by Moore et al. (2006) in the Okatee estuary (South Carolina, USA), who found SGD rates and associated nutrient fluxes in summer higher by a factor of 3-4 than those in winter, and by Charette (2007), who describes an increase of saline SGD during dry periods in the Pamet River estuary (Massachussets). Because SGD is regulated by the hydraulic gradient at the land-ocean interface, this implies that the piezometric head of saline groundwater is positive relative to mean sea level during the drier periods.
Our site is located on the lagoon side of the barrier island strip complex, but during months of positive net recharge (Fig. 3a), the freshwater component of SGD is isotopically linked to the continental aquifer system , which implies that the freshwater lens on the island is connected to the continental aquifer system, either through the Holocene permeable continuum underlying the marsh facies within the lagoon or via a connection to the underlying Miocene unit (Engelen and van Beers, 1986). As shown theoretically by Vacher (1988a) and discussed with case studies elsewhere by Vacher (1988b), Urish and Ozbilgin (1989), Nielsen (1999) and Masterson et al. (2014), the higher saltwater head and greater upper beach storage on the side facing the ocean will cause an asymmetry in the shape of the freshwater lens within the barrier island. The annual dynamics of freshwater lens morphology in these systems is largely controlled by recharge (e.g., Comte et al., 2010). Consequently, two effects of the annual recharge dynamics (Fig. 3) over the Ria Formosa on the morphology of the subterranean estuary (the transition zone) and associated freshwater lens are expected. At the barrier island scale (e.g., Gulley et al., 2016), the thickness (vertical extension) of the transition zone, or the STE, is expected to widen drastically during drought periods, roughly in a 2:1 proportion to the relative reduction in effective precipitation (White and Falkland, 2010). Periods of evapoconcentration (net negative potential recharge, e.g., April/ March to September/October in the Ria Formosa, viz. Fig. 3a) will ultimately lead to the disappearance of the 'pure' freshwater lens and expansion of the transition zone (the STE) into space previously occupied by freshwater. This will imply that the hydraulic gradient established across the barrier island from the ocean to the lagoon will increase salty SGD to the lagoon. It is interesting to note in this regard that the potential recharge, in negative figures (Fig. 3a), which can be seen as moisture demand from soil, increases from~− 10 to~− 145 mm from April to July, with 45 % of this excursion taking place from June to July. These changes would imply an expansion of the STE by a factor of~2 between April and June (White and Falkland, 2010) and roughly the same again between June and July, i.e., in a single month. The hydraulic gradient would hence increase during the whole period by a factor of~3-4. This increase would result in SGD intensifying by the same proportion between April and July, which is in fact what we observe: SGD rises between 1.15 ± 0.12 and 3.54 ± 1.1 m 3 m −1 d −1 , i.e., by a factor of 3 ± 1, between April and July. Apart from the dampening of short-term variability on the data that comes necessarily form the monthly analysis we conduct, that this uptick in SGD is more evident between June and July might be explained by additional factors that increase demand on the groundwater locally, including tourism. Visitors to the Faro region are~4.6 million a year,~45 % of which stay between late June and early September, so the peak season starts effectively in July, placing immediate demand on shallow groundwater in the barrier islands to satisfy demand for water to employ in surface cleaning, dishwashing as well as waste disposal by local hospitality industry. July is also the driest month of the year on average (57 % mean humidity) and the month with the longer daily period of sunshine (average 12 h daily) and could see the most extreme high temperatures during the day (up to 45°C). All these factors contribute to an upconing of the salty groundwater lens and hence in SGD observed. In addition, at the local beach scale, hypersaline beach porewater will sink through the upper beach (Geng et al., 2016) accelerating forced convective circulation (Stringer et al., 2010) and thus enlarging the upper saline plume, or USP (Li et al., 2008;Robinson et al., 2007). The combination of these effects of annual recharge dynamics provides a mechanistic explanation for our observations that may be conceptually illustrated as in Fig. 11.
Briefly, during the 'wet' season, the saltwater wedge toe moves seaward, the subterranean estuary contracts, net recharge contributes to the formation and expansion of a freshwater lens under the barrier island, the positive ocean-lagoon hydraulic gradient promotes freshwater SGD into the lagoon, porewater becomes less saline at the higher beach profile and the USP contracts with the expansion of the freshwater channel (Fig. 11a). Conversely, during the 'dry' season (March/April to September/October), the saltwater wedge toe moves inland, the soil moisture deficit resulting from evaporation leads to the contraction and disappearance of the freshwater lens under the barrier island, the subterranean estuary expands, the positive ocean-lagoon hydraulic gradient drives an increase in salty SGD into the lagoon, porewater becomes more saline at the higher beach profile enhancing density driven circulation in the beach aquifer, and the USP expands with the retraction of freshwater (Fig. 11b). The biogeochemistry of the STE responds as expected (Magaritz and Luzier, 1985) to the perturbations caused by the annual recharge cycle. The pH of STE outflows increases during the dry season and decreases during the wet season, given the major influence of seawater on the dissolved inorganic carbon and alkalinity balance of coastal aquifers (Cai et al., 2003;Mercado, 1985). This behaviour imprints a seasonal variability of STE biogeochemistry and SGD composition that maps groundwater recharge dynamics (i.e., the 'hydrological year'), rather than insolation and temperature (i.e., the Julian year). This realization in turn implies that the ecology of similar coastal systems, connected to nearby aquifers, might be more responsive to groundwater level variance as regulated by both natural and anthropogenic drivers than previously thought. Fig. 11. Conceptual diagram of the oscillation of the freshwater lens and transition zone (the STE) under the barrier islands of the Ria Formosa lagoon in response to effective precipitation. A: 'wet' season. Saltwater wedge toe moves seaward in response to (positive) recharge, the subterranean estuary contracts, freshwater lens expands under the barrier island, the positive ocean -lagoon hydraulic gradient enhances freshwater SGD into the lagoon, porewater becomes less saline at the higher beach profile and the USP contracts with the expansion of the freshwater channel. B: 'dry' season. The saltwater wedge toe moves inland, net evaporation leads to the contraction and eventual disappearance of the freshwater lens under the barrier island, the subterranean estuary expands, the positive ocean-lagoon hydraulic gradient drives an increase in salty SGD into the lagoon, porewater becomes more saline at the higher beach profile, and the USP expands with the retraction of freshwater. The changes in composition of the STE outflows in the form of SGD are highlighted for the 'wet' season (A) and the 'dry' season (B) by boxes with arrowsarrows pointing upward describe an increase in fluxes.

Coupling between STE oscillation and SGD-driven nutrient fluxes
On an annual timescale, fluxes of nitrate and phosphorus to the Ria Formosa lagoon via SGD are highest in summer and fall (July-December) and generally much lower in winter and spring (January-June). The magnitude of fluxes and their composition (Fig. 8a-c) can be explained by the annual transgressive oscillation of the subterranean estuary around its point of equilibrium (Fig. 3c, d) and its drivers (Figs. 2, 3a) rather than by production-respiration processes, which would impact the oxidationreduction potential of SGD, or by the magnitude of the freshwater component of SGD, but is not supported by our observations (Figs. 5b, 6). As the saltwater wedge moves inland (Figs. 3c, 5a) during periods of net negative recharge (Fig. 3a), the pH of SGD increases (Fig. 4c), and SGD-driven fluxes of N and P increase dramatically following the maximum pH attained in July (Fig. 8b, c). This observation is coherent with our conceptual understanding of the dynamics in response to recharge of freshwater and brackish lenses under barrier islands (Fig. 11). Because the expansion of the transition zone beneath the island during periods of potential negative recharge leads to the salinization of aquifer solids and the disappearance of a pure freshwater lens follows periods of soil moisture deficit, the pH of groundwater will increase on par with salinity given the enrichment of saltwater in bicarbonates (Appelo, 1994), while saline SGD will increase driven by the saline groundwater table gradient. This has clear implications for SGDdriven fluxes of P and N to surface waters, and, eventually, also DRSi, since the dissolution of biogenic silica in these environments is expected to be regulated by salinity and pH (Loucaides et al., 2008).
The mechanics of P flux to coastal waters is usually described in association with the dissolution of Fe-oxides, and eventually conversion of these to Fe sulphides in anoxic environments (Kemp et al., 2005;Slomp and Van Cappellen, 2004). However, in oxic environments like ours, benthic nutrient release is less influenced by redox state (Fig. 5b), which has been shown to be controlled essentially by O 2 availability. Spiteri et al. (2006) also showed that the pH gradient (and not the E H ), regulated the speciation of iron oxides and controlled phosphorus sorption dynamics in the STE of Waquoit Bay (Massachussets, USA). Once pH reaches peak annual values (>8.5) in July-September (Figs. 4c, 5a), P release from the Ria Formosa STE is clearly enhanced (Fig. 8c). The elevation of porewater pH, especially when high pH is attained (>8) enhances the mobility of phosphorus by breaking surface Fe\ \P bonds (Gao et al., 2012). This effect has been shown to support P demand in lakes (Xie et al., 2003) as well as tidal regions (Andersen, 1975;Seitzinger, 1991), and is also here shown to be important in supporting benthic primary production (Fig. 9c).
Because exchangeable NH 4 + is weakly bound to negatively charged surface adsorption sites in natural sediments, an increase in the alkalinity and ionic strength of STE fluids with salinization of the coastal aquifer will also necessarily drive changes to the mobility and speciation of nitrogen, since cation exchange reactions can promote NH 4 + mobilization (Rosenfeld, 1979) and the ratio between NH 3 and NH 4 + increases by~one order of magnitude for each pH unit increase (Emerson et al., 1975). Hence elevated pH further accelerates the desorption of exchangeable NH 4 + pools by promoting the formation of NH 3 , which reduces NH 4 + concentrations in both porewater and mineral surfaces and further enhances desorption. Under oxic conditions in situ, an increased availability and mobility of NH 3 , which is the primary substrate for bacterial nitrification (Suzuki et al., 1974), combined with increased salinity and pH driven by salinization, would stimulate nitrification (Isnansetyo et al., 2014;Jones and Hood, 1980;Kemp and Dodds, 2002). Potential enhancement of nitrification in the coastal aquifer driven by the changing chemistry brought by salinization is consistent with the increase in NO 3 − fluxes from the STE verified from July onwards (Fig. 8b) and explains the apparent bi-modal dependence of DIN concentrations on salinity of SGD throughout the year (Fig. 7a, d). The results therefore support the hypothesis whereby an expansion inland of the STE brought by net negative recharge will result in an increase of SGD-borne transport of P and N into the lagoon, because of the stimulating effect of pH and salinity increases on nitrification and phosphorus desorption.
Somewhat in contrast, dissolved reactive silicate fluxes (DRSi) are significantly correlated with the fraction of freshwater in SGD. This merits more research but indicates that an important part of the reactive silicate reaching the Ria Formosa lagoon has a continental, rather than marine, origin. This is not a new idea: Conley (2002) proposed that the fixation of silica in the form of phytoliths by land plants could be comparable to that wrought by marine diatoms. Biogenic silica produced by land plants and freshwater diatoms and mobilized by interactions with seawater was shown to be potentially underrepresented in global silicate balances (Loucaides et al., 2012(Loucaides et al., , 2008, and this topic has received increased attention by the SGD research community (Oehler et al., 2019) following the seminal paper of Kim et al. (2005), which brought to light the potential role of saline SGD (i.e., seawater recirculated through coastal aquifers) in the global Si biogeochemical cycle.

SGD-driven nutrient fluxes and benthic primary production
The evidence here supports the hypothesis whereby the seasonality of SGD in terms of magnitude and composition (Fig. 8) controls the annual dynamics of primary production in the lagoon (Fig. 9). Microphytobenthos production would be the overwhelmingly dominant form of primary production in the system, assimilating most of the sediment-water nutrient fluxes (Fig. 9) by way of its prime location at the sediment-water interface and indeed MPB was estimated to represent up to 99 % of the total Chl a in the system by (Brito et al. 2011). Our hypothesis is also consistent with the absence of correlations between water column Chl a levels and nutrient concentrations from 1991 to 2010 (Brito et al., 2012), between MPB Chl a and tidal range, wind speed, water temperature, salinity or nutrient concentrations (Brito et al., 2009), and, moreover, would be fully consistent with the cyclic components explaining the periodicity of MPB in the lagoon derived from Fourier analysis by Brito et al. (2009). This study shows that 91 % of the variance of MPB in the Ria Formosa can be explained by three components: within-day variance, corresponding to spatial patchiness (61 %), variance with a period of 14-18 days (25 %) and finally a cyclic component with a periodicity of 3 yr −1 . These periodicities are entirely coherent with the annual variability of SGD magnitude and composition and its drivers. Patchiness of benthic Chl a is an inherent result of the within-day temporal variance of SGD cross and along-shore as explained before. The periodicity of 14-18 days is coherent with the oscillation of tidal pressure gradient acting on the beach aquifer over a neap-spring cycle. Acting in opposition to inland hydraulic head, these induce regular fluctuations on the magnitude and composition of SGD at the sediment-water interface, affecting the distribution and production of MPB on a fortnightly timescale. The 3 yr −1 cyclicity is clearly explained by the seasonality of groundwater recharge (Fig. 3a), and its effects on the oscillation of the STE (Fig. 3d), that drives pH in groundwater (Figs. 4,5a) and SGD composition in terms of its salinity (Fig. 6). The interaction between these drivers actually divides the Julian year into three distinctive periods: winter-spring (i.e., January-March), and fall-winter (i.e., September-November) when both fSGD rates and the fSGD/rSGD ratio are relatively high, intercalated by the spring-summer months (April to August) when these were at their annual lows (Fig. 6), with considerable effects on the variability of nutrient structure in SGD (Fig. 8), which evidently conditions MPB (Fig. 9). Given that it is the recharge cycle on land that drives the oscillation of the STE, the production cycles within the lagoon should be analysed against the backdrop of the hydrological year (i.e., 1st October to 30th September) and not the Julian year (1st January to 31st December) as usually done by ecologists. This phase difference between the perceived (seasonal temperature and insolation) and the actual drivers of primary production (effective precipitation) in the Ria Formosa might have contributed to the weak understanding to date of the main drivers of productivity in a lagoon that shows some uncommon traits, including nitrogen limitation of phytoplankton growth in the water column during summer, in spite of luxury consumption of silicate by the local species assemblage (Domingues et al., 2015), which is consistent with the N:P and Si:P ratios of <16 and < 15 respectively found in SGD reaching the lagoon this month (Fig. 8d, f). Groundwater recharge dynamics, leading to a wet season with the STE moving seaward over the period October/November to March/April and a dry season with expansion of the STE inland from March/April to October/November (Fig. 11) may therefore fully explain locally expressed symptoms of eutrophication that have to date confused assessment of the water quality status in the lagoon based on the traditional insolation/production relationships (Newton et al., 2003;Barbosa, 2010). These include occasional occurrence of nuisance macroalgae blooms, mainly in winter (Sprung et al., 2001), low physiological status of marsh plants (Padinha et al., 2000), localized decreases in oxygen saturation levels (Mudge et al., 2008(Mudge et al., , 2007 and episodic water column nutrient enrichment (Cabaço et al., 2008;Newton and Mudge, 2005).

The biogeochemistry and nutrient stoichiometry of SGD in the Ria
Formosa lagoon is controlled by the continental hydrological cycle, or, more precisely, by the oscillation of the groundwater table in the coastal aquifer and its effects on the dynamics of the STE. This makes the dynamics of primary production in the lagoon, and consequently its ecology, dependant on the hydrological year and not on insolation. This is a cause-effect relationship that will prevail in other systems subject to the same drivers elsewhere in the world, and will imply a significantly more direct and short-term impact of climate and environmental changes on ecosystem function as precipitation and groundwater recharge patterns are disrupted by ongoing warming, as well as a more direct dependence of the stability of coastal lagoon ecosystems on groundwater mining, which tends to scale linearly with warming and rising demand on coastal ecosystem services. 2. The dynamics of the groundwater freshwater lens and transition zone is mapped onto the biogeochemistry of STE outflows. The main physicochemical driver is pHwhich increases in SGD with salinization and expansion of the STE into aquifer areas that were previously occupied by the freshwater lens during 'dry' periods and decreases as the saltwater wedge retreats in response to increased groundwater recharge during 'wet' periods. This oscillation of the STE is likely to occur in similar systems throughout the world. The broader implications include stronger links between coastal ecosystem function and precipitation dynamics on the one hand and groundwater resource management on the other, making the issue of true holistic management of the coastal zone considering the subterranean connection in the land-water continuum ) more pressing. 3. Nitrogen and phosphorus respond differently to silicate when the STE oscillates, with Si availability linked more tightly to the freshness of SGD and the saltwater wedge retreat stage while N and P surge primarily as a function of saline intrusion that promotes nitrification and P desorption. This suggests a mechanism by which the biogenic silicates originating in land plants can effectively contribute to the oceanic silicon cycle via the link provided by oscillating STEs. The combined effect of annual STE dynamics and coastal aquifer biogeochemistry results in large increases of the N:Si ratio of SGD during the saline intrusion phase, while the drier months of the year could see N limitation of primary production occurring.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.