Transport and transformations of groundwater-borne carbon discharging through a sandy beach to a coastal ocean

Abstract

Submarine groundwater discharge (SGD) is now recognized as an important coastal process affecting local and regional sources of solutes to coastal oceans. The objective of this study was to compare SGD estimates through Martinique Beach (Îles-de-la-Madeleine, QC, Canada) using different sampling techniques in order to quantify SGD fluxes of dissolved organic and inorganic carbon (DOC, DIC) exported to a coastal embayment. These fluxes were estimated using hydrogeological and geochemical methods, including direct measurement via seepage meters and hydrogeological estimation using Darcy’s law. Darcy estimates led to a fresh SGD flow ranging from 0.020 m³/s at the shoreline to 0.030 m³/s at the beach face. Direct measurements by seepage meters revealed higher SGD flows, with a mean value of 0.090 m³/s. These SGD flows were mainly composed of recirculating seawater flushing by the falling tide. The unconfined aquifer of Martinique Beach is a biogeochemically reactive zone, where DOC and DIC are nonconservative. The groundwater residence time (~32 days), low oxygen conditions (~20%) and high DOC concentrations (~ 2 mmol/L) represent a suitable environment for biogeochemical reactions to occur, and subsequently alter groundwater-borne carbon concentrations. We calculated a fresh-groundwater carbon load of 147 kg/d of DIC and 27 kg/d of DOC at the beach face. Even though recirculating seawater dominated the volume of total SGD, fresh SGD was an important carbon pathway, accounting for 12 and 20% of total DIC and DOC, respectively. This site-specific study is the first attempt to estimate volumetric and chemical groundwater fluxes to a coastal Canadian ocean, and demonstrates the strong interaction between fresh groundwater and coastal systems.

Les décharges d’eau souterraine (SGD pour submarine groundwater discharge) sont courantes en milieu côtier et sont maintenant reconnues comme des sources locales et régionales de solutés à l’océan côtiers. L’objectif de cette étude était de comparer différentes techniques pour estimer les SGD à travers la plage de la Martinique (Îles-de-la-Madeleine, QC, Canada) et ainsi pouvoir évaluer les flux de carbone organique et inorganique dissous (DOC et DIC) exportés à la côte. Des données aussi bien hydrogéologiques que géochimiques ont été utilisées pour évaluer ces flux, en utilisant notamment des mesures directes avec des chambres d’exfiltration (seepage meter) et des mesures hydrogéologiques en se basant sur la loi des écoulements de Darcy. Les estimations faites à partir de Darcy donnent des flux d’eau douce entre 0,020 m³/s à la ligne de rivage et 0,030 m³/s dans la zone de décharge de la plage. Les mesures directes donnent des valeurs plus élevées, autour de 0,090 m³/s mais ces flux incluent l’eau de mer de recirculation advectée par la pompe tidale. L’aquifère non confiné de la plage est aussi une zone biogéochimiquement active où le DOC et DIC montrent un comportement non-conservatif. Le temps de résidence de l’eau douce (autour de 32 jours), les faibles concentrations en oxygène dissous (DO ~20% de saturation) et les fortes concentrations en DOC (~2 mmol/L) permettent le développement de conditions idéales aux processus de respiration et à la transformation du carbone apporté par les écoulements souterrains. Ces transformations mènent à des flux souterrains de 147 kg DOC/j et 27 kg DIC/j dans la zone de décharge. Et même si l’eau de mer de recirculation domine le volume total de SGD dans ce site, l’écoulement d’eau douce est une importante voie de transport de carbone à l’océan et représente respectivement 12 et 20% des apports en DIC et DOC. Bien que ces résultats soient spécifiques au site d’étude, ils représentent la première tentative d’estimation des flux, qu’ils soient volumétriques ou chimiques, dans l’océan côtier canadien et soulignent la forte interaction entre les eaux souterraines et les milieux côtiers.

Introduction

Coastal aquifers constitute the major source of potable water for many regions throughout the world (Bear et al. 1999). Their importance as a fresh water source for coastal populations as well as the impacts of their submarine discharges to coastal zones highlight the need to implement an integrated water resource management plan in coastal zones. This is particularly important in island systems, where the connectivity between groundwater and seawater is strong. Submarine groundwater discharge (SGD) has been increasingly identified as an important process that coastal managers should consider in any water resources management plan (Moore 1996; Taniguchi et al. 2002; United Nations Educational, Scientific and Cultural Organization (UNESCO) 2004; McCoy and Corbett 2009). This is particularly important where booming coastal populations induce numerous changes in landscapes and ecosystems as well as increase pressure on coastal resources. In locations around the world, deteriorations in coastal water quality (e.g. bacterial proliferation, harmful algal blooms, hypoxia, acidification, fish and shellfish mortality) have been attributed to the quality and volume of groundwater discharge to the ocean (Valiela et al. 1978, 1980; Johannes 1980; Laroche et al. 1997; Paerl 1997; Gobler and Boneillo 2003; Hwang et al. 2005; Y.W. Lee and Kim 2007; Y.W. Lee et al. 2009; G. Kim et al. 2011). These disturbances alter coastal ecosystems and endanger species that have ecological, cultural and economic importance.

Moore (2010, p. 60) defined SGD as “any and all flow of water on continental margins from the sea bed to the coastal ocean, with scale lengths of metres to kilometres, regardless of fluid composition or driving force.” SGD includes both fresh groundwater flowing directly to the sea and seawater recirculating through the seabed. SGD occurs in many environments along continental margins, usually as submarine springs or scattered seeps through permeable seafloor and beach sediments. This seeping occurs whenever the water table remains above sea level and coastal aquifer recharge is efficient (see Taniguchi et al. 2006 and references therein).

