Improved Approach for the Investigation of Submarine Groundwater Discharge by Means of Radon Mapping and Radon Mass Balancing

Abstract

The paper presents an improved approach for investigating submarine groundwater discharge (SGD) based on radon mapping and radon mass balancing in the coastal sea. While the use of radon as an environmental tracer in SGD studies is well-established, we identified based on our longstanding experience six methodical shortcomings of the conventional approach and suggest corresponding developments. The shortcomings include: (1 and 2) inadequate consideration of both detection equipment response delay and influence of tidal stage; (3 and 4) incorrect quantification of radon losses, due to offshore mixing and degassing resulting in a potentially incorrect radon mass balance; (5) inaccurate determination of the terrestrial groundwater endmember, due to inhomogeneous radon distribution in the coastal aquifer; and (6) difficulties in distinguishing between discharged fresh groundwater and recirculated seawater. The improved approach is practically demonstrated in a step by step manner in a large-scale field study, which was carried out in False Bay (South Africa) and which consisted of two parts, namely (i) qualitative SGD localization along the entire False Bay coastline based on coastal radon distribution patterns and (ii) quantitative SGD investigation within a defined coastal area of interest (AOI) based on a radon mass balance (RMB). The plausibility of the AOI related results was evaluated by a hydrogeological model, used for qualitative SGD localization, and a hydrological model, applied for estimating groundwater recharge within the AOI catchment.

1. Introduction

The increasing urbanization of worldwide coastlines results in a general increase in environmental stress on coastal ecosystems [1,2]. A key aspect regarding the retention of the related ecosystem services is adequate water resource management in heavily populated coastal zones. That in turn requires a profound understanding of all associated matter budgets and related processes.

Submarine groundwater discharge into the coastal sea (SGD) is a key process in this regard [3,4,5]. Its investigation is of major interest (i) because of the detrimental impact of nutrient- or contaminant-laden groundwater on the coastal sea and (ii) because of the loss of freshwater to the ocean. While the first aspect is of relevance along all coastlines that are impacted by intense urbanization [6,7,8,9], the latter is of particular interest in coastal areas with (seasonally) limited freshwater availability.

Besides air-borne thermal imaging [10,11] and hydrogeological modelling [12,13,14,15], most SGD related studies rely on naturally occurring “environmental” tracers [16,17,18,19,20]. One of the most suitable tracers for the purpose is the radionuclide radon (²²²Rn; half-life 3.8 d). While radon is constantly produced in every aquifer matrix by radioactive decay of its parent nuclide ²²⁶Ra, its production in surface waters is almost negligible. As a consequence the radon activity concentration in groundwater is generally three to five orders of magnitude higher than in surface waters, which allows (i) the use of radon distribution patterns in the coastal sea for the localization of SGD zones and (ii) the use of radon mass balances within defined SGD affected stretches of coastline for the quantification of water (and associated matter) fluxes.

Hence, radon mapping in the coastal sea has become a well-established method for SGD investigation. However, there are a few limitations to this conventional approach, which detrimentally impact both the quality of the raw field data and the reliability of the radon mass balance that is based on this data. Regarding (i) field data collection and (ii) mass balancing six shortcomings (1–6) can be named.

Field data that is collected with the aim of locating SGD zones includes radon distribution patterns in the coastal sea. Onsite radon mapping is generally done from a slowly cruising boat. The (setup specific) response delay of the employed radon detection equipment requires a time correction (time shift) of the recorded radon time series. However, in a considerable number of studies this response delay is either not taken into consideration at all or only allowed for by associating each recorded radon concentration with the geographical position that was passed a certain time before (e.g., 25 min) (shortcoming 1). Furthermore, the effect of the tidal stage, which might change during the radon survey, is often not quantitatively considered in spite of its potentially substantial impact on the recorded data sets (shortcoming 2).

Radon mass balancing for flux quantification within a previously localized coastal SGD zone relies on the precise determination of the mass balance endmembers, i.e., on radon concentrations detected in (i) the terrestrial groundwater, (ii) the seafloor sediment pore water, (iii) the coastal seawater and (iv) the offshore seawater. In the radon mass balance (RMB) all measurable radon inputs (i.e., diffusion from the seafloor sediment and production within the water column) and radon losses (i.e., decay, offshore mixing, and degassing) within the investigated coastal patch are quantified and put in a box model. The excess inventory that is needed to balance the inputs against the losses is then attributed to radon input associated with SGD [16].

While diffusive input, in-situ radon production and radon decay can easily be quantified, the quantification of the offshore mixing loss is more involved. In coastal lagoons or estuaries it can be calculated based on the tidal prism, i.e., the difference in the water volumes of lagoon or estuary between mean high tide and mean low tide [20]. However, this approach does not apply for open shorelines. Here the mixing loss is usually deducted from the maximum decrease of the radon concentration in the coastal sea at a fixed location close to the shoreline recorded during high-tide (i.e., during the expected maximum mixing loss) [21,22]. Still, this approach relies on the (potentially incorrect) assumption of zero SGD during high tide (shortcoming 3). Besides, this approach is only applicable if a long-term single spot time series measurement (over at least two tidal cycles) in the coastal sea is practically feasible at all.

Quantification of the degassing loss is not trivial either. It is mostly quantified as a function of the wind speed recorded during the sampling campaign [23]. However, this approach does not account for any storm events prior to the campaign even though their aftermath may also have a significant impact on the mapped radon inventory (shortcoming 4).

