Reactive Solute Transport Through Two Contrasting Subterranean Estuary Exit Sites in the Ría de Vigo (NW Iberian Peninsula)

Subterranean estuaries (STEs), where continental groundwaters and saltwaters meet, are zones of intense biogeochemical reactivity. As such, STEs significantly modify groundwater-borne nutrient fluxes to the coastal zone. Thus, evaluating their reactive role is crucial to anticipate impacts of submarine groundwater discharge (SGD) over coastal ecosystems. Here, we studied the nitrogen biogeochemistry of two STEs with contrasting wave-exposure and redox conditions in Panxón and Ladeira beaches (Ría de Vigo, NW Iberian Peninsula). Seasonal surveys were performed at the permanently saturated zone of both beaches during low tide in February, May, July, and October 2019. Sediment was sampled and porewater samples collected using push-pull piezometers. Salinity, 222Rn and 226Ra activities were used to trace water circulation inside each beach. Porewater nitrate, ammonium, nitrite and dissolved oxygen were used to evaluate the role of these STEs as reactive sinks or sources of inorganic nitrogen. Our results showed a marked seasonal variability of water circulation inside both beaches, with strong salinity gradients in February and May and weakened circulation in July and October. The presence of a gravel layer in Panxón beach completely altered the typical structure of STEs by increasing porewater transport and mixing through the beach interior. As a result, Panxón beach profiles were highly enriched in nitrate and oxygen. Conversely, suboxic, and anoxic conditions were prevalent in Ladeira beach during the study period, with ammonium being the prevailing inorganic nitrogen form. High nitrate concentrations occurred associated to the tidal circulation cell during February and May, being the only effective mechanism of sediment oxygenation in Ladeira beach. Although nitrate reduction and production were observed in both STEs, comparison with averaged conservative mixing porewater profiles showed that Ladeira beach acted as a net nitrogen sink whereas Panxón beach acted as a net nitrogen source. The presence of a gravel layer oxygenates the interior of Panxón beach, thus limiting nitrate reduction and promoting the amplification of groundwater-borne nitrogen fluxes to the coast.