Sandy sediments make up more than 40% of the world’s coastlines (Riedl et al. 1972; Brown and McLachlan 2002) and are characterized by high to very high biogeochemical reaction rates despite their low particulate organic carbon content (Boudreau et al. 2001; Anschutz et al. 2009). Multiple physical forces that drive subsurface flow and residence times control the reaction rates and thus the transformations in this active biogeochemical zone (e.g. Boufadel 2000; Horn 2002, 2006; Michael et al. 2005, 2011; Charette et al. 2007; Robinson et al. 2007; Cable and Martin 2008). Santos et al. (2012) identified at least 12 independent drivers that include both terrestrial (e.g. hydraulic gradient, seasonal oscillation of the water table) and marine (e.g. wave and tidal pumping, ripple and bed form migration, bioirrigation) processes that interact in complex ways. All of these processes force flow across the water–sediment interface, albeit at different spatial and temporal scales, and consequently influence the export of groundwater-borne, recirculated and newly formed compounds to the sea. The chemistry of SGD changes during its passage through the aquifer and permeable coastal sediments, but how these transformations affect the concentrations and exported fluxes of groundwater-borne carbon is not yet well understood. In particular, carbon is a key element of coastal eutrophication, and the exported fluxes of dissolved inorganic and organic carbon (DIC and DOC, respectively) via SGD are still not well quantified (Caï et al. 2003; Goñi et al. 2003; Beck et al. 2007; Dorsett et al. 2011; Moore et al. 2011; Cyronak et al. 2013. Kim et al. 2013; Maher et al. 2013). Such quantifications are vital for managing coastal ecosystems, which include the groundwater–ocean continuum. This is particularly important in northern systems because these regions are experiencing – and are expected to continue to experience – important changes in hydroclimatic and terrestrial carbon storage (Hinzman et al. 2005; Roulet et al. 2007).

In the Îles-de-la-Madeleine archipelago, which is located in the southern Gulf of St. Lawrence (QC, Canada; Figure 1a, b), groundwater is the only source of fresh and drinking water. The absence of major streams limits the discharge of groundwater to continental surficial ecosystems, and groundwater rapidly flows to coastal ecosystems (residence time < 1 year; Madelin’Eau 2004). This is a typical coastal aquifer system where the quality of groundwater can affect the ecological health of adjacent embayments.

Figure 1. (a) Location of Quebec, (b) the Îles-de-la-Madeleine and (c) the archipelago’s main island (Cap-aux-Meules). The study area (from le Chemin des Chalets to Gros-Cap) and Martinique Beach are shown. The dashed line corresponds to the cross section presented in Figure 2a. The simplified geology was adapted from Brisebois (1981) (data source ^©Gouvernement du Québec).

The objective of this study was to compare SGD estimates using two different sampling techniques in order to quantify SGD fluxes of DOC and DIC exported to the coastal zone through a sandy beach. More specifically, hydrogeological and geochemical data are used to estimate carbon fluxes. Data were collected by direct measurement via seepage meters, and hydrogeological estimations were made using Darcy’s law. This study combines knowledge of marine chemistry and hydrogeology, and proposes an integrated view of water resources in the groundwater–ocean continuum.

Material and methods

Study area

The Îles-de-la-Madeleine are part of the Maritimes Permo–Carboniferous Basin included in the Northeastern Appalachian Geological Province (Figure 1; Brisebois 1981). The archipelago comprises a group of 15 islands, seven of which are joined as tombolos. The main island is part of the Mississippian–Windsor group, which is an interbedded assemblage of volcanic and volcanoclastic rocks flanked by sandstones of the Permian Inferior period (the Cap-aux-Meules Formation; Brisebois 1981). This formation constitutes the main aquifer in the archipelago; the mean transmissivity (~4 × 10⁻³ m²/s) and mean recharge (~230 mm/y, i.e. ~25 to 30% of the annual snow and rain precipitation; Leblanc 1994; Madelin’Eau 2004) support high pumping rates. Groundwater flows from the central volcanic core recharge through the unconfined Permian sandstone aquifer (Figure 2) and discharges to the sea, both directly and through overlying Quaternary deposits.

Figure 2. (a) Schematic cross section through the main geologic units of the main island (Cap-aux-Meules). The municipal (P4, P5, PU6) and private (PM) water wells are located along the geological transect. The position of the salt water/fresh water interface (shaded line) was estimated according to the Ghyben–Herzberg relation based on mean annual piezometric levels recorded in the municipal and private wells (corrected to mean sea level). (b) View of the sampling transect in Martinique Beach showing locations of piezometers (Pz_1–3), multi-level samplers (M_1–7), and seepage chamber deployments. The different morphology units of the beach also are reported.

Most of the archipelago’s coasts are sandy barrier systems with fine-grained sands resulting from erosion of Holocene deposits and eolian transport. Sand grain size is about 300 μm, and sediments are mainly composed of silicate fragments (Chaillou et al. 2014). The underlying sandstone aquifer consists of fine red-orange sands (~100 μm) composed of silicate and aluminosilicate with iron (Fe)-coated silicate grains. The archipelago is located near an amphidromic point (mean and maximum tidal ranges are 0.7 m and 1.1 m, respectively), with mixed semi-diurnal and diurnal tides. The wave regime is characterized by locally generated waves with a mean height of 0.87 m (Owens and Frobel 1977).