Another shortcoming of the standard RMB approach is the general statistical uncertainty associated with the endmembers, which is often not adequately allowed for [24] (shortcoming 5). The groundwater, seawater and sediment pore water endmembers of the RMB calculations are in most cases deduced from the mean values of the associated radon concentrations detected on site. However, in particular in groundwater and sediment pore water the mean values do not account for the non-Gaussian parameter distribution that is typical for natural systems.

A further limitation that can be named for the conventional RMB approach concerns the fact that two forms of SGD can be distinguished: (i) Pure freshwater discharge (fresh SGD—“FSGD) from the terrestrial aquifer driven by a positive hydraulic gradient and (ii) discharge of seawater that is periodically pumped in and out of the coastal sediment driven mainly by the tides (recirculated SGD—“RSGD) [25]. The RMB alone does not allow a separation of these two SGD components (shortcoming 6). Still, their distinction is essential because the respective chemical compositions and thus the individual effects on ecosystems are different [26,27].

In this paper we suggest and present an improved approach for using radon as an environmental tracer for SGD assessment addressing the named six methodical limitations. Improvement of the coastal radon mapping data is achieved by inverse modelling of the recorded time series, allowing quantitative correction for (1) the response delay of the equipment and (2) the impact of the tides [28]. The suggested improvements regarding the RMB calculations include the quantification of both (3) offshore radon mixing loss based on seawater age estimates derived from radium (²²⁴Ra/²²³Ra) concentration ratios and (4) radon degassing loss by including (potential) degassing prior to the actual sampling campaign. Moreover, we demonstrate how to increase the statistical reliability of the RMB endmembers by (5) conducting stochastic simulations (Monte Carlo) that are based on probability density functions. Finally (6) we show how to independently estimate the shares of FSGD and RSGD based on a salt mass balance. For demonstration the steps of the improved approach were consecutively applied at an exemplary large scale study site in False Bay, South Africa. The gained results were evaluated by independent hydrological and hydrogeological models.

2. Materials and Methods

2.1. Study Area

False Bay, Western Cape Province, South Africa, with a coastline of 110 km was chosen as study area for two reasons: (1) The greater False Bay area suffers from both water shortages during the dry season and (2) signs of contamination of the coastal sea originating from untreated sewage disposal that is often transferred to the sea as storm water runoff from informal settlements.

The overall False Bay catchment (ca. 1400 km²) (Figure 1) shows an altitude ranging from the sea level up to 1511 m. The climate is Mediterranean with an annual average temperature of 16.3 °C, monthly mean summer temperatures of 21 °C, monthly mean winter temperatures of 13 °C [29] and annual precipitation varying between 400 and 800 mm (for Cape Town; [30]). The amount and spatial distribution of rainfall and temperature changes significantly within the catchment, due to topographic effects; precipitation may double in the mountainous parts of the catchment.

The study area belongs to the Cape Fold Belt and its lithology is characterized by exposed bedrocks along the Cape Peninsula in the west, as well as along the eastern flank of the bay that mainly consists of quartzitic sandstone, siltstone and shale of the Table Mountain Group (TMG; Ordovician to Silurian), as well as underlying greywacke, phyllite and quartzitic sandstone of the Malmesbury Group (Neoproterozoic) [31]. The coastal plain of the Cape Flats in the north is made up by unconsolidated sediments of the Sandveld Group (Tertiary to Quaternary) [30]. The mountainous hinterland is widely associated with granitic intrusions (Neoproterozoic). Several streams, e.g., Lourens River and Sir Lowry’s Pass River, are discharging attached sub-catchments. On the Cape Flats, small groundwater fed coastal lakes or “vleis” (e.g., Sandvlei, Zeekoevlei, Rondevlei) with partly poor water quality exist [32].

A large part (40%) of the False Bay catchment is covered with Fynbos vegetation, the natural shrubland and heathland. Agricultural uses cover 29% of the area; the urbanized area encompasses 24%. The population of the wider Cape Town metropolitan area was estimated at around 3.7 million inhabitants. The land use pattern in the sub-catchments of Sir Lowry’s Pass River and Lourens River, which discharge into the coastal area of interest (AOI”, cf. Section 2.2), is comparable to the overall land uses of the False Bay catchment.

False Bay itself covers a water area of ca. 900 km² with an average water depth of 40 m and a maximum depth of 100 m at the bay’s outlet. The tides have a semidiurnal character with a tidal range of up to 2 m at spring tide.

2.2. Methods

We first conducted a radon mapping along the complete False Bay coastline for approximate localization of SGD zones (110 km; four days). The results indicated major SGD occurrence in the north-eastern part of the bay along a ca. 2.5 km stretch of coastline close to Gordon’s Bay and south of Strand (Figure 1). Subsequently a more detailed survey (2.5 km; one day) that aimed at setting up an RMB for this area of interest (AOI) was carried out. The Lourens River catchment and the Sir Lowrys Pass River catchment discharge into the AOI. However, there was no noteworthy surface water discharge from these rivers immediately before and during the field campaign. This observation is reasonable, since our study was carried out at the end of the dry season (March).

2.2.1. Sampling and Analytical Techniques

Radon Mapping and RMB Offshore Endmember Determination

All radon mapping in the coastal sea was done from a vessel cruising at an average speed of about 4 km h⁻¹. GPS was used for recording all cruise tracks. For the four-day False Bay coastal survey the vessel was cruising as close to the shoreline as possible, i.e., in water depths as shallow as possible in order to minimize radon signal dilution. The water depth was constantly monitored by means of an echo sounder. Generally the distance to the shore line was shorter along the eastern and western rocky flanks of the bay (with water depth between about 6 and 10 m) and larger along its northern boundary (the Cape Flats with water depth between about 2 and 4 m). A limiting factor in this regard was the swell along the Caps Flats. However, the overall geographical setting allowed rather constant distances and water depths along the Cape Flats and along the eastern/western flanks of the bay, respectively.