INTRODUCTION
Submarine groundwater discharge (SGD) is defined as any and all water flow on continental margins from the seabed to the coastal ocean, including both the discharge of meteoric groundwater from coastal aquifers (fresh SGD) and the recirculation of seawater through the sea floor (Burnett et al., 2003). The importance of fresh groundwater discharge lies in its usual enrichment in nutrients, trace metals, nitrous oxide, methane, or CO 2 as compared with surface waters (e.g., Leote et al., 2008;Liu et al., 2014;Rodellas et al., 2014;Yang et al., 2015). This general enrichment of biogeochemically active solutes is originated in the degradation of organic matter infiltrated into the aquifer, the percolation of nutrients or the contact with metalcontaining aquifer substrates (Slomp and Van Cappellen, 2004;Knee and Paytan, 2012) during the prolonged residence times of groundwaters compared with surface waters. Furthermore, the indiscriminate use of fertilizers in intensive agriculture from the beginning of the twentieth century have provoked an enrichment of nitrogen in aquifers verified worldwide (Valiela et al., 1990;Vitousek et al., 1997;Rivett et al., 2008). This anthropogenic groundwater-borne nitrogen can reach the ocean via different transport processes, including SGD (Beusen et al., 2013), where it has been linked with coastal eutrophication processes (Capone and Slater, 1990;Valiela et al., 1990;Hwang et al., 2005) or harmful algal blooms (Hu et al., 2006;Lee et al., 2010).
Solute flux estimations associated to SGD need to take into account the high biogeochemical reactivity commonly found at the exit sites that can significantly modify the composition of the circulating waters prior to discharge (Slomp and Van Cappellen, 2004;Spiteri et al., 2008;Santos et al., 2008Santos et al., , 2009aRocha et al., 2009;Ibánhez et al., 2011Ibánhez et al., , 2013Ibánhez and Rocha, 2014Jiang et al., 2018Jiang et al., , 2020. This zone is called subterranean estuary (STE), by analogy with surface estuaries (Moore, 1999;Burnett et al., 2003).
In most tidally driven STEs, three water bodies mix as a result of the interaction between coastal and terrestrial forcings (Robinson et al., 2006;Santos et al., 2012). The infiltration of seawater at the upper intertidal and its exfiltration generally at the low tide mark driven by tidal and wave forcing develops a high salinity zone called tidal circulation cell or upper saline plume (Robinson et al., 2006(Robinson et al., , 2007b. Secondary cells could also develop backshore due to berm overtopping and the inundation of backshore depressions (Heiss and Michael, 2014). The intrusion of seawater into the seafloor forms a second high salinity zone seaward called salt-wedge, driven by density differences between fresh groundwater from the coastal aquifer and saltwater (Robinson et al., 2006(Robinson et al., , 2018. Between these two saline zones, the fresh SGD from the coastal aquifer forms a freshwater discharge "tube" (FDT; Robinson et al., 2006). Mixing, flux, and the location of the FDT exit point depend on the hydraulic gradient of the coastal aquifer but also on tidal pumping, wave set-up, density and sediment permeability (Burnett et al., 2003;Santos et al., 2012;Heiss and Michael, 2014;Robinson et al., 2018). The expansion of the tidal circulation cell can make the FDT shallower or non-existent as it causes vigorous saltwaterfreshwater mixing due to unstable density distributions inside the beach leading to convective circulation (Robinson et al., 2007a). Spatial sediment heterogeneity within the intertidal zone can alter porewater circulation, creating zones of preferential flux, and complex local-scale mixing rates (Geng et al., 2020a,b).
Biogeochemical reactivity within STEs has been linked to different porewater circulation patterns in numerous studies (Spiteri et al., 2008;Santos et al., 2009a;Ibánhez et al., 2013;Heiss and Michael, 2014;Rocha, 2016, 2017;Couturier et al., 2017). During high tide, seawater infiltration resupplies the sediment with fresh organic matter, and dissolved oxygen (Ibánhez and Rocha, 2016). Particulate organic matter is retained in the shallower depths, whereas dissolved organic matter penetrates the beach interior driven by the recirculation of seawater though the sediment (Rusch et al., 2000;Santos et al., 2009a;Huettel et al., 2014). The degradation of organic matter can lead to the attenuation of terrestrial groundwaterborne N through nitrate reduction processes within the STE, such as denitrification (Kroeger and Charette, 2008;Thamdrup, 2012). Autotrophic processes, such as anaerobic ammonium oxidation (ANAMMOX), are also important sinks of N and they have been reported in STEs (Kroeger and Charette, 2008;Sáenz et al., 2012;Thamdrup, 2012). In anoxic conditions, the dissimilatory nitrate reduction to ammonium process competes with denitrification for both electron donors and acceptors, although the reduction of nitrate to ammonium seems to be favored at strong reducing conditions (Thamdrup and Dalsgaard, 2008). On the contrary, in oxic STEs, limited N consumption together with the production of ammonium from organic matter mineralization can enhance N loads to surface water bodies (Ibánhez et al., 2013;Ibánhez and Rocha, 2017). The mixing between seawater and meteoric groundwater, together with enhanced transport caused by the permeable nature of STEs, creates a strong biogeochemical zonation that enables the co-occurrence of all these processes in close proximity. However, the modulation role of STEs over landderived N fluxes to the coast and its controlling factors are still highly uncertain.
For the first time in the Ría de Vigo (NW Iberian Peninsula), two contrasting STEs were examined (Figure 1). The evaluation of the mechanisms that control the final SGD N loads to the receiving waters is particularly important in this embayment, as its valuable natural resources depend on the good environmental status of these coastal waters. Both STEs are placed next to each other and yet, they present differential redox conditions. Environmental factors underlying these differences are explored and related to the cycling of the redox-sensitive inorganic N compounds within these STEs. The final objective of this study is to evaluate the factors determining the capacity of STEs to mitigate or amplify groundwater-borne N loads to the receiving waters.

Study Site
The Ría de Vigo (NW Iberian margin; Figure 1) is a large coastal inlet at the northern limit of the Eastern Boundary Upwelling Ecosystem of the Iberian-Canary current (Arístegui et al., 2009). FIGURE 1 | Study area. Sampling sites are identified by red dots. The location of the wells and boreholes used in this study are represented with green triangles, while the location of the Miñor river sampling site is shown with a yellow triangle. These data are obtained from Ibánhez et al. (2019). Orthophoto (taken in 2018) and the mapping of local fluvial courses were obtained from the Instituto Geográfico Nacional (www.ign.es). BH denotes borehole.
Its morphology resembles a funnel oriented SE-NW and protected at its mouth by the Cíes Islands. Depending on the dominance of oceanic or riverine inputs, the Ría can be divided in two zones: an inner estuarine zone, formed by San Simón Bay and Baiona Bay, and an external oceanic zone, which covers the rest of the Ría. Intermittent upwelling of Eastern North Atlantic Central Water during March-April to September-October fertilizes the Ría (Álvarez-Salgado et al., 2003), sustaining intensive mussel aquaculture, shellfish harvesting, and small-scale fisheries in the region (Fernández et al., 2016). The coastal inlets in this area supports 15% of the world's mussel production (around 250,000 tons), billing 112 million € in 2015 and supporting more than 8,000 full-time employees (Labarta and Fernández-Reiriz, 2019).
Using radioisotope tracers, Ibánhez et al. (2019) recently identified and quantified SGD in the Ría de Vigo. San Simón and Baiona Bays (Figure 1) were among the most affected zones by this diffusive source of continental water. This study was placed in two wave-dominated beaches located in Baiona Bay: Ladeira and Panxón (Figure 1). Baiona Bay is exposed to westerly and southwesterly winds, allowing the most energetic waves to enter the bay during winter storms. During the summer, the waves are less energetic due to the protection of the bay from northerly winds, making oceanic swell the most important energy source most of the year (Alejo et al., 1999). Tides are semidiurnal, with a range of 3-3.5 m during spring tides and 1.5-2 m during neap tides (Martín Míguez, 2002).
Ladeira beach is located in the inner part of Baiona Bay, protecting the A Ramallosa beach-barrier-lagoon complex. Its sedimentary composition is dominated by siliciclastic sands with carbonate content below 5%, mean grain size of 0.245 mm (Queralt et al., 2002) and little variation in grain size along the beach profile (González-Villanueva et al., 2007). It is influenced by fluvial sediments transported by the Miñor, A Groba and Belesar rivers. These rivers converge at the northern edge of the littoral arrow that connects the A Ramallosa complex and Baiona Bay (Figure 1). Backshore and foreshore zones appear well-defined in this beach, with vegetated dunes and beach cusps degraded by human action. Seaweed is accumulated in the intertidal zone during autumn and spring and prevents wave energy dissipation, creating an erosive vortex that develops steeps up to 70 cm high. This type of morphology could be preserved during winter, but it appears covered by sediment the rest of the year (Alejo and Vilas, 1987).
Panxón beach is located in the outer part of Baiona Bay, bounded by Punta do Castro to the North and by the Muiños River to the South (Figure 1). This beach is more exposed to wave action compared to Ladeira beach. The backshore ridge of dunes has been altered by human action, and the construction of a pier created a shadow zone in terms of sediment accumulation. Sediment is composed by unimodal fine sands with a mediumhigh selection, coarser and with poorer selection in more exposed zones of the beach (Álvarez-Vázquez et al., 2003). Mean grain size in Panxón is 0.236 mm with high content in carbonates (60-65%; Queralt et al., 2002).
Both beaches are surrounded by the Monteferro-O Rosal geological complex, which is a metasedimentary formation dominated by shales with intrusions of igneous rocks such as granite, pegmatite, pegmaplite and quartz. This complex is bounded by granitic igneous rocks from the Hercynian age with alkaline affinity to both the East and the West. These crystalline rocks are traditionally considered to have low permeability (Gustafson and Krásný, 1994), but the high fracturing present in the zone alongside with the high weathering of the surface layer (regolith) give them the capacity to store groundwater. Nevertheless, Holocene alluvial and colluvial sedimentary deposits appear associated to the studied beaches and the Monteferro-O Rosal geological complex (Pliego Dones et al., 1978), which further increases groundwater storing capacity in the region. These two beaches were selected attending to their close proximity and yet contrasting characteristics such as wave exposure, organic matter content, oxygenation status, and sediment composition.

Sampling Strategy
Seasonal surveys in Panxón and Ladeira beaches took place on the 20th and 22nd of February, 6th and 8th of May, 2nd and 4th of July and 29th and 30th of October 2019, respectively. These dates were selected to coincide with spring tides and thus, minimize biogeochemical differences caused by changes in the tidal regime. Three sampling stations were selected in the upper, middle and lower limits of the permanently saturated sediments of the intertidal zone (Figure 2). At each station, push-pull piezometers (M.H.E. Products, United States; 4 cm of screen-length) were buried at 178, 116, 86, 56, 29, and 10 cm depth within the sediment (the 10 cm depth piezometer was not installed in February). The presence of a gravel layer in Panxón (Álvarez-Vázquez et al., 2003) at a variable depth between 90 and 150 m complicated the installation of the deepest piezometers (116 and 178 cm). As a consequence, the installation depth of the longest piezometer was variable there. Beach profiles in February were measured with a GPS and referenced to the Alicante (SE Spain) datum. The remaining beach profiles from May to October were taken with a scaled rod as the slope from a reference point set at the backshore.
During the four surveys performed at each beach, water sample collection began after the installation of the piezometers during low tide. Porewater samples were collected with air-tight syringes after discarding three times the internal volume of the piezometers. Water samples collected from each piezometer were used to determine salinity, dissolved oxygen (O 2 ), radon ( 222 Rn), and Ladeira (lower panel). Panxón average beach slope was 3.00% (February: 2.90%, May: 2.44%, July: 3.31%, and October: 3.38%). Ladeira average beach slope was 6.86% (February: 6.71%, May: 7.10%, July: 6.54%, and October: 7.08%). Beach slopes were calculated with the data from the intertidal zone. Black dots represent the averaged location of the three sampling stations. The dark gray area represents the permanently saturated zone. HT, approximate position of high tide; LT, approximate position of low tide. radium ( 226 Ra), and dissolved inorganic nitrogen (DIN; NO 3 − , NO 2 − , and NH 4 + ). Temperature was measured onsite with a thermometer with an accuracy of ± 0.1 • C.
Additional sediment samples were taken from the permanently saturated zone of the beach with a PVC minicorer (10 cm long, 2.5 cm of diameter) from the lower and upper sampling stations in May, July and October. In February, sediment samples were taken from the middle station. At each point, two minicorers were collected, sliced onsite into three different depths (0-3, 3-6, and 6-9 cm in February and 0-2, 2-6, and 6-10 cm in the rest of the surveys) and equivalent depths were mixed to reduce spatial variability. The homogenized sediment samples were immediately transferred to the lab and stored at −20 • C until further processing. Sediment samples were used to determine water content, organic matter and total nitrogen and organic carbon content.

Analytical Methods
Sample volume collection was limited as much as possible to avoid potential spatial interference between piezometers. Thus, analytical techniques were adapted to low volume water samples. A Guildline Portasal salinometer was used to measure salinity. Nitrate and nitrite were measured following standard colorimetric methods (Grasshoff et al., 1999) on an Alliance Futura segmented flow autoanalyzer. In the case of ammonium the fluorimetric method of Kérouel and Aminot (1997) was used. Dissolved oxygen was measured following the spectrophotometric Winkler method (Labasque et al., 2004) adapted for low sample volumes (25 mL).
Radon activities were measured within 2 days after collection using a Rad7 Radon detector with a RadH2O accessory (Durridge Company, Inc.). The collected samples (70 mL) were transferred into 250 mL flasks and diluted with 222 Rn-free water. Readings were corrected for the internal 222 Rn decay from the time of collection to analysis. To obtain 226 Ra activities, purged samples were stored until 222 Rn and 226 Ra reached secular equilibrium (∼30 days, i.e., more than 5 times the 222 Rn halflife) and measured again. The equipment was purged with an activated carbon 222 Rn trap to minimize background activities. The detection limit of our procedure determined from the regular analysis of blanks (n = 36) made with 222 Rn-free water (250 mL) was 36 Bq m −3 . Thus, the detection limit of 222 Rn and 226 Ra in our diluted samples is estimated at 128 Bq m −3 . Porewater 222 Rn equilibrium activities were measured with sediment incubations (∼200 g dry weight; n = 9 for each beach) run with 222 Rn-free water (Colbert and Hammond, 2008). After 30-40 days of incubation, 222 Rn activities were determined following the same procedure as for field water samples. Results are expressed as the average and the standard error of the 9 equilibrium activities determined.Porosity was approximated to the water content in the sediment as follows: where V pw is the volume of porewater and V sand the volume of sand. Sediment subsamples were dried at 60 • C until constant weight was reached. V pw was then calculated using the salinities determined in the shallower depths of each beach. Local dry sediment density was determined (2.599 ± 0.001 kg/L in Panxón and 2.640 ± 0.038 kg/L in Ladeira) and used to calculate V sand . The percentage of organic matter present in the sediment was determined by loss-on-ignition (Ball, 1964). To obtain total N and organic C content, sediment samples were powdered using a laboratory ball mill with natural agate tank during 10-15 min at medium intensity. Then, samples were exposed to HCl flumes during 6 h to remove inorganic C. Afterward, samples were analyzed in a CHN elemental analyzer Perkin Elmer 2400 A.

Statistical Analyses
Statistical analyses presented in this study were carried out in R 4.0.0 with the Rstudio 1.2.5019 editor. Normality of the measured variables was tested with the Shapiro-Wilk test. Spearman correlation analysis was performed using the PerformanceAnalytics 2.0.4 package for R, distinguishing between beaches and among surveys. The Kruskal-Wallis test was applied to verify whether the differences between groups were statistically significant. Principal Component Analysis (PCA) was performed with the "princomp" function from the stats 4.0.0 package and represented using the ggbiplot 0.55 package in R. Variables were scaled before performing the PCA.

Meteorological and Oceanographic Data
Meteorological and oceanographic conditions in the study area were evaluated for a period of 30 days prior to each survey (winter with data from the 20th January to the 22nd February, spring with the data from the 6th April to the 8th May, summer from the 2nd June to the 4th July and autumn from the 30th September to the 30th October; Supplementary  Table 1). Daily water balance (Supplementary Table 1) was obtained from the repository of Meteogalicia 1 using the data from the meteorological platform of Baiona (autumn; 42 • 6.93 N 8 • 50.23 W) and Oia (remaining seasons; 41 • 59.67 N 8 • 51.97 W). Water balance integrated for the month prior to each survey was calculated as the sum of daily water balances for that period. Water balance is used as a qualitative indicator of groundwater table fluctuations due to the fast aquifer rechargedischarge rates and fast response to environmental conditions verified in the local crystalline aquifers (Raposo et al., 2012). Minimum and maximum wave heights during the studied periods (Supplementary Table 1

Sampling Site Characterization
Autumn showed the most positive water balance, followed by spring and winter. On the contrary, summer was characterized by low precipitation and a resulting negative water balance (Supplementary Table 1). Storm events were recorded in the days prior to the February and May surveys, whereas in July and October calm conditions dominated. Although being the rainiest period, autumn showed the highest averaged salinity within Baiona Bay, followed by summer and winter (Supplementary Table 1). Dissolved oxygen levels in surface waters of the bay remained close to saturation throughout the studied period (Supplementary Table 1).
The beach profile of Panxón is typical of a dissipative beach, with a leveled morphology caused by high-energy conditions (Figure 2). The maximum slope was reached in October (3.38%) and the minimum slope in May (2.44%). Conversely, the beach profiles of Ladeira presented characteristics of a reflective beach, with a well-developed supra-tidal and a marked berm throughout the year (Figure 2). This is in accordance with its lower exposure to the dominant waves compared to Panxón beach. The maximum beach slope was reached during May (7.10%) and the minimum in July (6.54%). Following the beach berm, a 50 cm high slope appeared during October and probably caused by storm waves developed in the days prior to this survey (Supplementary Table 1). The abnormal storm frequency occurring in summer 2019 caused a switch in the optimum season for the maximum sediment accumulation at both beaches, which occurred in autumn rather than during the summer.
Panxón and Ladeira beaches presented similar porosity (average: 0.47 ± 0.02 in Panxón and 0.46 ± 0.06 in Ladeira; Kruskal-Wallis test p-value > 0.05) but significantly different sedimentary organic matter content (average: 1.33 ± 0.10% dry weight in Panxón and 0.93 ± 0.40% dry weight in Ladeira; Kruskal-Wallis test p-value < 0.01) and C:N molar ratio (average: 3.23 ± 0.81 in Panxón and 10.71 ± 3.48 in Ladeira; Kruskal-Wallis test p-value < 0.01). Organic matter content decreased with depth in both beaches (Figure 3), with the exception of the lower station in May and October at Panxón and in February at Ladeira. Organic matter content decline with depth was higher in Ladeira than in Panxón, especially in July (from 1.78 to 0.76% in the lower station and from 1.34 to 0.63% in the upper station;  Seasonal differences in the measured sediment properties were only observed in the organic matter content at Panxón beach (Kruskal-Wallis test p-value < 0.05). The lack of organic matter peaks with depth in the collected 10 cm length sediment corers suggests that these differences may be related to seasonal changes in the supply of pelagic organic matter from surface waters rather than seasonal sediment transport along the beach profile.
Porosity and organic matter content appeared significantly correlated in both Panxón (p < 0.05, r = 0.46, n = 21) and Ladeira (p < 0.001, r = 0.76, n = 21). Moreover, Ladeira showed high correlations between porosity and both organic C content (p < 0.001, r = 0.72, n = 21) and total N content (p < 0.001, r = 0.72, n = 20); and between organic matter and both organic C content and total N content (p < 0.001, r = 0.95, n = 21, and r = 0.95, n = 20, respectively). In Ladeira, total N content and total organic C content also appeared significantly correlated (p < 0.001, r = 0.92, n = 20) although the C:N molar ratio did not show any relationship with the other variables. On the contrary, in Panxón, organic matter content was significantly correlated with total N content (p < 0.05, r = 0.47, n = 21).

Tracing Water Circulation: Salinity, 222 Rn and 226 Ra Profiles
Porewater salinity was lower than that of Baiona Bay (Figure 4 and Supplementary Table 1) revealing the existence of a freshwater influx to the beach interior. Moreover, high porewater salinity gradients were observed, even in the same month and beach. The high porewater salinity values found in October (30.7 ± 0.7-31.4 ± 1.8) in Panxón did not occur in Ladeira, and the maximum average salinity in Ladeira reported at the lower station of May (32.1 ± 1.1) did not appear reflected in the parallel survey performed in Panxón. This suggests that seasonal salinity changes are specific of each beach, probably due to differences in circulation patterns and freshwater supply.
Porewater 222 Rn equilibrium activities determined with sediment incubations were 3698 ± 251 Bq m −3 from Panxón sediments and 1610 ± 329 Bq m −3 from Ladeira sediments. Porewater 222 Rn activities were, in most cases, greater than these 222 Rn equilibrium activities, suggesting external inputs of 222 Rn to the beach porewaters. Exceptions are found in some sampled porewaters near the sediment surface (Figure 4). Average activities of dissolved 226 Ra were one order of magnitude lower than the average activities of 222 Rn and seasonal and spatial differences were also found.
Vertical gradients of salinity and radioisotopes were observed throughout the four sampling periods (Figure 4). Decreases in salinity with depth were recorded in most porewater profiles (e.g., between 29 and 56 cm depth at the middle station of February in Panxón), followed in some of them by new increases (e.g., between 86 and 116 cm depth in the same porewater profile). Temperature variations along the profiles were too small (1-2 • C) to cause an increase in density that would support the stability of the salinity profile. Consequently, unstable density profiles were maintained by the flow of the freshwater and seawater bodies that interact in the STE. Only in a few porewater profiles (at the upper station in May, and at the lower and middle station in July, in Panxón) salinity was mostly vertically homogeneous.
Tracing Nitrogen Cycling: Nitrate, Nitrite, Ammonium, and Dissolved Oxygen Porewater Profiles Panxón beach was more oxygenated, with higher profile-averaged concentrations of dissolved oxygen (from 19.8 ± 22.9 µM to 214.2 ± 21.1 µM) than Ladeira (from 0.0 ± 0.0 µM to 22.9 ± 15.2 µM), where suboxic (≤ 4.5 µM) and anoxic conditions dominated (Figure 5). These differences in oxygenation were strongly linked to the differences found in inorganic nitrogen speciation in porewaters. Nitrate concentrations were between one and two orders of magnitude higher in Panxón compared to Ladeira. On the contrary, ammonium concentrations were higher in Ladeira (Figure 5).
The largest variations in porewater dissolved oxygen concentrations occurred between stations sampled in the same month and beach (Figure 5). Upper stations presented, in most cases, higher levels than middle and lower stations. Moreover, the oxygenation in Panxón beach occurred in depth (up to 240 µM), and lowest oxygen levels were found close to the sediment surface (∼ 10 µM; Figure 5). In Ladeira, increases in dissolved oxygen concentrations appeared both at shallow and at intermediate depths, as shown in the upper porewater profiles of February and May (Figure 5).
Nitrate concentration peaks appeared linked to increases in oxygen levels in most porewater profiles. However, nitrate was also present at low concentrations in some of the shallower depths during July and October in Panxón. Dissolved oxygen was significantly correlated with nitrate in February (p < 0.05, r = 0.69, n = 14), May (p < 0.01; r = 0.66, n = 17) and July (p < 0.01; r = 0.66, n = 18) in Panxón beach and in May (p < 0.05, r = 0.59, n = 18) and July (p < 0.05, r = −0.52, n = 18) in Ladeira beach. On the contrary, ammonium generally appeared at shallow depths negatively correlated with oxygen in Ladeira (p < 0.01, r = −0.51, n = 69) and with nitrate in both Panxón (p < 0.001, r = −0.42, n = 67) and Ladeira (p < 0.001, r = −0.24, n = 69). During February and May high ammonium concentrations (>30 µM) were present throughout the entire vertical porewater profile at the lower station in Ladeira. Nitrite was associated to shallow depths similarly to ammonium, particularly in Panxón beach (Figure 5).

Relationships Between Circulation Tracers and Nutrient Profiles
Spearman correlation analysis showed significant relationships between nutrient and internal porewater circulation tracers. In Panxón beach, nitrate and 222 Rn (p < 0.001; r = 0.47, n = 67) as well as oxygen and salinity (p < 0.05, r = −0.25, n = 67) were significantly correlated. Furthermore, in October, the month with the most positive water balance (Supplementary Table 1), oxygen and 222 Rn (p < 0.05, r = 0.57, n = 18), as well as salinity and nitrate (p < 0.01, r = −0.67, n = 18) showed strong correlations. This adds to the correlations found between 222 Rn and salinity there and suggests a relationship between nitrate enriched, oxygenated porewater and 222 Rnenriched, low salinity porewater. By contrast, in Ladeira ammonium was significantly correlated with salinity (p < 0.05; r = 0.24, n = 69) and showed negative correlation with 222 Rn (p < 0.001; r = −0.53, n = 69). Thus, ammonium seems to be related to infiltrated seawater rather than fresh groundwater in Ladeira beach.
A principal component analysis (PCA) of the measured variables in Panxón and Ladeira is presented combined ( Figure 6A) and individually for each beach (Figures 6B,C). Components 1 (C1) and 2 (C2) account for 31.2 and 17.1% of the total variance in the PCA performed with the data from both beaches. There, component C1 is related to the redox conditions, as it is dominated by ammonium (0.51), nitrate (−0.50) and dissolved oxygen (−0.48); whereas Component C2 is more related to transport processes, as it was dominated by 222 Rn (−0.52) and salinity (0.47), in addition to dissolved oxygen (0.48) and nitrate (0.46; Figure 6A). The data from each beach appeared clustered and well separated from each other in this PCA, giving a different sign of the slope of a linear regression of PCA components 1 and 2 ( Figure 6A). A more detailed analysis is provided by the individual PCAs. In the PCA performed with the data obtained in Panxón (Figure 6B), Components C1 and C2 account for 28.4 and 21.0% of the total variance. Component C1 is dominated by nitrate (−0.59) and oxygen (−0.49), while Component C2 is dominated by salinity (−0.54), 226 Ra (−0.45) and ammonium (−0.42). In Ladeira (Figure 6C), PCA components C1 and C2 account for 30.3 and 27.1% of the total variance. C1 is dominated by oxygen (−0.59) and nitrate (−0.56) and C2 is dominated by 222 Rn (−0.59), salinity (0.52) and ammonium (0.47). Both PCAs confirmed the relationships observed through Spearman correlation analysis.

The Source of 222 Rn-Rich Porewater
Highly fractured granites and shales dominate the geology of the NW Iberian Peninsula (Llerena et al., 2013) and are part of the Monteferro-O Rosal geological complex (Pliego Dones et al., 1978) that surrounds Baiona Bay. These geological substrates are naturally enriched in isotopes from the 228 U decay chain (Tollefsen et al., 2016), such as 226 Ra and 222 Rn. Therefore, groundwaters contained within these rocks are commonly enriched in these radioisotopes, particularly 222 Rn since 226 Ra shows very low solubility at low salinities .
Porewater 222 Rn equilibrium activities estimated from sediment incubations together with the measured dissolved 226 Ra activities (∼130-1,600 Bq m −3 ) are not sufficient to explain the observed 222 Rn enrichment in the sampled porewaters. Furthermore, the negative correlation between 222 Rn and salinity found in some of the surveys and the opposite PCA scores between 222 Rn and salinity in Ladeira (Figure 6C) point toward a freshwater origin for this 222 Rn enrichment. Ibánhez et al. (2019) measured radon activities at various freshwater bodies near Baiona Bay (Table 1). There, 222 Rn activities found in the Miñor River (Figure 1; 2.7 ± 0.1 ×10 3 Bq m −3 in winter/spring and 2.5 ± 0.1 ×10 3 Bq m −3 in summer) are one order of magnitude below the maximum 222 Rn activities observed in Ladeira. On the other hand, 222 Rn activities found in boreholes and wells near the study area (Camos borehole and well, Ramallosa borehole and Panxón well; Figure 1) were very high (between 32 ± 2 ×10 3 Bq m −3 and 802 ± 28 ×10 3 Bq m −3 ; Table 1). Thus, considering that rainfall has low 222 Rn content, the high radon levels found in beach porewaters can only be sourced from groundwater from the coastal aquifers.
Assuming that 222 Rn is in secular equilibrium within the coastal aquifer and that no net 222 Rn loss occurred during the transit to the studied beaches except from 222 Rn decay, the difference between the radon activities found in the proximate groundwaters (Rn aquifer ; Table 1) and those found in the beach porewaters could provide a first estimate of the transit time of groundwater from the terminus of the coastal aquifer to the beach surface. 222 Rn activities in excess (Rn excess ) are defined as the difference between measured porewater 222 Rn activities and those supported by porewater 222 Rn equilibrium activities plus porewater dissolved 226 Ra activities. We further assume that measured porewater 222 Rn equilibrium activities in shallow where λ is the decay constant of 222 Rn (0.18 d −1 ) and Rn extrapolated is the 222 Rn activities in excess observed at each beach and extrapolated to 0 salinity. This extrapolation was performed using the percentage of fresh-groundwater (% fresh ) and Rn excess in the sample: Rn extrapolated = 100 · Rn excess % fresh The percentage of fresh-groundwater (% fresh ) was obtained using as saltwater endmember the sample with the maximum salinity measured at each beach and survey. Only samples with a percentage of fresh-groundwater above the median of each beach and survey were used to guarantee a Rn excess signal sufficiently high, thus limiting our extrapolations (Supplementary Table 2). Transit times between 3.5 ± 2.6 and 22.8 ± 1.6 days were obtained in Panxón beach and between 1.6 ± 1.6 and 21.4 ± 4.4 days in Ladeira beach, depending on which groundwater endmember is used (Supplementary Table 2). These values are in the range of other transit times found elsewhere (Michael et al., 2005;Bokuniewicz et al., 2008). Transit times were higher in summer than winter, probably because higher aquifer residence times result from lower rainfall and higher evaporation, which might have resulted in lower groundwater table during the summer (Supplementary  Tables 1, 2). The shortest apparent transit times were calculated when assuming that the coastal aquifer supports Rn levels as those measured in the wells from Camos and Panxón (between 3.4 ± 1.6 and 8.9 ± 4.4 days when using the Camos well, and between 1.6 ± 1.6 and 6.0 ± 1.6 days when using the Panxón well; Supplementary Table 2). Radon activities measured in boreholes and wells spans one order of magnitude (Table 1) and this strongly affects the calculated transit times.
Transit time calculations performed here assume a steadystate, advection-dominated and homogeneous flow system (Michael et al., 2011). The hydrological Peclet number (P e ) was calculated for each of the transit times to evaluate the assumption of an advection-dominated system (Supplementary Table 4): where L is the characteristic pore length, approximated to the averaged grain size of each beach, D molec is the diffusion coefficient of 222 Rn in water and v the groundwater velocity. A D molec of 1.102 × 10 −5 cm 2 /s obtained in seawater at 17 • C (Dueñas et al., 1985), i.e., the averaged measured temperatures observed in both STEs, was used. Groundwater velocities (Supplementary Table 4) were obtained using the calculated transit times (Supplementary Table 2) and the shortest path between each sampling site and the nearest U-enriched body rock (300 m for Panxón beach and 580 m for Ladeira beach; Pliego Dones et al., 1978). In all cases, P e numbers were higher than 1 (Supplementary Table 4), confirming that advection dominates over diffusion regardless of the groundwater endmember used.
A sensitivity analysis was performed to determine the relative importance of Rn aquifer and Rn extrapolated in the calculated transit times as follows (Ibánhez et al., 2011): where x is the transit time and y the variable tested (Rn aquifer or Rn extrapolated ) and x and y the respective ranges of variation. A 10% of variation was imposed to each variable tested ( y). Results showed that transit times are very sensitive (S > > 1) only at low Rn extrapolated (i.e., at low unsupported 222 Rn enrichment in the beach porewaters) and Rn aquifer values (< 50 × 10 3 Bq m −3 ; i.e., those measured in wells surrounding the study sites). Conceptually, the local aquifers can be approximated to a two-layer aquifer: a surface layer contained within the regolith and the highly fractured surface basement rock where wells are screened, and a deeper layer contained within the fractured basement rock where boreholes are screened (Raposo et al., 2012). Aquifer specific storage and water transmissivities in the study region are orders of magnitude higher in the surface aquifer layer compared to the fractured rock underneath (Raposo et al., 2012 and references therein). These characteristics determine contrasting matrix porosity and groundwater residence times that may explain the differences in 222 Rn activities observed between both aquifer layers. Furthermore, these contrasting hydrological characteristics of both aquifer layers suggest the surface layer as the most probable source of continental groundwater to the studied STEs. Nevertheless, further analysis is needed to accurately identify the groundwater endmembers and transit times of both STEs.

Porewater Circulation in the Ladeira and Panxón STEs
The vertical distribution of Rn excess can be used to identify the location of the radon-enriched FDT in the measured porewater profiles (Figure 7). For this, homogeneous porewater 222 Rn equilibrium activities with depth are assumed. The use of Rn excess to identify the FDT is also limited by internal mixing, that could promote high salinity porewaters associated to high 222 Rn activities.
Radon activity "excess" is observed in most of the porewater profiles of Ladeira and Panxón beaches (Figure 7). A peak in Rn excess was observed between 86 and 116 cm in the upper station of July and in the three stations of October in Panxón (Figure 7). There, the gravel layer detected between 90 and 150 cm likely represents a zone of preferential porewater transport. This modifies the typical porewater circulation structure of STEs, contracting the salt-wedge to shallower depths and increasing mixing, especially in October (Figure 8). Secondary peaks of Rn excess (Figure 7) could be related to density instabilities (Simmons et al., 2001;Robinson et al., 2006) or to zones with higher permeability with respect to adjacent areas (Geng et al., 2020a,b).
The tidal circulation cell and the salt-wedge can be detected in the measured porewater profiles associated to high salinities (>30; Figures 4, 8), but also to 222 Rn activities lower than porewater 222 Rn equilibrium (Figure 8). This is an indication of porewaters with low residence time, i.e., recently exchanged with surface waters. Porewaters with high salinity and radon activities below equilibrium were found in the upper part of the sampled beaches, identified as the tidal circulation cell, and below the FDT, identified as the salt-wedge (Figures 4, 8).
The seasonal variation of the continental hydraulic gradient and the occurrence of wave events are related with the position of the FDT exit point and the size of the tidal circulation cell in STEs (Robinson et al., 2007b(Robinson et al., , 2014(Robinson et al., , 2018Heiss and Michael, 2014;Xin et al., 2015). High inland water table could cause the narrowing of the tidal circulation cell and the displacement of the FDT exit site landwards (Heiss and Michael, 2014). This seems to be the case during February and May in both beaches (Figures 4, 7), which agrees with the positive water balance observed during these months (Supplementary Table 1). Increased wave forcing has been related to strengthen flows in the tidal circulation cell (Robinson et al., 2014;Xin et al., 2015). Apparently, wave events prior to the surveys of February and May (Supplementary Table 1) might have influenced the STE (Figures 4, 8), where high salinity and 222 Rn activity gradients showed a well-developed tidal circulation cell (Figures 4, 8). Furthermore, the tidal circulation cell penetrated deeper in the upper station of Panxón than in the upper station of Ladeira (Figure 4). This suggests that Panxón beach is more sensitive to wave forcing than Ladeira beach, which agrees with its higher wave-exposure. The lower salinities of the FDT in February compared to May might be related to a higher inland hydraulic gradient during February.
July showed vigorous mixing in both beaches, probably caused by a low inland hydraulic gradient and calm wave conditions (Figure 4). However, in Panxón, the reduced salinities of the FDT in July suggest a higher fresh groundwater flux, which contrasts with the negative water balance found during the summer (Supplementary Table 1) and suggests a time lag between aquifer recharge and discharge. October, which showed the highest water balance of the studied period (Supplementary Figure 1), showed a radon-enriched and high salinity FDT in Panxón beach and lower salinity FDT in Ladeira beach (Figures 4, 8). Increased residence times and seawater intrusion into the aquifer due to reduced inland hydraulic gradient during the summer might explain the high salinities and 222 Rn activities found in October at Panxón beach. Adsorbed 226 Ra at low salinities can be mobilized due to ionic exchange as salinity increases (e.g., during seawater intrusion at low aquifer pressure periods; Dulaiova et al., 2008;Gonneea et al., 2008). With the recovery of the continental groundwater table during autumn, this mobilized 226 Ra will then exit the STE as observed during the October survey in Panxón beach. This seems to further support that the high salinity, 222 Rn-enriched waters observed during October in Panxón are originated from seawater previously intruded into the aquifer during the summer.

Oxygen and Organic Matter Transfer and Reactivity in the Panxón and Ladeira STEs
In STEs, the tidal circulation cell is, together with the inflow of oxygenated terrestrial groundwater, the main driver of oxygenation of the beach interior (Robinson et al., 2007a;Ibánhez and Rocha, 2016). This is well illustrated by the measured dissolved oxygen porewater profiles at Ladeira beach in February and May. High porewater oxygen concentrations (>20 µM) coincided with radon activities near equilibrium (Figures 5, 7) and with the high salinity zones at the upper station (Figure 4), where the tidal circulation cell was active. Conditions became hypoxic in July and October, probably due to the weakening of the tidal circulation cell and/or higher benthic oxygen consumption rates. FIGURE 7 | Measured porewater 222 Rn activities (black) and equilibrium 222 Rn activities (those produced by the sediment in equilibrium + 226 Ra activities; gray) per month and beach. When measured 222 Rn activities are over the equilibrium activities it is considered an indication of continental groundwater intrusion (dark gray), and when equilibrium activities are higher than those measured, it is considered an indication of recent seawater intrusion (light gray).
The opposite situation occurred in Panxón beach. Very high dissolved oxygen concentrations (> 100 µM) were found in the interior of the beach at the middle and lower stations, showing that beach oxygenation was not only caused by tidal and wave forcing. The inverse correlation between oxygen and salinity and the proximity of oxygen to 222 Rn in the PCA (Figure 6B) suggest that continental groundwater was, in part, the source of the observed high dissolved oxygen concentrations (Figure 5). Moreover, the highest inverse correlation between oxygen and salinity was observed in October (p < 0.05, r = 0.57, n = 18), when the local water balance was the highest of the studied period (Supplementary Table 1). This is consistent with the high oxygen levels observed in the nearby continental aquifers ( Table 1). The saltwedge, typically assumed anoxic (e.g., Robinson et al., 2007a), presented high oxygen concentrations in Panxón (Figure 5). The enhanced water transfer through the gravel layer seems to be the main cause. In Ladeira, low oxygen concentrations appeared at the bottom depths in most profiles (Figure 5), suggesting that the FDT was not completely anoxic. Oxygen (and nitrate) consumption in these groundwaters during their transit to the sampled zone of the STE may explain these observations, as none of the surrounding wells had low oxygen concentrations (Table 1). Additionally, a certain contribution of the nearby surface freshwater courses to the brackish waters found in the beach interior cannot be discarded, although these would not explain the observed porewater unsupported 222 Rn enrichment.
Organic matter to fuel heterotrophic processes in beaches can be sourced by either seawater infiltration (Ibánhez and Rocha, 2016) or fresh groundwater, in particular from contaminated aquifers (Yang et al., 2015). In this study, most of the sediment profiles of particulate organic matter decreased with depth (Figure 3). Thus, recirculated seawater may be the main source of particulate organic matter into the beach. The rapid infiltration of seawater during high tide and the accumulation of particulate organic matter in the shallower depths caused a higher consumption of oxygen near the surface than at the interior of the beach, where the only available form of organic matter is dissolved organic carbon (Rusch et al., 2000;Huettel et al., 2014;Ibánhez and Rocha, 2016). In consequence, the consumption of oxygen was reduced in the interior of the beach. This explains the preservation of the oxygen supplied by groundwater in the interior of Panxón and its decrease near the surface (Figure 5). The organic C: total N elemental ratio was lower than the Redfield ratio (106 C: 16N) in all the profiles of Panxón and higher than it for all the profiles taken in Ladeira (Figure 3). In Panxón, this suggests the presence of a large pool of adsorbed N, enhanced by the large content in carbonates of this beach (Rosenfeld, 1979).

Chemical Reactivity of the STEs
Conservative mixing profiles of nitrate, DIN and dissolved oxygen were calculated using salinity and 222 Rn as mixing tracers (Supplementary Figure 1) to support the analysis of the reactivity of Panxón and Ladeira beaches. For each porewater profile, the samples with the maximum and minimum values of the tracers (salinity or 222 Rn activity; Figure 4) and their associated nitrate, DIN and dissolved oxygen concentrations ( Figure 5) were used to construct mixing lines. The slope (m eq ) and the Y-intercept (n eq ) of these mixing lines were used to calculate the expected values of nitrate, DIN and dissolved oxygen (C mix ) for each porewater profile from the distribution of each conservative tracer by using the following equation: C mix = m eq C tracer + n eq where C tracer is the measured salinity or 222 Rn activity. The relative change between measured and conservative mixing nitrate, DIN and dissolved oxygen profiles ( Table 2 and  Supplementary Table 3) was calculated from the vertically integrated areas of each porewater profile ( Supplementary  Figure 2), using the following equation: where [C measured ] int is the vertically integrated nitrate, DIN or dissolved oxygen concentration for each measured profile and [C mix ] int is the vertically integrated nitrate, DIN or dissolved oxygen concentration for each conservative mixing profile. These integrals were calculated using the trapezoid approximation. Large differences were observed whether we used 222 Rn or salinity as mixing tracer (Table 2), mainly related to mixing processes (Spiteri et al., 2008;Robinson et al., 2018). Moreover, the relative changes of dissolved oxygen (Supplementary Table 3) were parallel to the relative changes of DIN ( Table 2). This could hinder the analysis since these positive values might be partly related to lateral transport, because oxygen cannot be produced in the beach interior. In the following subsections, conservative mixing profiles (Supplementary Figure 1), together with circulation tracer profiles (Figure 4) and porewater nutrient profiles (Figure 5) are examined individually to shed some light on the processes acting in the STEs of Panxón and Ladeira beaches. Errors are calculated from the analytical errors of measured profiles and standard errors of mixing lines.

Tidal Circulation Cell-Induced Nitrate Profiles in Ladeira Beach
Peak nitrate and dissolved oxygen concentrations appeared at the upper stations of February and May in Ladeira beach, with concentrations above the conservative mixing profiles (Supplementary Figure 1). There, the tidal circulation cell was actively pumping oxygen and organic matter into the sand (Figures 5, 8). However, neither PCA nor correlation analysis associated nitrate and salinity in Ladeira. Furthermore, DIN concentrations in the surface waters of the Ría de Vigo or in A Ramallosa Lagoon (Nogueira et al., 1997;Viana and Bode, 2013;Alonso-Pérez and Castro, 2014) are too low to explain the observed peak nitrate concentrations within the tidal circulation cell. Its origin seems to be the ammonification of the organic matter introduced into the sediment and fueled by the parallel transport of oxygen, followed by nitrification of the ammonium produced (Ibánhez and Rocha, 2017). Incomplete nitrification would thus explain the accumulation of nitrite observed near the surface in the upper profile taken in May (Figure 5). The absence of nitrate in depth further suggests that the significant nitrate levels measured in continental groundwaters in the studied area ( Table 1) might be consumed before arriving to Ladeira beach.

Nitrate Reduction and Groundwater Supply in Panxón Beach
The terrestrial fresh groundwater transported through the gravel layer was the main source of nitrate to Panxón beach. This is supported by the high correlation between 222 Rn activity and nitrate concentrations and the confluence of nitrate peaks in zones of enhanced transport (Figures 5, 7). Solute transfer was maintained even during the reduced continental groundwater flux of July, observable as an increase of nitrate and oxygen above conservative mixing porewater profiles (Supplementary Figure 1B). Scarce availability of fresh (labile) organic matter in the beach interior (Hartog et al., 2004;Slomp and Van Cappellen, 2004) would limit heterotrophic processes and may explain the presence of brackish waters with high oxygen and nitrate content in depth in Panxón. However, at shallower depths, the progressive availability or particulate organic matter trapped in the sediment matrix caused a decrease of oxygen and nitrate over conservative mixing profiles (Supplementary Figure 1B). Peaks of nitrite accompanied nitrate consumption (Figure 5), which seems to show active nitrate reduction in place (Ibánhez and Rocha, 2017). In the porewater profiles from July to October, organic matter seems to be depleted before nitrate was used as electron acceptor, which resulted in high nitrate concentrations under suboxic conditions (Figure 5). Furthermore, the ammonium released by organic matter remineralization was accumulated in the anoxic, shallower layers of some of the vertical porewater profiles (Figure 5).

Ammonium Distribution in the Suboxic Sediments of Ladeira Beach
Several processes (terrestrial fresh groundwater and anoxic seawater transport, dissimilatory nitrate reduction to ammonium, desorption from sediments, organic matter decomposition. . .) can support ammonium content in anoxic STEs (Slomp and Van Cappellen, 2004;Santos et al., 2008). In Ladeira beach, the slight increase in oxygen with depth and the accumulation of ammonium in the shallower depths ( Figure 5) suggest that ammonium is primarily originated from the remineralization of particulate organic matter trapped near the surface (Huettel et al., 2014;Ibánhez and Rocha, 2016), rather than carried by the FDT or infiltrating seawater. The vertical distribution of ammonium in the porewater seems to be also linked to porewater circulation. The increased mixing observed in July could have been the cause of the nearly homogenous vertical distribution of ammonium in the porewater (Figure 5 and Supplementary Figure 1C). Increased residence times and anoxic conditions could be linked to DIN consumption in the middle station of May, as conservative mixing profiles above measured DIN profiles suggest (Supplementary Figure 1C).
In the lower sampling station during February and May, high ammonium concentrations (>30 µM) were observed throughout the entire porewater profile (Figure 5). The accumulation of seaweed in the low intertidal, observed during these surveys, has been pointed to fuel benthic heterotrophic processes (Santos et al., 2008), leading to the leaching of ammonium into the interior of the beach.

Porewater Circulation and Biogeochemical Reactivity of STEs
The suitability of calculating conservative mixing profiles for the study of the chemical reactivity in STEs has been previously discussed (Spiteri et al., 2008;Rocha et al., 2009). This approach assumes vertical, steady state porewater transport (Rocha et al., 2009). As a result, lateral transport cannot be differentiated from chemical reactivity when following this approach and thus, it is not recommended if the previous assumptions are not fulfilled (Spiteri et al., 2008). In this study, the two endmembers used to define the conservative mixing porewater profiles are local, i.e., the maximum and minimum values of 222 Rn and salinity at each porewater profile (Figures 4, 5). Besides, to reduce the possible misinterpretation of deviations of nitrate and DIN concentrations from the calculated mixing lines, conservative mixing profiles of dissolved oxygen were also calculated to identify lateral transport, as oxygen cannot be produced in the interior of the beach. The results discussed in the previous sections showed that oxygen, DIN and nitrate conservative mixing profiles agreed well with the circulation processes defined in section "Porewater Circulation in the Ladeira and Panxón STEs." Although the apparent reaction magnitudes shown in Table 2 should be taken with caution due to the aforementioned limitations, the approach followed permitted to unveil the sharp biogeochemical zonation of the studied STEs. This study establishes a necessary link between internal beach hydrodynamics and the biogeochemical processes that regulate the export of nutrients from STEs to the coast. Porewater circulation is the main driver of the strong redox differences found between the studied STEs, thus strongly determining the reactivity of redox-sensitive N compounds (Slomp and Van Cappellen, 2004;Spiteri et al., 2008;Ibánhez et al., 2013;Rocha, 2016, 2017). Our results highlight the role of the tidal circulation cell over the biogeochemical reactivity of STEs and the seasonality of N speciation observed. Nevertheless, our results further suggest that the different redox conditions found in the FDT are still the main factor to determine N transformations within the STEs. Sediment heterogeneity has been suggested to be particularly important in the establishment of these redox conditions in the FDT (Geng et al., 2020a,b), as it was the main cause of the biogeochemical differences between the two STEs studied here.

CONCLUSION
The biogeochemical reactivity of two contrasting STEs has been studied for the first time in the Ría de Vigo. The three water bodies that typically shape tidally driven STEs were identified in both beaches using salinity, 222 Rn and 226 Ra as circulation tracers. Large gradients of nitrate, oxygen, and ammonium were observed at both STEs, partially related to porewater circulation. The tidal circulation cell acted as a pump of oxygen and organic matter into the interior of the beaches, inducing the aerobic respiration of organic matter and the nitrification of the ammonium produced. Thus, the tidal circulation cell regulates the transport of regenerated nutrients to the coast. Terrestrial fresh groundwater transported nutrients from the coastal aquifer into the interior of the beaches. Its role in nitrogen processing in Panxón and Ladeira STEs depended principally on the dissolved oxygen levels present at the FDT, which also depended on the different porewater circulation patterns observed in those beaches. Suboxic conditions in the terrestrial fresh groundwater of Ladeira beach caused the dominance of nitrogen reduction processes and the accumulation of ammonium. On the contrary, oxic conditions prevented the consumption of nitrate in the interior of Panxón beach. The layer of increased permeability in depth, besides of having completely altered the circulation structure of the STE, seems to play an important role in enhancing the transport of landderived nitrate to the coast by limiting nitrate reduction along the flow path but also by fueling nitrification in the beach interior. The degradation of the particulate organic matter trapped near the surface during surface water infiltration (i.e., high tide) caused the consumption of oxygen and groundwater-borne nitrate and the accumulation of ammonium in the shallower depths. As a consequence, our results suggest that the preferential transfer of solutes through the gravel layer together with the oxidation of organic matter within the tidal circulation cell causes Panxón beach to be a net source of N to the coast. In opposition, the suboxic conditions of Ladeira might make this beach a net sink of N.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
JSPI and CR conceived the study. JSPI, EC-M, and XÁ-S participated in the field and laboratory work. EC-M and JSPI elaborated the first draft of the manuscript, which was enriched by contributions of XÁ-S and CR. All authors revised and approved the manuscript.

FUNDING
This research was supported by the SUBACID project (SUBmarine groundwater discharge (SGD) impact on coastal ACIDification processes in contrasting Atlantic Shores: towards securing ecosystem services and food production) that has received funding from the Irish Research Council and the European Union's Horizon 2020 Research and Innovation