Martinique Beach is a typical sandy barrier system on the southwest coast of the archipelago. It has been classified as a low-energy beach with a low relative tidal range, indicating wave control of beach morphology (Masselink and Short 1993; Jackson et al. 2002). A high-definition elevation profile resulting from differential global positioning system (DGPS) surveys revealed three morphologic zones (see beach morphology in Figure 2b). These zones include a vegetated beach terrace on the backshore (~10 m wide) composed of fine-grained sands, a beach face composed of fine-grained sands with a mean slope of 6° (~15 m wide), and a large low-tide terrace (~50 m) composed of very fine sands with a lower gradient (~2°) and persistent micro-topographic features under non-storm wave conditions; all of these are characteristic of low-energy beaches (Jackson et al. 2002). The discharge zone is not well characterized, but a recent radon-222 (²²²Rn) survey suggests a narrow discharge zone (~100 m) from the beach face (Lemay-Borduas et al. 2014). In the absence of surface streams, a recently developed groundwater flow model shows that ~96% of the fresh groundwater volume seeps to the coastal ocean via SGD each year, leading to discharges of around 4300 m³/d in Martinique Bay (Madelin’Eau 2004; Chaillou et al. 2013).

Approach

Two approaches were used to investigate SGD through the sandy beach of Martinique Bay to quantify exported fluxes of DIC and DOC. These techniques included groundwater flow estimations using Darcy’s law and direct measurements via seepage meters installed in the low-tide terrace between 20 May and 7 June 2013. Hydrogeological data from municipal and private water wells were used to estimate the regional groundwater flow perpendicular to the shoreline (Figure 2a; Madelin’Eau 2004, 2007, 2011). Based on Darcy’s law, and assuming the hydraulic gradient component parallel to the coastline is null, groundwater fluxes in the Permian aquifer (q_inland in cm/d) were calculated as:(1)

where K is the hydraulic conductivity of the Permian aquifer and gradh is the component of the hydraulic gradient perpendicular to the shoreline. The fresh groundwater fluxes were also estimated by Darcy’s law in the beach aquifer (q^beach in cm/d) using data from three monitoring piezometers placed perpendicular to the shoreline (Figure 2b). Seepage meter experiments were conducted in the low-tide terrace of the beach (see locations in Figure 2b) to measure the volume of water that seeps from the seabed. Despite criticisms concerning pressure-induced flow (see Shinn et al. 2002 and references therein), this technique provides a flux at a specific time and site from a limited seabed area (Corbett and Cable 2003; Burnett et al. 2006). The two techniques have advantages and disadvantages, but their combination provides valuable insight into the SGD regime, which includes both fresh groundwater and recirculating seawater (see Burnett et al. 2006 and references therein).

Sample collection

Three piezometers were inserted, in the beach terrace (P_Z1), the beach face (P_Z2) and the low-tide terrace (P_Z3; Figure 2b). Piezometers were constructed of 2.5-cm polyvinyl chloride (PVC) pipes (Ø = 5 cm) and were installed to a depth of 2.2 m using a manual auger; a 10-cm solid PVC riser was used for the casing. Water levels were monitored manually and with automated level loggers (HOBO models). Tides and waves were recorded at Pz₃ and Pz₂, and water table levels were corrected to obtain mean values of the water table over the sampling period and to estimate the component of the hydraulic gradient perpendicular to the shore. Guelph permeameter measurements (N = 15) were conducted to estimate local hydraulic conductivity of the beach’s surficial sandy aquifer. Values are reported in Table 1.

Table 1. Summary of hydrogeological values used to estimate groundwater fluxes and submarine groundwater discharges assuming a 1200-m sandy beach shoreline. K is the hydraulic conductivity (in m/d); gradh is the change in head per unit length in the direction of flow considering an isotropic aquifer. The groundwater fluxes are reported in term of q (in cm/d) and Q (in m³/s) based on Darcy’s law. Mean values and standard deviations (mean ± SD) are presented. Data from Madelin’Eau (2004) were used to estimate groundwater fluxes in Permian aquifer, in inland water wells.

			K (m/d)		gradh	q (cm/d)		Q (m^3/s)
		Method	Range	Mean		Range	Mean	Range	Mean
Inland
	Well	Darcy	n.d.	1.8	0.003	0.7–3.6	2.1	0.007–0.036	0.021
Beach
	Beach face	Darcy	7–20	11.4 ± 4.4	0.037	14–42	30.4 ± 11.8	0.008–0.046	0.029 ± 0.009
	Low-tide terrace				0.012
Seabed
	Seepage meter	Direct measurement				3–10	7.5 ± 2.0	0.041–0.130	0.085 ± 0.046

n.d.: no data.

The seepage meter design followed D.R. Lee (1977), using the top 30 cm of a standard 55-US gallon (208-L) plastic drum fitted with a sample port and a plastic collection bag. Seepage meters were carefully pushed to a depth of 15 cm into the seabed and allowed to equilibrate for ~4 hours. Seepage bags (2000 mL) were prefilled with 100 mL of seawater and the air was removed; these bags were deployed for 1.5 to 4.0 hours and then carefully removed. Collected water volumes were measured to the nearest millilitre (using 1000-mL graduated cylinders) to determine the groundwater seepage fluxes (i.e. the collected volume divided by the time interval and the seepage meter area). Because seepage meters can be deployed only in calm wave and wind conditions, only seven experiments were performed (22, 24, 28, 31 May; 2, 3, 5 June 2013). Seepage meters were deployed at different locations in the low-tide terrace of the beach, a few metres from the shoreline (Figure 2b).