For the one day AOI mapping the vessel track covered the AIO in loops with varying water depth. Again, the water depth was constantly monitored by means of an echo sounder. Additionally, high resolution bathymetry data (South African Council for Geoscience) was used for further data processing.

While cruising, radon was recorded as time series by running two mobile radon-in-air monitors simultaneously (RAD7, Durridge Company, Inc.). Each of the monitors was connected to a radon extraction unit (MiniModule^®, Membrana GmbH) [34] allowing radon extraction from a continuous water pump stream. The detected radon-in-air concentrations were converted into radon-in-water concentrations by allowing for the temperature and salinity dependent radon partitioning coefficient between water and air [35]. For the False Bay coastal survey (and for measuring the RMB offshore endmember 3 km off the shoreline) a 10 min counting cycle of the detectors was applied; a counting cycle of 5 min was chosen for the more detailed AOI mapping.

RMB Endmember Determination for Terrestrial Groundwater and Marine Pore Water

For RMB endmember determination samples were taken from groundwater and marine sediment pore water. Regarding the first, we assumed that beach sediment pore water represents the SGD tracer signature better than groundwater that is taken from wells further in the hinterland. Hence, beach sediment pore water was sampled as terrestrial groundwater endmember. The samples were taken adjacent to the AOI shoreline from depths between 0.5 m and 2 m by means of a piezometric push-point sampler. Marine pore water samples were taken with the help of scuba divers. Lances that were connected to syringes were inserted into the seafloor sediment down to 30 cm. Radon concentrations in all groundwater and sediment pore water samples were measured with a mobile radon-in-air monitor (AlphaGuard PQ 2000; Saphymo) after degassing the radon from the water samples as described in Reference [36]. The AlphaGuard detector works in a different way than the RAD7 utilized on the coastal sea. While the RAD7 counts ²¹⁸Po decays with a solid-state ion-implanted planar silicon alpha detector, the AlphaGuard counts ²²²Rn decays by means of a pulse ionization chamber. The most important difference between the two monitors from the practical point of view is that, on the one hand, the RAD7 requires a decay equilibrium between ²¹⁸Po and ²²²Rn, which takes about 15 min, and on the other hand, the AlphaGuard allows immediate data availability without any decay equilibrium related delay. While this delay was not debilitating on the coastal sea, we used the AlphaGuard for measuring groundwater samples because the sampling and measurement time that was available in particular at wells that were located on private land was on many occasions rather limited. The concentration ranges and detection limits of both monitors are comparable and both suitable for the discussed application.

²²³Ra and ²²⁴Ra Detection in Groundwater and Seawater

The ²²³Ra/²²⁴Ra ratio in the seawater was used as an indicator for the apparent age of the coastal water, i.e., the time that any water parcel has spent in the AOI since entering it through the sediment interface. For quantitative determination of radium concentrations in distinct water samples (incl. ²²⁶Ra) six seawater samples of about 60–100 L and two groundwater samples of 2.5 L were collected. Radium was extracted from the water samples by means of MnO₂ coated acrylic fibers as described in Reference [37]. Aditional radium sampling was performed by placing three mesh bags containing MnO₂-fiber within the AOI water body for one day. No absolute radium-in-water concentrations (and hence no ²²⁶Ra activities), but only ²²⁴Ra/²²³Ra ratios resulted from the three mesh bag samples. All MnO₂-fiber samples were measured using a Radium Delayed Coincidence Counting System (RaDeCC) as described in Reference [38]. For the efficiency calibration of the RaDeCC, ²²⁷Ac and ²³²Th standards were used [39,40]. Counting errors were propagated following [41].

2.2.2. Processing of the Radon Mapping Data

All radon time series recorded during the mapping cruises by means of the RAD7 monitor required a response delay correction (cf. shortcoming 1). To meet this requirement a method developed earlier by two of the authors [28] was applied. The method accounts by inverse time series modelling for (i) the water: Air transfer kinetics of radon and (ii) the delayed decay equilibrium between ²²²Rn and its progeny ²¹⁸Po (that is actually being counted by the RAD7 monitor). The data processing approach requires a mass transfer coefficient, which depends on both the detection set-up and water pump stream flow rate. The mass transfer coefficient for the set-up applied in this study had previously been quantified in the laboratory to be 1.3 × 10⁻³ s⁻¹ [28].

In order to achieve comparability of all radon data recorded during both the four-day radon cruise along the False Bay coastline and the one day radon mapping within the AOI all radon time series were corrected for the varying tidal influence (cf. shortcoming 2). The effect of the varying tidal stage on the radon concentration in the coastal seawater was continuously monitored at a fixed location in Gordons Bay harbour.

2.2.3. Radon Mass Balance for the AOI

General RBM Setup

Radon mass balances (RMB) generally assume a steady-state, i.e., a constant radon inventory within the AOI water body (I_Rn). This radon inventory is expressed as Becquerel over a square meter [Bq m⁻²], i.e., the product of radon concentration and water depth. The AOI radon inventory is balanced by the radon sinks and sources mentioned above (cf. sec. 1) [16,42,43]. The sinks include radon decay (F_dec), offshore radon mixing loss (F_mix) and atmospheric radon degassing (F_atm). The sources are (apart from the not directly detectable parameter of interest: The radon input associated with SGD) diffusive radon input from the sea floor sediment (F_dif) and radioactive production within the water column (F_prod). In our case study any other radon input was assumed negligible, since discharge of the two rivers into the AOI was negligible during and prior to the field campaign (<0.01 m³ s⁻¹). After quantifying the directly measurable sources and sinks [Bq m⁻² d⁻¹] the positive radon flux that is required to sustain a balanced AOI radon inventory was attributed to SGD (F_SGD) (Equation (1)). Based on the radon concentration of the SGD endmember (i.e., the radon concentration of beach sediment pore water; Rn_SGD; [Bq/m⁻³]) the SGD related radon flux was finally converted into the actual SGD water flux [m d⁻¹].

