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). 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 medium-high 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.

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₂), radon (²²²Rn), radium (²²⁶Ra), and dissolved inorganic nitrogen (DIN; NO₃^–, NO₂^–, and NH₄⁺). 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.

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 ²²²Rn-free water. Readings were corrected for the internal ²²²Rn decay from the time of collection to analysis. To obtain ²²⁶Ra activities, purged samples were stored until ²²²Rn and ²²⁶Ra reached secular equilibrium (∼30 days, i.e., more than 5 times the ²²²Rn half-life) and measured again. The equipment was purged with an activated carbon ²²²Rn trap to minimize background activities. The detection limit of our procedure determined from the regular analysis of blanks (n = 36) made with ²²²Rn-free water (250 mL) was 36 Bq m^–3. Thus, the detection limit of ²²²Rn and ²²⁶Ra in our diluted samples is estimated at 128 Bq m^–3. Porewater ²²²Rn equilibrium activities were measured with sediment incubations (∼200 g dry weight; n = 9 for each beach) run with ²²²Rn-free water (Colbert and Hammond, 2008). After 30–40 days of incubation, ²²²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 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 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¹ 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 recharge-discharge 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) were obtained from the SIMAR Point (3014003; 42°7.8′N 8°55.2′W) available from Puertos del Estado². Average and standard deviation of sea surface salinity and dissolved oxygen data (Supplementary Table 1) were calculated with the data from the INTECMAR repository³, from the buoy located at the port of Baiona (42°7.26′N 8°50.77′W).

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; Figure 3). C:N molar ratio showed higher average and range in Ladeira (from 7.0 to 21.2; Figure 3) than in Panxón (which had a range of 2.1–4.8; Figure 3).

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).

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 ²²²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 ²²²Rn activities were, in most cases, greater than these ²²²Rn equilibrium activities, suggesting external inputs of ²²²Rn to the beach porewaters. Exceptions are found in some sampled porewaters near the sediment surface (Figure 4). Average activities of dissolved ²²⁶Ra were one order of magnitude lower than the average activities of ²²²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.

In most porewater profiles ²²²Rn activities were the lowest at shallow depths. Conversely, ²²²Rn peaks of up to 11234 ± 353 Bq m^–3 in Panxón and 13551 ± 732 Bq m^–3 in Ladeira (both at the lower station during October) were registered in depth (Figure 4). Besides, increases in ²²²Rn activity were followed by decreases in salinity at various depths (e.g., at the lower station in October at Panxón). This inverse relationship between salinity and ²²²Rn activity was significant in May (p < 0.05, r = −0.50, n = 17) and October (p < 0.01, r = −0.62, n = 18) in Panxón and in May (p < 0.01, r = −0.68, n = 18) in Ladeira beach. Peaks of ²²⁶Ra activity (> 600 Bq m^–3) appeared in some porewater profiles (e.g., at the lower station in July, in Ladeira; Figure 4). ²²⁶Ra activities were correlated with salinity in July (p < 0.05, r = 0.54, n = 18) in Panxón, and in May (p < 0.05, r = 0.54, n = 18) and July (p < 0.001, r = −0.82, n = 18) in Ladeira.

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).

Spearman correlation analysis showed significant relationships between nutrient and internal porewater circulation tracers. In Panxón beach, nitrate and ²²²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 ²²²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 ²²²Rn and salinity there and suggests a relationship between nitrate enriched, oxygenated porewater and ²²²Rn-enriched, 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 ²²²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 ²²²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), ²²⁶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 ²²²Rn (−0.59), salinity (0.52) and ammonium (0.47). Both PCAs confirmed the relationships observed through Spearman correlation analysis.

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 ²²⁸U decay chain (Tollefsen et al., 2016), such as ²²⁶Ra and ²²²Rn. Therefore, groundwaters contained within these rocks are commonly enriched in these radioisotopes, particularly ²²²Rn since ²²⁶Ra shows very low solubility at low salinities (Gonneea et al., 2008).

Porewater ²²²Rn equilibrium activities estimated from sediment incubations together with the measured dissolved ²²⁶Ra activities (∼130–1,600 Bq m^–3) are not sufficient to explain the observed ²²²Rn enrichment in the sampled porewaters. Furthermore, the negative correlation between ²²²Rn and salinity found in some of the surveys and the opposite PCA scores between ²²²Rn and salinity in Ladeira (Figure 6C) point toward a freshwater origin for this ²²²Rn enrichment. Ibánhez et al. (2019) measured radon activities at various freshwater bodies near Baiona Bay (Table 1). There, ²²²Rn activities found in the Miñor River (Figure 1; 2.7 ± 0.1×10³ Bq m^–3 in winter/spring and 2.5 ± 0.1×10³ Bq m^–3 in summer) are one order of magnitude below the maximum ²²²Rn activities observed in Ladeira. On the other hand, ²²²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³ Bq m^–3 and 802 ± 28 ×10³ Bq m^–3; Table 1). Thus, considering that rainfall has low ²²²Rn content, the high radon levels found in beach porewaters can only be sourced from groundwater from the coastal aquifers.

Assuming that ²²²Rn is in secular equilibrium within the coastal aquifer and that no net ²²²Rn loss occurred during the transit to the studied beaches except from ²²²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. ²²²Rn activities in excess (Rn_excess) are defined as the difference between measured porewater ²²²Rn activities and those supported by porewater ²²²Rn equilibrium activities plus porewater dissolved ²²⁶Ra activities. We further assume that measured porewater ²²²Rn equilibrium activities in shallow sediments (0–10 cm) are representative of all sampling depths (0–180 cm), since sediment composition is very homogeneous in the studied beaches, with the exception of the gravel layer of Panxón beach (Álvarez-Vázquez et al., 2003; Perez-Arlucea et al., 2005; González-Villanueva et al., 2007). The ²²²Rn activities in excess recorded during February and July were used to calculate groundwater transit times (t) as follows:

where λ is the decay constant of ²²²Rn (0.18 d^–1) and Rn_extrapolated is the ²²²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:

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 steady-state, 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 ²²²Rn in water and v the groundwater velocity. A D_molec of 1.102×10⁻⁵ cm²/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 ²²²Rn enrichment in the beach porewaters) and Rn_aquifer values (< 50 × 10³ Bq m⁻³; 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 ²²²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.

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 ²²²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 ²²²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 ²²²Rn activities lower than porewater ²²²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, 2014, 2018; Heiss 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 ²²²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 ²²²Rn activities found in October at Panxón beach. Adsorbed ²²⁶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 ²²⁶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, ²²²Rn-enriched waters observed during October in Panxón are originated from seawater previously intruded into the aquifer during the summer.

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.

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 ²²²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 salt-wedge, 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 ²²²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).

Conservative mixing profiles of nitrate, DIN and dissolved oxygen were calculated using salinity and ²²²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 ²²²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:

where C_tracer is the measured salinity or ²²²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 ²²²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.