The flux of exported carbon associated with groundwater discharge was calculated based on DIC and DOC concentrations measured during the transit of groundwater to the ocean. Water samples were collected during the low spring tide period in May and June 2012 and 2013. A ~35-m transect of multi-level samplers (M_1–7; Figure 2) was installed on the beach. The multi-level samplers consisted of 2.5-m PVC pipes with eight pores connected to flexible Tygon® tubing, as described in Martin et al. (2003). Multi-level samplers were designed to collect water at 10, 30, 50, 80, 110, 150, 190 and 230 cm below the beach surface. In 2012, beach groundwater samples were collected at M₁, M₂, M₃ and M₆ (Figure 2b). In 2013, the focus was the intertidal zone, where beach groundwater and recirculating seawater most likely discharge (M₃ to M₇; Figure 2b). Beach groundwater was pumped with a peristaltic pump into an on-line flow cell where specific conductivity (SpCond) and dissolved oxygen saturation (DO) were monitored using a calibrated multiparametric probe (600QS, YSI Inc.). After these parameters stabilized, groundwater samples were collected and preserved for DOC and DIC analysis. DOC samples were taken using 60-mL acid-clean syringes and rapidly filtered (Whatman GF/F 0.7-μm filters), acidified to pH = 2, and preserved in acid-clean glass tubes. All filters and glass tubes had been pre-combusted at 500°C for 4 h. Samples were kept in the dark at 4°C until analysis. Samples for δ¹³C-DIC were collected by filling 250-mL glass bottles with no headspace, preserved with 0.5 mL of saturated HgCl₂ and stored at 4°C. δ¹³C-DIC was analyzed for the 2012 and 2013 samples while DIC concentrations were measured only from samples collected at M₃, M₄ and M₇ in 2013.

The fresh and marine end-member concentrations of DOC and DIC were systematically measured. Samples from the Permian sandstone aquifer were collected as described above from private and municipal water wells located 50 to 2000 m inshore from the multi-sampler transect (see PM, P5, P4 and PU6 wells in Figure 2a). Deep samples of bay water (~50 cm above the seabed) were collected 1 km offshore in Martinique Bay. The water samples were collected from a small boat using a submersible pump connected to an on-line flow cell.

DOC samples were analyzed by high temperature combustion (HTC) using a Total Organic Carbon (TOC) analyzer (TOC-V_cpn; Shimadzu) based on the method proposed by Wurl and Sin (2009). The detection limit is usually 0.05 mg/L and the analytical uncertainties were < 2% for concentrations higher than 1 mg/L. Dissolved inorganic carbon (DIC or ΣCO₂ = CO₂ + HCO₃^– + CO₃^2–, also called total carbonate) were analyzed using the flow injection analysis (FIA) technique based on the method proposed by Hall and Aller (1992). A gas-permeable membrane is used to move CO₂ from acidic reagent streams to a basic receiving stream and conductivity detector. The calibration was done using a reagent-grade sodium bicarbonate (NaHCO₃) solution, and the conductivity response was linear for samples over DIC concentrations from < 0.2 to 20 mmol/L with a precision of ~1%. The δ¹³C-DIC samples were analyzed using a gas chromotography coupled to an isotope ratio mass spectrometry (GC-IRMS) (GV Instruments). An internal calibration was performed using two internal standards to obtain a two-point calibration. Results are reported in parts per thousand in the δ notation relative to relative to Vienna Pee Dee Belemnite (VPDB) (standardized to TS-Limestone [NEBS19] and lithium carbonate [LSVEC] International Atomic Energy Agency (IAEA) referenced materials) with a precision of ± 0.1‰.

Results

Groundwater flux and SGD estimates

The survey of water table levels in the different transect piezometers revealed that the hydraulic head difference perpendicular to the shore is different in the inland part of the transect and in the beach. The multiannual average component of the hydraulic gradient from P5 to PU6 has been estimated at 0.0035 (Madelin’Eau 2004, 2007). During the survey in June 2013, the head gradient gradh (i.e. the change in head per unit length in the direction of flow considering an isotropic aquifer) is one order of magnitude higher in the beach aquifer, with values between 0.037 and 0.012. The mean conductivity value K of the unconfined Permian sandstone aquifer is 1.8 m/d (Madelin’Eau 2004, 2007). In surficial beach sand sediments, the K values measured with the Guelph permeameter varied between 7 and 20 m/d, with a mean value of 11.4 ± 4.4 m/d. The spatial variability of K in the beach sediments results from the heterogeneity of surficial sediment: Chaillou et al. (2014) reported the presence of an old forest soil layer on top of the sandstone aquifer caused by a recent sea-level rise. The range of values is consistent with those reported for similar lithology elsewhere on the Îles-de-la-Madeleine archipelago (5–9 m/d; Sylvestre 1979).

In the landward part of the studied transect, in the Permian aquifer, q_inland varied between 0.7 and 3.6 cm/d, with a mean value of 2.1 cm/d. As fresh water approaches the coastline, groundwater flux is expected to increase because the groundwater flows through a smaller area. Using hydraulic conductivity from the sandy beach, higher q_beach fluxes were calculated: values ranged from 14 to 42 cm/d, with an average of 30.4 ± 11.8 cm/d. Here residence time is defined as the time required for groundwater to transit horizontally through the system, in a manner similar to surficial estuaries. If the beach material is assumed to be homogeneous and isotropic (as is often assumed in flow modeling studies; see e.g. Raubenheimer et al. 1999 and Abarca et al. 2013), the maximum residence time in the surficial sediment is 32 days for a beach transect of 35 m (from M₁ to M₇; Figure 2b).