$${\mathrm{F}}_{\mathrm{SGD}}\left[{\mathrm{Bq} \mathrm{m}}^{-2}{ \mathrm{d}}^{-1}\right]={\mathrm{F}}_{\mathrm{dec}}+{\mathrm{F}}_{\mathrm{mix}}+{\mathrm{F}}_{\mathrm{atm}}-{\mathrm{F}}_{\mathrm{dif}}-{\mathrm{F}}_{\mathrm{prod}}$$

(1)

Calculation of the Offshore Radon Mixing Loss

The net offshore radon mixing loss (cf. shortcoming 3) was calculated from estimations of the mean water apparent age of the AOI water body derived from the ²²⁴Ra/²²³Ra concentration ratio in the seawater (RaR_SW) and the related SGD endmember (RaR_SGD). The concept of the approach is comparable, but not similar to the one applied by Reference [43]. The calculation of the water apparent age (τ) assumes that there are no other radium inputs to the discharged groundwater while moving offshore with time (cf. [43]). Consequently, after groundwater discharge into the sea, the ratio changes only as a function of time, due to the difference of the two decay constants (λ_Ra223 = 7.01 × 10⁻⁷ s⁻¹ and λ_Ra224 = 2.20 × 10⁻⁶ s⁻¹). The ²²⁴Ra/²²³Ra activity ratio “decays” with a combined half-live of 5.4 days. The time that elapsed, since a parcel of water was last in contact with the aquifer matrix can be thus derived using Equation (2) [44].

$$\mathsf{\tau }=\frac{\ln \left(\frac{{\mathrm{RaR}}_{\mathrm{SW}}}{{\mathrm{RaR}}_{\mathrm{GW}}}\right)}{\left({\mathsf{\lambda }}_{{\mathrm{Ra}}_{224}}-{\mathsf{\lambda }}_{{\mathrm{Ra}}_{223}}\right)}$$

(2)

Based on τ we simply defined the daily net mixing rate (mix) as the reciprocal of the mean water apparent age of the entire AOI water body (Equation (3)). If, for instance, the radium ratios indicated four days as mean apparent age of the AOI water body, i.e., of the groundwater that discharged into the AOI, we concluded that the complete AOI water body is exchanged by offshore mixing every four days resulting in a net mixing rate of mix = 1/4 d⁻¹. The mixing loss (F_mix) is finally calculated based on the assumption that the intensity of the mixing between offshore water and coastal water increases as an exponential function with the distance to the shoreline (Equation (3)).

F_mix [Bq m⁻² d⁻¹] = I_Rn [Bq m^-2](1 − e^−mix) with mix[d⁻¹] = τ [d]⁻¹

(3)

Calculation of the Radon Degassing Loss

Our approach to quantifying radon degassing is based on the empirical model developed by Reference [23], which assesses the loss rate as a function of the actual wind speed [16,22,43]. Still, the approach does not take into account any intense storm events that occurred prior to the actual sampling campaign (cf. shortcoming 4). Since the radon inventory needs time to regenerate after such a degassing event, a time span that starts several days prior to the sampling should be included in the calculations.

Since the effect of a previous degassing event on the detected radon concentrations is the more influential the closer the degassing event has occurred prior to the actual sampling campaign we defined a weighting factor that quantifies the influence of the previous degassing for the time step t [d]. The weighting factor w_t was calculated based on the radon decay constant (λ_Rn = 0.181 d⁻¹) and the offshore water mixing rate given in Equation (3) (Equation (4)).

$${\mathrm{w}}_{\mathbf{t}} [-]={\mathrm{e}}^{\left(-({\mathsf{\lambda }}_{\mathrm{Rn}}+\mathrm{mix}\right)\times \mathrm{t})}$$

(4)

The effective radon degassing that affects the radon inventory during a sampling campaign is then calculated following Equation (5).

$${\mathrm{F}}_{\mathrm{atm}} \left[{\mathrm{Bq} \mathrm{m}}^{-2}{ \mathrm{d}}^{-1}\right]=\frac{{{\displaystyle \sum }}_{\mathrm{t}=0}^{{\mathrm{t}}_{\max }} \left[\mathrm{w} \left(\mathrm{t}\right)\times {\mathrm{F}}_{\mathrm{atm}}\left(\mathrm{t}\right)\right]}{{{\displaystyle \sum }}_{\mathrm{t}=0}^{{\mathrm{t}}_{\max }} \mathrm{w} \left(\mathrm{t}\right)}$$

(5)

Calculation of Radon Decay, Production and Diffusive Flux

Radon production (F_prod) and decay (F_dec) in seawater are solely functions of the concentrations and decay constants of ²²⁶Ra and ²²²Rn, respectively. The diffusive radon flux (F_dif) is mainly driven by the radon concentration gradient between seawater and marine pore water and was calculated from a depth independent model first introduced by Reference [45]. More detailed descriptions on calculations of F_prod, F_dec and F_dif can be found elsewhere [16,22,46].