To convert these results to volumetric flux, which is the discharge per unit width of the aquifer (in m³/s), the cross-sectional flow area was determined using global positioning system (GPS) measurements of 1200 m of shoreline, from the sandy beach east of Le Chemin de la Martinique to the rocky shore of Gros Cap (Figure 1c). Furthermore, based on the Ghyben–Herzberg and Glover relationship (Cooper 1964), the fresh water/salt water interface is reached about 73 m below the water table of the aquifer at the PM well (the well nearest the coastline; Figure 2a). Hence, a flow depth of 73 m is used to estimate the freshwater flux at the coastline. The volumetric groundwater flux from the Permian aquifer to the shoreline (Q_inland) is estimated to be between 0.007 and 0.036 m³/s, with a mean flux of 0.021 m³/s. Also based on the Ghyben–Herzberg and Glover relationship (Cooper 1964), the top 3.2 m of the aquifer at the top of the intertidal zone (Pz₂) is fresh (except for the narrow surficial saltwater lens). Hence, we used a depth of 3.2 m to estimate the volumetric groundwater flux through the beach face in May 2013. This led to a beach flux (Q_beach) between 0.008 and 0.046 m³/s, with a mean value of 0.029 ± 0.009 m³/s.

Direct measurements revealed a mean seepage flux (q_seep) of 7.5 cm/d across the low-tide terrace seabed (Table 1). Seepage fluxes are not uniform over time: they tend to increase with the falling tide to values > 10 cm/d, and to decrease to ~3 cm/d at low and rising tides. However, given the limited number of measurements (N = 7), no clear relationship was observed between seepage fluxes and tidal levels. Assuming a 1200-m shoreline and a 100-m discharge area (Lemay-Borduas et al. 2014), SGD is estimated to be between 0.041 and 0.130 m³/s, with an average value of about 0.085 ± 0.046 m³/s. Table 1 summarizes these results.

Groundwater chemistry

End-members

Conductivity of the shoreline ocean water, i.e. the marine end-member, was measured at ~47 mS/cm, and the temperature was around 11°C (Chaillou et al. 2014). The dissolved oxygen (DO) saturation was > 100%, indicating high saturation due to wave and wind activity. Dissolved carbon concentrations were low, with DOC between 0.13 and 0.16 mmol/L and DIC between 2.1 and 2.5 mmol/L. The δ¹³C-DIC was slightly more variable, with values from –4.2 to –7.1‰ (mean –5.65‰). Groundwater conductivity in municipal and private wells, i.e. the freshwater groundwater end-member, was around 0.3 mS/cm, and temperature was measured at 7°C (Chaillou et al. 2024). DOC measured in well water was low, 0.30 ± 0.08 mmol/L, and minimal values of DIC were between 0 and 0.4 mmol/L. The δ¹³C-DIC was variable, with values from +8.0 to –1.9‰ (mean +3.05‰).

Beach groundwater

Two-dimensional profiles of beach groundwater chemistry are presented in Figure 3. The sample depths in these profiles are relative to mean sea level (i.e. 0 m depth), and contour lines were derived by the linear spatial interpolation (kriging method) of data points. The interpolation model adequately fits the empirical data set (95% confidence level). The two-dimensional conductivity profiles revealed the presence of fresh groundwater (conductivity < 5 mS/cm) in the landward part of the transect. The fresh groundwater lens diminished gradually towards the seaward discharge region below the intruding narrow saltwater lens (conductivity > 20 mS/cm) located high in the intertidal zone. Vertical conductivity profiles revealed limited mixing between the fresh groundwater layer and the surficial saltwater lens, which is characteristic of a microtidal system dominated by fresh groundwater discharge (Moore 1999; Robinson et al. 2007; Heiss and Michael 2014). While groundwater was suboxic, with DO around 20% saturation, the surficial saltwater lens was the most oxygenated layer, with DO saturation > 60%, indicating inputs of oxygen from seawater due to mixing.

Figure 3. Conductivity (mS/cm; a, c), and dissolved oxygen saturation (%; b, d) in the beach groundwater. The left panels (a, b) present data measured in 2012 along a 35-m transect of multi-level samplers; the right panels (c, d) present data measured in 2013 in the intertidal zone (dashed rectangle in a and b panels). The beach morphology has been done by differential global positioning system (DGPS) measurements. The depth is relative to mean sea level (i.e., 0 m depth) and the contour lines were derived by linear spatial interpolation (kriging method) of the data. The interpolation model reproduced the empirical data set very well (95% confidence level). Data were modified from Chaillou et al. (2014).

The concentrations of DOC and, DIC, and the δ¹³C-DIC, are plotted against conductivity (Figure 4a–c). As in surficial estuaries, carbon distribution in the beach transect can be examined considering the theoretical conservative mixing (or theoretical dilution line, TDL), where salinity (or conductivity) is used as the conservative property and the composition of the two mixing end-members is assumed to be constant (Boyle et al. 1974). At the inland end of the groundwater transect, i.e. samples collected at M₁ and M₂, DOC was > 1.6 mmol/L (maximum 5.5 mmol/L; Figure 4a), which is an order of magnitude greater than the fresh groundwater end-member. DOC concentrations generally decreased toward the discharge, although the distribution was patchy with some hotspots. For example, high concentrations of DOC were measured in the mid-region of the transect between M₃ and M₄ at 1 m below surface (DOC ~1.6 to 2.5 mmol/L; Figure 4a). In the intertidal zone (i.e. M₃ to M₇), DOC concentrations (~0.2 to 1.2 mmol/L) were also elevated relative to the marine end-member. There was clearly no relationship between DOC concentration and measured conductivity, and most of the samples were higher than the linear mixing model would predict (i.e. the TDL between the two end-members; see Figure 5a). For DIC, most of the measured concentrations were higher than the end-member concentrations (Figure 4b), with peaks in the first 30 cm below the beach surface in the intertidal zone (DIC > 10 mmol/L, Figure 4b). This region corresponded to the position of the surficial saltwater lens. Farther below the surface, concentrations decreased drastically to minimal values (DIC ~0–0.9 mmol/L), which were similar to those of the fresh groundwater end-member. The beach groundwater δ¹³C-DIC was greatly depleted relative to end-members, although values were highly variable through the transect (Figure 4c). δ¹³C-DIC ranged from –34‰ in the upland end of the groundwater transect to values > –5‰ in the top of the intertidal zone. In the bottom part of the transect, at depths deeper than 80 cm below the surface, δ¹³C-DIC was around –25‰. Data obtained in 2012 and 2013 are consistent, with similar values through the beach transect. Because the DIC concentrations of the two end-members differ substantially, the mixing line is curvilinear (details of the mixing equation are given in the caption of Figure 4). In this system, δ¹³C-DIC values are consistently more negative than would be expected with conservative mixing.