RMB Endmember Uncertainty Propagation

The RMB approach requires radon endmember concentrations as input parameters. Calculation of the individual endmembers is often based on an only limited number of samples. Considerable uncertainty arises hence in particular in case of spatial heterogeneous radon distributions in groundwater (cf. shortcoming 5). In order to decrease these uncertainties, we conducted stochastic Monte Carlo simulations of the RMB that were based on probability density functions (pdf) for propagating the uncertainties of the input parameters into the final SGD flux uncertainty. In detail, a large number of model realizations (in our case 10,000) were run with varying input parameter settings. The parameter setting for each realization was randomly sampled from the pdf of the respective parameter. Hence, a range of possible results was computed allowing for uncertainty quantification of the model output. While the uncertainties of the on-site radon measurements were assumed to be Gaussian distributed, the RMB input parameters (in particular the radon concentrations detected in the groundwater) were assumed to be lognormal distributed in space as it is typical for natural systems [47].

2.2.4. FSGD/RSGD Differentiation

With the aim to distinguish between the fractions of freshwater (f_FSGD) and recirculated seawater discharge (f_RSGD) within the AOI water body (cf. shortcoming 6) an independent salt mass balance was set up for the AOI. We assumed that FSGD is the only freshwater source to the AOI. Precipitation and surface water runoff, as well as evaporation were assumed to be negligible during the sampling campaign and immediately prior to it. Further, we assumed that RSGD has no effect on the AOI water salinity. So f_FSGD was determined from the mean salinity of the AOI water body (Sal_AOI) and the salinity of the offshore seawater endmember (Sal_SW) (Equation (6)).

$${\mathrm{f}}_{\mathrm{FSGD}}=1-\frac{{\mathrm{Sal}}_{\mathrm{AOI}}}{{\mathrm{Sal}}_{\mathrm{SW}}}$$

(6)

In a subsequent step, we calculated the flux of FSGD (F_FSGD) that is required to sustain f_FSGD assuming that the only process that removes freshwater from the AOI is offshore mixing. Based on the mixing rate and the AOI water volume (V_AOI) the FSGD flux can be estimated based on f_FSGD as follows (Equation (7)).

$${\mathrm{F}}_{\mathrm{FSGD}}=\left(1-{\mathrm{e}}^{-\mathrm{mix}}\right)\times {\mathrm{V}}_{\mathrm{AOI}}\times {\mathrm{f}}_{\mathrm{FSGD}}$$

(7)

2.2.5. Validation of the RMB Results Using Hydrogeological and Hydrological Models

The results of the RMB calculations were validated by means of two model approaches, a hydrogeological numerical groundwater flow model and a hydrological model. Both models were run for a 210 km² model area that included the 9.3 km² AOI catchment (Figure 1).

The groundwater flow model was set up as transient, density-coupled, three-dimensional finite element model with the aim (i) to indicate locations along the AOI coastline were FSGD is likely to occur and (ii) to define the AOI catchment. The model was based on hydrogeological structural data [48] containing information about the general hydrogeological setting, as well as the spatial extend of bedrock units, weathered materials and younger sediments. Additional geological data was provided by the Council of Geoscience, South Africa, and has been used for compiling a differentiated hydrostratigraphy for the model area. The data allowed the flow model parameterization in a horizontal and vertical scale with regard to the actual hydraulic conductivity parameter distribution. All available water level data was applied for a steady-state model calibration and a sufficient model calibration was reached even though existing spatial and temporal information on groundwater levels were limited. The resulting groundwater flow model was transferred into a time-variant flow model for considering probable seasonality of FSGD fluxes. Furthermore sea level variations were implemented by using the local tidal characteristics. The transient model was run for a modelling period of 365 days. Results of individual model run served repetitively as initial conditions for a subsequent simulation run until “stable” inter-annual variabilities of the state variables were reached.

For setup of the hydrological model we applied the WebGIS based software STOFFBILANZ [49,50,51]. The hydrological model aimed at providing estimates of groundwater recharge and surface water runoff rates for the model area and the AOI catchment [51,52,53,54]. The input data include daily rainfall, daily temperature, daily reference evapotranspiration, soil hydraulic properties, land cover and the growing state of crops.

3. Results

3.1. SGD Localization

3.1.1. Radon Mapping in the Coastal Sea

Tide Correction of Radon Time Series