Figure 4. Mixing plot of beach groundwater dissolved organic carbon (DOC, a), dissolved inorganic carbon (DIC, b), and isotopic sigmature of dissolved inorganic carbon (δ¹³C-DIC, c) relative to beach groundwater conductivity for 2012 (dots) and 2013 (triangles). The value ranges for the two end-members (i.e., fresh groundwater and marine) and the theoretical dilution line (TDL; black line) are also reported. Concerning the isotopic signature of dissolved inorganic carbon (δ¹³C-DIC): because DIC concentrations of both end-members differ substantially, the conservative mixing between fresh groundwater and seawater is not linear. It follows a curvilinear mixing relationship based on the following equation:

where S is the sample salinity, is the DIC concentration and stable isotope composition of the freshwater end-member, and is the DIC concentration and stable isotope composition for the marine end-member.

Figure 5. Schematic diagram of the hydrogeologic flow at the Martinique Beach site, indicating the potential end-member contribution to volumetric (in bold, unit: m³/s) SGD. Chemical fluxes of dissolved organic carbon (DOC) and dissolved inorganic carbon (DIC) are reported in kg/d. Diagram is not to scale. The tidal range is reported. The dashed line indicates the position of the fresh/saltwater lens. The dotted line is the position of the water table level. Weight of arrows is representative of relative flow strength.

Discussion

A number of techniques have been developed for quantifying SGD (Burnett et al. 2006; Moore et al. 2006; Martin et al. 2007). Here, SGD fluxes have been calculated in two ways: (1) the hydrological approach based on Darcy’s law, which accounts only for the freshwater discharge component in the Permian aquifer and in the unconfined shallow beach aquifer; and (2) using seepage meters, which measure total SGD fluxes at the discharge face of the beach. In general, direct measurements using seepage meters are higher than estimates based on water balance and Darcy’s law calculations. Most of this difference has been attributed to the integration of recirculating seawater. Our hydrological analysis led to a freshwater SGD estimate, Q^inland, of 0.021 m³/s across the 1200 m of sandy beach shoreline of Martinique Bay. This rate assumes a uniform hydraulic conductivity (K) and an isotropic shallow aquifer. Based on the same approach, Q_beach was estimated at 0.029 m³/s. These results agree quite well with the fresh groundwater flux estimate based on the mass-balance approach in the same area (i.e. 0.020 m³/s; Madelin’Eau 2004). The mass balance was based on an average recharge rate to the water table, assuming 36% of the annual average precipitation. Thus, using both Darcy’s law and a mass-balance approach, the fresh groundwater fluxes to the shoreline are equivalent.

Groundwater flux was also estimated by direct measurements during the spring tide to capture maximum outflow from the beach. Because of the microtidal conditions, the effect of tide on SGD would be expected to be minor. The volumetric fluxes obtained by seepage meters at the Martinique Beach site are the highest fluxes measured over the sampling period (Table 1). High seepage rates have also been recorded in several other microtidal sites, including western Australia (Loveless and Odham 2010), the Mediterranean Sea (Weinstein et al. 2011), Cape Cod, Massachusetts (Mulligan and Charette 2006) and Turkey Point, Florida (Santos et al. 2009). At the Martinique Beach site, freshwater SGD estimated using Darcy’s law represents ~25 to 35% of total volumetric additions, indicating that continental drivers (likely water table level) may control SGD at this site. During seepage meter measurements, conductivity was measured in water that had seeped through the seabed. In addition, ambient coastal seawater was collected at the start and end of the deployment. The goal was to capture freshening events and to estimate the volume of seeped fresh groundwater. Conductivity measurements of water collected within seepage bags showed essentially the same values as those from ambient seawater samples. This is not surprising considering the time needed for the exchange of seawater in the seepage meter: with an average Q_seep of 7.5 cm/d and assuming a seepage volume of ~65 L (0.24 m × 0.27 m² = 0.065 m³), it would take more than 3 days. This flushing time seems consistent with our observations and can explain why we observed no change in the salinity. Thus, discrimination of the fresh and marine end-members was not possible in these samples.

Figure 5 proposes a schematic diagram of the hydrologic flow at Martinique Beach, that includes fresh groundwater flows estimated using Darcy’s law and marine SGD derived from tidal pumping and recirculating seawater through the 35-m transect. To compare our SGD estimates with other studies, the volumetric fluxes were divided by the length of shoreline (1200 m). This leads to shore-normal estimates ranging from 1.5 m³/m/d (Q_inland) to 2.1 m³/m/d (Q_beach) for fresh SGD and 6.1 m³/m/d (Q_seep) for total SGD. These shoreline fluxes show magnitudes similar to those at other sites, such as Waquoit Bay (Cape Cod, Massachusetts, 4–5.6 m³/m/d, Mulligan and Charette 2006), the Gulf of Thailand (3.7–7.8 m³/m/d, Burnett et al. 2007), western Australia (2–30 m³/m/d, Loveless and Oldham 2010; and 2.5–8 m³/m/d, Smith and Nield 2003), Shelter Island (Long Island, New York, 0.4–17 m³/m/d, Burnett et al. 2006) and Turkey Point (Florida, 0.7–22 m³/m/d, Santos et al. 2009). At these locations, however, fresh groundwater discharges are a minor component, accounting for only ~5% of total SGD fluxes.