All radon time series that were recorded from the cruising boat both along the False Bay coastline and within the AOI were tide corrected. The corrections were based on simultaneous radon time series that were recorded at a fixed location in the Gordon’s Bay harbour area. The harbour time series revealed a high inversely correlated dependence of the radon concentration on the tidal stage (r = −0.84, R² = 0.71, p < 0.001) (Figure 2A). Based on this data we selected an exponential regression model for tidal correction of the actual radon mapping time series, since it allows for relative rather than for absolute corrections of the individual concentrations, thereby avoiding negative concentrations as result of the correction. As a result tide-corrected radon concentrations (Rn_w; [Bq m⁻³]) were derived from the measured concentrations (Rn_obs; [Bq m⁻³] by making allowance for the associated tidal stage [m] i.e., for the associated deviation from the mean water level (Equation (8)).

The tide data was recorded in Simon’s Town on Cape Peninsula (data derived from the University of Hawaii, Sea Level Center). An exponential regression coefficient of −0.61 produced the best data fit and was thus used for tidal correction of the radon mapping time series (Figure 2B). We assume this relationship to be valid along the entire False Bay coastline. Since the coastal bathymetry doesn’t show any distinctive features, such as isolated steep seafloor slopes or shallow ridges close to the coast, we assume tidal range and tidal onset to be approximately the same all over the complete coastline.

$${\mathrm{Rn}}_{\mathrm{w}}=\left({\mathrm{Rn}}_{\mathrm{obs}}\right)/{\mathrm{e}}^{\left(-0.61\times \mathrm{tide}\right)}$$

(8)

Radon Distribution along the False Bay Coastline

The radon cruise along the entire False Bay coastline revealed low radon concentrations (<20 Bq m⁻³) along the bedrock sections of both the Cape Peninsula and the eastern flank of False Bay and along the western part of the Cape Flats. In contrast, significantly higher concentrations (potentially indicating SGD) were observed along the middle and eastern part of the Cape Flats reaching highest concentrations south of Strand and close to Gordon’s Bay (up to 55 Bq m⁻³). All recorded time series of the four-day mapping cruise were response delay and tide corrected.

Radon Distribution within the AOI

Based on the results of the radon mapping along the entire False Bay coastline an AOI was defined in the Strand/Gordon’s Bay area and a subsequent radon mapping campaign with a higher spatial resolution (5 min counting cycle) was performed there. Figure 3 displays the radon distribution pattern mapped within the AOI. In order to illustrate our approach for stepwise data correction Figure 3A–D show the results of the consecutive steps of data processing.

For the plot in Figure 3A the raw radon time series was not corrected at all but only shifted 2.5 min backwards in order to position each data point to the middle of the cruise track section that was covered during the corresponding 5 min counting cycle. Figure 3B shows the same concentration distribution pattern, but shifted backwards by an additional 25 min in order to roughly allow for the response delay of the equipment (reflecting the often-used conventional approach). Figure 3C shows the radon time series as in Figure 3A, but (in contrast to the approximate correction shown in Figure 3B) quantitatively response delay corrected. The radon anomaly is now located between the anomalies that can be identified in Figure 3A,B, respectively. Furthermore, the radon concentrations increase as a result of the inverse modelling that removes the radon detectors’ signal smoothing effect (as discussed in Reference [28]). Figure 3D shows finally the radon time series both response delay and tide corrected. Since the AOI mapping was conducted mostly at ebb tide, concentrations increase (with varying intensity) further after tide correction. As a result, the spatial radon patterns changed slightly in comparison to Figure 3C, due to the correction of different tidal stages during mapping.

The results of the AOI mapping confirm and improve the findings of the cruise along the entire coastline of the bay that indicated SGD close to Gordon’s Bay. Furthermore, it can be seen that radon concentrations within the AOI are decreasing with increasing water depth, due to radon dilution.

3.1.2. Hydrogeological Modelling

Hydrogeological modelling was applied for validation of the AOI radon survey results. The average hydraulic heads calculated by the model reveal insight into the local groundwater flow pattern (Figure 4). The hydraulic gradient is steepest south and east of the AOI indicating FSGD potential there. In contrast, distances between the hydraulic contour lines are much wider north of the AOI indicating FSGD to be less likely there.

The AOI catchment, i.e., the recharge area of groundwater that discharges into the AOI, was delineated by applying backward particle tracking. The capture zone, of which parts can be seen in Figure 4, has a spatial extent of 9.3 km².

Based on this information the spatial distribution of potential FSGD occurrence along the AOI coastline was simulated. The AOI coastline was divided into ten subsections of equal length. For each subsection the relative FSGD contribution in relation to the overall FSGD along the AOI coastline was calculated. Figure 4 illustrates clearly that the model indicates the main share of FSGD in the southern part of the AOI. One single subsection receives 57% of the total FSGD, whereas each of the subsections in the northern part receives less than 1%. However, input data is scarce. Therefore, the exact localization of FSGD only based on the hydrogeological model data is associated with high uncertainty. Still, the modelled data confirms the radon survey results.

3.2. SGD Quantification within the AOI

3.2.1. Radon Mass Balance

AOI Definition and Radon Inventory Calculation

The extent of the AOI for which the RMB was calculated was derived from the initial radon mapping data. The AOI covers an area of 0.85 km² with 2500 m of coastline. The average water depth of the AOI is 1.4 m; its water volume amounts to approximately 1.19 × 10⁶ m³.

All radon data that resulted from the AOI mapping (5 min counting cycle) was spatially inter-/extrapolated for the extent of the entire AOI area using universal kriging. Subsequently, the radon inventory of the AOI was calculated under the assumption of a vertically well-mixed water body. The mean radon concentration in the entire AOI water body was 41 Bq m⁻³; the mean radon inventory was 57 Bq m⁻². The radon inventory of the whole AOI was derived as a spatially integrated value from this data (ca. 48,800 kBq).

Calculation of AOI Radon Losses and Radon Inputs

Mixing Loss

The mean water apparent age of the AOI water body was determined from the ²²⁴Ra/²²³Ra ratios in groundwater (RaR_GW, n = 2) and AOI water (RaR_SW, n = 9). The mean ratios were found to be RaR_GW = 16.0 ± 3.1 and RaR_SW = 9.7 ± 2.2. The data revealed a mean apparent age of the AOI water body of 3.86 d (25th/75th perc.: 2.4/5.6 d) with an uncertainty of ±2.5 d (considering the uncertainty of the radium endmembers) (Equation (2)). The mean mixing rate derived thereof amounts to 0.26 d⁻¹ (25th/75th perc.: 0.17 and 0.37 d⁻¹). Thus, the median radon loss, due to offshore mixing is 13.1 Bq m⁻² d⁻¹ (25th/75th perc.: 9.0/18.0 Bq m⁻² d⁻¹) (Equation (3)).

Degassing Loss

The effective radon degassing loss was calculated not only based on the wind speed during the day of the AOI cruise, but including the wind speed history for a five day period prior to the AOI sampling campaign. Historic and recent wind speed data was provided in hourly resolution by the South African Weather Service (SAWS) as detected at the observatory in Strand located only 2 km north of the AOI. An average atmospheric radon concentration of 30 Bq m⁻³ was detected during the AOI sampling campaign and assumed similar for the preceding five days as well.

Figure 5 illustrates wind speed, resulting in associated atmospheric radon flux, historic weighting factor and final weighted atmospheric radon loss for the prior five-day period. The weighting factor shows an exponential decrease with increasing temporal distance to the AOI campaign assuming a constant offshore radon mixing rate during this five day period, i.e., no changing influence of the changing wind speed on the mixing rate (Equation (4)). The latter assumption is a simplification, but a sound reconstruction of the mixing rate history is hardly possible. Recorded wind speed and the resulting calculated atmospheric radon flux are inversely correlated during the five-day period.

The atmospheric radon flux, which was finally used for the RMB setup, results from the sum of all products of (historic) atmospheric fluxes and corresponding weighting factors divided by all weighting factors (Equation (5)). The resulting radon loss for the RMB was 32.2 Bq m⁻² d⁻¹ (25th/75th perc.: 25.2/39.2 Bq m⁻² d⁻¹).

Decay Loss

The daily radon decay loss used for the RMB results from the detected radon inventory of the AOI water body (57 Bq m⁻² d⁻¹) and the radon decay constant (2.1 × 10⁻⁶ s⁻¹) and amounts to 9.5 Bq m⁻² d⁻¹.

Diffusive Input from Sediment

An average radon concentration of only 100 Bq m⁻³ (n = 29) was detected in the marine pore water at ca. 30 cm sediment depth. The associated concentration gradient between seawater and marine pore water within the AOI resulted in a diffusive radon flux from the sediment of 0.1 Bq m⁻² d⁻¹.

In order to evaluate the rather low pore water concentrations that were measured on site, marine sediment samples were taken and analysed for ²²⁶Ra activity concentrations using gamma-spectroscopy. The detected mean ²²⁶Ra activity was 5.1 kBq kg⁻¹. With a measured sediment dry density of 1700 kg m⁻³, an assumed emanation coefficient of 0.15 and measured pore space of 0.38 this in-situ radium activity yields theoretically a radon pore water concentration of about 3.4 kBq m⁻³. The substantial difference between this potential radium supported radon concentration and the concentration that was actually detected in the shallow pore water indicates intense flushing of the upper sediment layer, due to tide-driven seawater circulation (tidal pumping). The detected radon concentration in the shallow pore water, if compared to the potential radium supported radon concentration, indicates a mean residence time of infiltrated seawater within the shallow sediment layer (i.e., the time it takes for any water parcel between entering and leaving the sediment) of about 4 h, a value that is consistent with the tidal pumping concept.

Input by Seawater in-situ Radon Production

The average ²²⁶Ra concentration in the coastal seawater was 1.8 Bq m⁻³ (n = 6). The daily radon production within the AOI water body that is supported by this activity amounts to 0.42 Bq m⁻² d⁻¹.

Terrestrial Radon Endmember in Beach Sediment Pore Water

The SGD radon endmember was derived from beach sediment pore water samples (n = 17) taken in the vicinity to the AOI coastline. The samples revealed a salinity range from 0 to 34‰ (i.e., from fresh groundwater to seawater) and a radon concentration range from about 1000 Bq m⁻³ to about 50,000 Bq m⁻³. Samples with salinities above 5‰ were immediately excluded from considerations regarding the SGD endmember definition because they were presumably taken too close to the seawater/groundwater interface and were hence mainly composed of infiltrated seawater. 13 selected samples were finally used for SGD endmember definition. However, instead of simply using their mean value, we rather produced an empirical probability density function (pdf) from the data set, including the associated analytical uncertainty. Under the assumption of a lognormal distribution of radon in groundwater we fitted a lognormal pdf to the empirical pdf. The terrestrial radon endmember derived from the fitted pdf was 4800 Bq m⁻³ (25th/75th perc.: 2100 Bq m⁻³/11,500 Bq m⁻³).

Radon and Water Flux Associated with SGD

The SGD bound radon flux was calculated by balancing all determined radon losses (offshore mixing, degassing, decay) and sources (diffusion, in-situ production). The resulting radon excess within the AOI water body, i.e., the required radon input that had to be associated with SGD amounted to F_SGD = 54.3 Bq m⁻² d⁻¹ (25th/75th perc.: 44.2/63.3 Bq m⁻² d⁻¹). On the loss side, degassing was the most important sink followed by offshore mixing and decay. On the source side, both diffusion and in-situ production are virtually negligible compared to the radon flux associated with SGD.

The SGD radon flux was converted into a total SGD water flux based on the terrestrial radon endmember. The resulting SGD estimate for the AOI is 11.3 mm d⁻¹ (25th/75th perc.: 3.8/24.8 mm d⁻¹) or 9615 m³ d⁻¹ (25th/75th perc.: 3300/21,300 m³ d⁻¹) or 3.8 m³ d⁻¹ per meter coastline (25th/75th perc.: 1.3/8.5 m³ d⁻¹ per m coastline).

3.2.2. Differentiation between FSGD and RSGD

During the AOI campaign we observed a mean salinity of 34.5‰ within the AOI water body and salinity of 34.8‰ in the offshore seawater detected ca. 3 km off the coast. Assuming a freshwater salinity of 0‰ this yields a freshwater fraction in the AOI water body of f_FSGD = 0.86% (which equals a freshwater volume of 10,259 m³). Considering the determined mixing rate of mix = 0.26 d⁻¹ and an AOI Volume of V_AOI = 1.19 × 10⁶ m³ the daily freshwater flux that is required to sustain the freshwater fraction was calculated to be F_FSGD = 2339 m³ d⁻¹ (25th/75th perc.: 1700/3400 m³ d⁻¹).

F_FSGD = V_AOI (1 − e^−mix),

(9)

or 0.9 m³ d⁻¹ per m coastline (25th/75th perc.: 0.6/1.3 m³ d⁻¹ m⁻¹). Subsequently, RSGD was estimated by subtracting the salt balance based FSGD flux from the RMB based total SGD flux (9615 m³ d⁻¹) resulting in an RSGD of 7276 m³ d⁻¹ (25th/75th perc.: 1100/19,500 m³ d⁻¹) or 2.9 m³ d⁻¹ per m coastline (25th/75th perc.: 0.5/7.9 m³ d⁻¹ m⁻¹).

In order to evaluate the plausibility of the calculated FSGD flux rate the annual groundwater recharge was modelled for the AOI catchment resulting in a rate of 218 mm a⁻¹, which is equivalent to 5555 m³ d⁻¹. Still, significant seasonal variations are present. For instance 54% of the total precipitation and 93% of groundwater recharge occur during the winter (June to August). Our study was carried out at the end of the dry season (March).

4. Discussion

4.1. Appropriateness of the Introduced Approach

Several improvements were applied alternatively to the established method for SGD investigation using radon as an environmental tracer. Modifications of the existing approach concerned both field data processing and RMB setup. Regarding field data processing the allowance for (1) an equipment specific response delay and for (2) a time-variant tidal impact on the radon time series resulted in more precise localization and quantification of radon anomalies within the coastal sea. Modification of the RMB setup resulted (3) in a more reliable quantification of offshore mixing loss and (4) degassing loss, as well as (5) in a more accurate terrestrial SGD endmember. Furthermore, it was shown that (6) approximate separation of FSGD and RSGD is possible based on an RMB independent salt mass balance.

4.2. General SGD Occurance in False Bay

The radon mapping results are generally consistent with the indications of the hydrogeological model. Both approaches point to major SGD occurrence in the southern part of the AOI. The two patterns do not perfectly match, though. This is most properly due to the limited data availability for the groundwater flow model.

The considerable SGD occurrence in the Gordon’s Bay area is a consequence of both relatively high groundwater recharge in the mountainous hinterland and the general hydrogeological setting that results in steep hydraulic gradient south-east of the AOI. Furthermore the outcropping Malmesbury Shale acts as aquitard and forces the overlying groundwater to discharge into the sea.

The western section of the Cape Flats shows a much lower SGD rate then the middle and in particular eastern part. This is probably due to groundwater withdrawal by agricultural activities (Philippi region; cf. Figure 6), and due to groundwater discharge into the coastal lakes. Strong evaporation in these lakes impedes the seaward groundwater flow.

SGD along the rocky western and eastern flanks of False Bay is very low due to the low hydraulic conductivity of the rocks and the resulting low groundwater recharge rate. Moreover, the small catchment areas along the two flanks limit groundwater recharge.

4.3. SGD Quantification for the AOI

The estimation of a total SGD flux of about 9600 m³ d⁻¹ (or 3.8 m³ d⁻¹ per m coastline) within the AOI is based on the RMB results and the subsequent conversion of the radon flux into a groundwater flux. The uncertainty of the SGD flux estimate results primarily from the terrestrial radon endmember uncertainty. For the differentiation between FSGD and RSGD based on a salt balance we considered FSGD as the only noteworthy source of freshwater input into the AOI. Since insignificant freshwater input might occur (through the two ephemeral rivers) the estimated FSGD of about 2340 m³ d⁻¹ (or 0.9 m³ d⁻¹ per m coastline) should be seen as an upper limit.

As a plausibility check, we compared the tracer based FSGD estimate (2340 m³ d⁻¹) to the groundwater recharge simulated by the hydrological model within the AOI catchment (5555 m³ d⁻¹). The comparison reveals that the modelled groundwater recharge rate is about twice as high as the tracer based FSGD estimate. While for a hydrological system in interannual steady state groundwater recharge should generally be equivalent to FSGD it has in our case to be kept in mind though, that the results of our tracer study reflect the conditions during the end of the dry season (March). In the False Bay area roughly 90% of groundwater recharge occurs from June to August. A study by Reference [55] revealed seasonal oscillations of SGD lagging three to five months behind peak recharge rates. Furthermore FSGD is presumably reduced by unmonitored groundwater abstraction (e.g., by private households in informal settlements) and groundwater discharge into the coastal lakes, as mentioned above. Hence, the tracer based FSGD estimate is still reasonable.

5. Conclusions

The key findings of the study are:

(1) Raw radon mapping data has been corrected for response delay and tidal effects in a way that improved data quality considerably. The identified SGD location is in good agreement with results from a numerical groundwater flow model;

(2) An alternative approach for calculating the radon offshore mixing loss was presented;

(3) An improved approach for the radon degassing loss, including previous degassing events was presented.

The applied improved approaches for radon parameter estimation and radon flux calculation contribute to the methodological progress in localizing and quantifying SGD by applying the environmental tracer radon and in determining the related uncertainties. The presented approaches are transferable to other groundwater/surface water systems, such as lakes, estuaries or lagoons.