The site-specific evaluation of SGD provides an opportunity to better assess steep spatial gradients in biogeochemical reactions and the resulting impacts SGD may have on biogeochemical fluxes to the coastal embayment. SGD in the unconfined aquifer of Martinique Beach comprises a combination of fresh groundwater and recirculating seawater; thus, the sources of carbon could be diverse. Carbon in the groundwater in Martinique Beach can be classified by its behaviour relative to theoretical conservative mixing between seawater carbon concentrations as one end-member and inland carbon concentrations from municipal and private wells to represent the other end-member. The measured concentrations of DOC and DIC (Figure 4a and b) revealed that the gradient in groundwater-borne carbon from land to ocean is not a simple dilution between the two end-members. In other words, both DOC and DIC showed clear nonconservative patterns in Martinique Beach, with significant internal inputs. These nonconservative behaviours along the groundwater flowpath are important in the interpretation of the ecological impact of SGD (Johannes 1980; Valiela et al. 1980), and the calculation of chemical fluxes using samples from the inland wells probably results in significant deviations from the true chemical flux. Local DOC production has been recorded at other sites, including Turkey Point Beach (Gulf of Mexico, Florida; Santos et al. 2009), Hampyeong Bay (Korea; T. Kim et al. 2012) and North Inlet (South Carolina; Goñi and Gardner 2003). Data on the transit of groundwater-borne DIC are scarce in the literature. Caï et al. (2003) demonstrated that a two end-members mixing model could not explain the high DIC values measured in North Inlet groundwater. More recently, Maher et al. (2013) also reported higher DIC along groundwater transit in an Australian mangrove tidal creek. Other studies have shown that tidal sands are important biogeochemical reactors within which organic matter is mineralized, producing metabolites and DIC (Boudreau et al. 2001; Anschutz et al. 2009; Moore et al. 2011; Charbonnier et al. 2013).

DOC was produced at the inland end of the transect, at the same location where organic-rich paleosol was found (Chaillou et al. 2014). Forest soils and the stumps of trees killed by seawater are commonly found on the US and Canadian Atlantic coast, on the shores of salt marshes, and on the seaward side of eroding barrier beaches (Grant 1970; Dredge et al. 1992; Dubois 1992; Gehrels 1994). The rapid rise in sea level along the coasts of Atlantic Canada during the middle-to-late Holocene buried terrestrial systems with tidal sediments, since sedimentation kept pace with the rising level of high tide. In Martinique Beach, this old enriched horizon of particulate organic carbon (POC) is characterized by a depleted ¹³C signature (~ –33‰; M. Couturier, unpub. data) and high carbon to nitrogen (C/N) ratio (> 25; M. Couturier, unpub. data). DOC is a key intermediate in the overall process of POC degradation (Kristensen and Hansen 1995; Robador et al. 2010; Komada et al. 2012), but little is known about the mechanism of porewater DOC accumulation in suboxic–anoxic sediment. Here, it is still unclear whether the DOC originated from hydrolysis of old POC from the soil horizon or if it was from other sources, such as marine particles or algal detritus. Isotopic investigations should provide further insight into the origin and reactivity of porewater DOC. Whatever its origin, DOC produced locally was transported seaward by groundwater flow. The vertical decline in oxygen saturation observed below the saltwater lens, the maximum residence time in the beach transect of ~32 days in the beach system and the high DOC concentrations were deemed suitable for biogeochemical reactions to occur. The high DIC concentrations measured in the discharge zone combined with high dissolved Fe concentrations (and other reduced compounds such as NH₄⁺ and Mn²⁺; M. Couturier, unpub. data) as well as the characteristic H₂S odour during sample collection previously reported by Chaillou et al. (2014) support the idea that suboxic to anoxic mineralization processes are likely to occur along the flowpath.

Based on total alkalinity and pH measurements, Chaillou et al. (2014) discussed the different ways that DIC may be produced in Martinique Beach. The introduction of fresh meteoric groundwater and the subsequent dissolution of carbonate minerals in addition to the chemical and biological processes of respiration may lead to local DIC production, as observed in other beach systems (Caï et al. 2003; Charbonnier et al. 2013 Maher et al. 2013). δ¹³C-DIC was depleted in the beach groundwater (–34 to –25‰) when compared to the marine and fresh groundwater end-members (–5.65 and +3.05‰, respectively). Since carbonate minerals tend to have δ¹³C values of ~2‰ (Weber and Woodhead 1969; Eadie and Jeffrey 1973), these depleted values may indicate an organic carbon source and may support the assumption that suboxic respiration, in the form of sulfate reduction, was the main DIC source. As seen in Figure 4c, the δ¹³C-DIC data at Martinique Beach largely confirm the nonconservative behaviour of DIC in this subterranean estuary. This pattern is consistent with the earlier hypothesis that intense mineralization in this system may increase DIC and explain DOC loss along the groundwater flow path. Given that the excess DIC should be isotopically similar to the organic matter from which it is derived (which could have various sources, including buried organic-rich horizons, with δ¹³C values between –32 and –34‰), this should result in the addition of ¹³C- depleted DIC.

The ratio of organic carbon to inorganic carbon exported to the coastal ocean could have a significant effect on the role as “source” or “sink” of atmospheric CO₂. Robust estimates of DIC and DOC fluxes are needed on a site-specific scale to obtain accurate regional and global estimations. In nonconservative systems, however, the determination of appropriate DOC and DIC concentration end-members for flux calculations is not straightforward (Beck et al. 2007, 2015; Dorsett et al. 2011; Cyronack et al. 2013). The fresh SGD-derived fluxes of carbon are calculated as the product of the concentration of the fresh groundwater end-member and the volume of fresh groundwater discharge (Q_inland). Assuming that the well closest to the site of discharge (i.e. PM well) represents the best estimate of the nearshore groundwater flux, the carbon load exported by fresh groundwater was 9 kg/d of DIC and 7 kg/d of DOC. In this calculation, biogeochemical transformation during transit from the aquifer to the ocean is completely neglected. Using the fresh groundwater flow estimated in the beach (Q_beach = 0.029 m³/s) and the mean DOC and DIC concentrations measured in the beach’s discharge zone (i.e. mean concentrations in the first 30 cm of M₇), the carbon load was 147 and 27 kg/d of DIC and DOC, respectively. These latter fluxes are 4 to 16 times higher than carbon fluxes derived from nearshore fresh groundwater, and highlight the importance of considering the geochemical behaviour of carbon along the flow path.

For estimating the total SGD carbon fluxes, mean carbon concentrations in the saltwater lens and surficial sediments were used as the recirculating seawater carbon end-member. Total DIC and DOC fluxes exported by total SGD (Q_seep = 0.085 m³/s) could be around 1242 and 147 kg/d, respectively. Fresh groundwater discharges in the beach are an important pathway of carbon and account for 12% of the total DIC and 20% of the total DOC exported to seawater. In contrast, groundwater-borne carbon accounted for less than 1 and 5% of total DIC and DOC fluxes, respectively. Fresh groundwater here is clearly not the source of carbon, but it is a pathway for DIC and DOC produced in the beach to reach the coastal ocean. The chemical SGD fluxes are reported in the schematic hydrogeological flow model in Figure 5.

Only a few papers have considered carbon fluxes in SGD, particularly as it transits through a sandy beach. In Turkey Point (Florida), Santos et al. (2009) reported SGD-derived DOC fluxes 2 times smaller than our values for similar volumetric SGD estimates. In the Okatee Estuary (South Carolina), Moore et al. (2006) measured high inputs of C – 1440 kg DOC/d and 16 800 kg DIC/d – based on flows greater than 1 m³/s. More recently, Liu et al. (2011) measured SGD contributions of DIC that were a hundred times greater than those in Martinique Bay, but their study site was located in a carbonate system on the northern South China Sea, where summer monsoons cause high fresh groundwater flows (~ 2.5 × 10³ m³/s). Thus, it is reasonable to assume that sandy-beach SGD-derived fluxes contribute to regional and global carbon budgets. Once reaching coastal waters, carbon could be consumed by biological production; thus, the net effect of SGD-driven DIC and DOC fluxes would be modulated. This modulation will depend on the origin and reactivity of the exported carbon. Further investigations are needed to determine the source of carbon in Martinique Bay, as this may clarify whether carbon in SGD is from recycled marine matter or degraded terrestrial organic carbon.

Conclusion

In the studied zone, Darcy estimates led to an SGD flow that ranged from 0.021 m³/s (or 1.6 m³/m/d) at the shoreline to 0.029 m³/s (or 3.8 m³/m/d) at the beach face. Direct measurements by seepage meters revealed higher SGD flows, with a mean value of 0.085 m³/s (or 6.1 m³/m/d). These fluxes were consistent with each other and with previous mass-balance estimates based on an average recharge rate to the water table and average annual precipitation. Based on total SGD volumetric fluxes, the total DIC and DOC fluxes exported could be around 1242 and 147 kg/d, respectively. Because the unconfined aquifer of Martinique Beach is a biogeochemically reactive zone, carbon is mineralized along the flow path. This explains why groundwater-borne carbon accounted for less than 1 and 5% of total DIC and DOC fluxes, respectively. Fresh groundwater here is clearly not the source of carbon, but it is a pathway for DIC and DOC produced in the beach to reach the coastal ocean. Our results demonstrate that sandy beaches are an important source of organic and inorganic carbon to the coastal ocean. However, to better quantify SGD-derived fluxes, the biogeochemical complexity of sandy beaches must be closely examined over larger areas and greater time scales. In northern regions where snow and ice cover beaches for a few months of the year, assessments of seasonal changes in the beach flow regime and the production and release of DIC and DOC are critical to determine the impact that northern systems have on regional and global carbon budgets. This site-specific study is the first attempt to estimate volumetric and chemical groundwater fluxes to coastal Canadian waters. It highlights the importance of combining knowledge of marine chemistry and hydrogeology in formulating an adequate management approach for coastal aquifers and adjacent coastal ecosystems.

Acknowledgements

The authors thank Gwendoline Tommi-Morin and Mathilde Gauthier-Pinsonneault for their assistance in the field and for DIC analysis; Tarik Toubal for DGPS measures; Marie-Andrée Roy for producing Figure 1; and Claude and Kathia Bourque for allowing us access to their beach and for the welcome they gave us. Marie Larocque (UQAM) provided valuable input throughout the study, and Laure Devine revised the English phrasing. I. Santos and two anonymous referees are gratefully acknowledged for their thoughtful and constructive comments on the original manuscript. This research was supported by a grant from the Natural Sciences and Engineering Research Council of Canada, the Canada Research Chair Program, and the Université du Québec à Rimouski (UQAR).