Next Article in Journal
Is Greenhouse Rainwater Harvesting Enough to Satisfy the Water Demand of Indoor Crops? Application to the Bolivian Altiplano
Previous Article in Journal
Predicting Urban Flooding Due to Extreme Precipitation Using a Long Short-Term Memory Neural Network
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Communication

Open-Source Code for Radium-Derived Ocean-Groundwater Modeling: Project Open RaDOM

Natural and Applied Sciences, College of Arts & Sciences, Lynn University, Boca Raton, FL 33431, USA
Submission received: 11 May 2022 / Revised: 4 June 2022 / Accepted: 10 June 2022 / Published: 14 June 2022
(This article belongs to the Section Marine Environment and Hydrology Interactions)

Abstract

:
Radium has been commonly used as a tracer of submarine groundwater discharge to the ocean and embankments, as radium activities are commonly input into box models to calculate a groundwater flux. Similarly, isotopes of radium (Ra224, Ra223, Ra226, Ra228) have been used to calculate water mass ages, which have been used as a proxy for residence times. Less commonly, radium and other tracers have been utilized in mixing models to determine the relative contribution of groundwater to a marine system. In the literature, all of these methods have almost exclusively been solved using analytical methods prone to large errors and other issues. Project Open RaDOM, introduced here, is a collection of open-source R scripts that numerically solve for groundwater flux, residence time, and relative contribution of groundwater to coastal systems. Solving these models numerically allows for over-constrained systems to increase their accuracy and force real solutions. The scripts are written in a way to make them user-friendly, even to scientists unfamiliar with R. This communication includes a description of the scripts in Project Open RaDOM, a discussion of examples in the literature, and case studies of the scripts using previously published data.

1. Introduction

Submarine groundwater discharge (SGD) is the discharge of fresh groundwater, sea water that has circulated through the coastal aquifer, or a mixture of the two water types to the coastal ocean [1]. SGD transports nutrients, gases, and/or pollutants to coastal systems that impact the ecology of those systems [2]. As such, the flux SGD and its associated solutes have been quantified worldwide [3,4]. SGD is identified and quantified in coastal systems thermally, via salinity differences, and by using tracers such as radon and radium [5,6,7,8]. One of the most common methods utilized is the quantification of SGD via radium activities and box models [4,9].
There are four radium isotopes typically used in the study of SGD. These isotopes are Ra224, Ra223, Ra228, and Ra226 with half-lives of 3.6 d, 11.4 d, 5.75 y and 1600 y, respectively [10]. In coastal aquifers, radium is a product of decay from its parent isotopes and is adsorbed onto sediment surfaces. Radium desorbs via ionic exchange in brackish and saline systems, and is therefore an indicator of brackish or saline SGD [11]. Once in the coastal ocean or an embayment, it is removed via radioactive decay or advection offshore [12].
The simplest form of SGD flux models are calculated with a single isotope of radium, assuming that the only source is SGD and that the only sink is advection offshore [13]. These models are calculated using the equation:
S G D   F l u x = ( R a B o x R a O f f ) × V τ × R a G ,
where SGD Flux is the flux in volume/time of groundwater to the coastal ocean, RaBox, RaOff, and RaG are the radium of activities of a given isotope in the coastal box that is receiving SGD, offshore waters, and the coastal groundwater, respectively. V is the volume of the coastal box, and τ is the residence time of the coastal box. The benefits of such box models are that they only require activities of a single radium isotope and are easy to solve analytically. The disadvantages of this model include that when different radium isotopes are used, wildly different SGD fluxes may be computed, leaving the researcher in somewhat of a conundrum (e.g., [14]), and sometimes a residence time can be hard to ascertain (e.g., [15]). Additionally, radium isotope ratios are often used to calculate the residence time, which can have errors of up to 100% when residence times are less than 5 days using the shorter-lived isotopes [16].
Another way to evaluate the importance of SGD to a coastal system is through mixing models. Using radium isotopes to determine SGD contribution was first pioneered by Moore [17] where he evaluated the contribution of two different groundwaters to a continental shelf using the equations:
R a O 1   F O + R a G 1 1 F G 1 + R a G 2 1 F G 2 = R a S 1 ,
R a O 2   F O + R a G 1 2 F G 1 + R a G 2 2 F G 2 = R a S 2 ,  
F O + F G 1 + F G 2 = 1 ,
where FO, FG1, and FG2 are the fractions of offshore water and two different groundwaters to each sampling location on the continental shelf, Ra1 and Ra2 are two radium isotopes observed, RaO, RaG1, and RaG2 are the activities of radium in offshore water and each of the groundwaters, and RaS is the radium activity in the sampling location on the continental shelf. The third equation forces a solution in which the fractions add up to 1. While useful as a tool to compare the relative contribution of two or more groundwaters, when solved analytically these equations require very conservative tracers and therein only the long-lived radium isotopes can typically be used. In addition, the analytical solution is rather large and burdensome to calculate, as shown in Moore [17].
While both these analytical models have served to quantify fluxes and groundwater contributions to marine systems worldwide, with the advent of free to use programming languages and software, improvements can be made upon their basic premises, namely numerical modeling. In contrast to analytical modeling in which equations are solved explicitly for an exact solution, in numerical modeling solutions are estimated based off algorithms. Advantages of numerical modeling include that models can be over-constrained, as more tracers than unknowns can be integrated together into a single solution. Once the model is built in a computer language/program it can effectively work as a plug and chug method in which only variables need to be altered and the model run, and solutions can be forced to solve for “real world” solutions. For example, an algorithm may be used that prevents negative residence times. The drawbacks of integrating numerical modeling into SGD studies based on using radium as a tracer is the barrier to entry, in that lack of knowledge of programming languages or the initial time cost to building such models in a programming language may be preventing the SGD community from utilizing such models. The goal of this paper is to alleviate that barrier by providing easy to use radium-based numerical models for calculating SGD flux and relative contribution (mixing models), through Open-source Code for Radium-Derived Ocean-Groundwater Modeling (Project Open RaDOM). Multiple variations of each model script are provided as supplemental material or can be run online. A how-to-guide for users not familiar with coding is also attached as supplemental material with a pictorial guide to editing and running a Project Open RaDOM model. Barnes [18] argues that a great way to move science forward is to post code online for other scientists to use, evaluate, and improve upon. The core sentiment that publishing code for SGD models online for other scientists to utilize is at the core of Project Open RaDOM.

2. Model Descriptions

2.1. General Approach

All Project Open RaDOM models are built in R. The reasons for using R as the programming language for this project are that R and RStudio are free to download and use worldwide, it is one of the more common programming languages in the natural sciences, and R scripts can be shared and run online, increasing the usability of these models as scientists need only to access the scripts via a link to use them instead of downloading and installing R and/or RStudio.
In both types of models presented here the equations are rewritten with the knowns on the right-hand side of the equals sign balancing the unknowns on the left-hand side, written in linear notation as:
Ax = b,
where A is a matrix of know activities, x is a vector of unknowns to be calculated (e.g., SGD flux or fraction of groundwater), and b is a vector of known terms (e.g., radioactive decay). All models are solved numerically using the Lawson–Hanson algorithm for non-negative solutions [19,20]. This ensures that only positive solutions or solutions of 0 can be returned by the model, eliminating unreal solutions (like negative residence times). This algorithm also allows for an over-constrained model in which the number of tracers exceeds the number of unknowns. This makes it possible to integrate all tracers into a single model instead of iteratively solving for pairs of tracers yielding multiple different results. The Lawson–Hanson algorithm takes an iterative approach in which the solution is updated for each iteration until the solutions begin to vary little from one iteration to another, converging on the solution and terminating the script, at which point the results are output. In this way the results are approximated.
Residuals are output each time the model is run so that goodness of fit can be evaluated. The closer the residuals are to 0 the better the fit of the model to the data, and the more the results can be trusted. All models are published with an open-source license, which allows other scientists to use, edit, and improve upon the Project Open RaDOM models so long as attribution is given. A request to users who modify these scripts (beyond basic plug and chug functionality) for their studies is that they attach any updated scripts as supplemental material to their publications so that other scientists can use the modified code.
Scripts are provided as supplementary material to this paper, but they can also be accessed and run remotely through this RStudio Online Workspace [21] 10 May 2022. Furthermore, for researchers unfamiliar with R, a how-to-guide on running an example script is provided in the Supplementary Material.

2.2. SGD Flux Models

The SGD flux models of Project Open RaDOM are a simplified version of a system of linear equations from Lecher et al. [22] originally based on the analytical model from Hwang et al. [23]. At their heart they are box models and the simplest version of the model calculates the flux of one source (SGD) and the residence time of the box based off of two radium isotopes.
R a G 1   G = R a B 1 V λ 1 + ( R a B 1 R a O 1 ) V τ ,
R a G 2   G = R a B 2 V λ 2 + ( R a B 2 R a O 2 ) V τ ,
where Ra1 and Ra2 are two radium isotopes observed, V is the volume of the ocean box receiving SGD, RaG, RaB, and RaO are the activities of radium in groundwater, the oceanic box, and offshore water (in dpm/100 L), respectively, τ is the residence time of the ocean box, λ 1 and λ 2 are the decay constants of the radium isotopes, and G is the flux of SGD. The terms of the equation calculate the flux of SGD (term to the left of the equals sign), the radioactive decay of the radium isotopes (first term to the right of the equals sign), and the loss of radium due to advection offshore (rightmost term in the equations). Generally, sources of radium are on the left side of the equals sign while sinks of radium reside to the right of the equals sign.
The model script outputs the flux of SGD in m3/day and the residence time of the oceanic box in days. The remaining terms (radium activities, box volume, and decay constants) are set by the user. The equations can be modified to add additional radium isotopes by adding more rows of equations or to add more sources (such as a river or a second groundwater source) by adding more terms to the left-hand side of the equation. In the Project Open RaDOM script collection, there are scripts for:
  • One source, two radium isotopes;
  • One source, two radium isotopes with example data;
  • One source, three radium isotopes;
  • Two sources, three radium isotopes.
The advantages to this approach for calculating SGD with radium activity data over the box model shown in Equation (1) is that two (or more) radium isotopes can be used in concert to calculate a cohesive SGD flux instead of having two disparate SGD estimates based off of two different isotopes. Additionally, a residence time estimate is not required for the model as it is estimated by the model itself. This model also accounts for radioactive decay as a sink of radium, unlike the model in Equation (1), which may be substantial for the shorter-lived isotopes. Non-radioactive, conservative groundwater tracers (such as barium, or in some cases silicon [24,25]) could also be integrated into the model by setting the radioactive decay constant to 0.

2.3. Mixing Models

The SGD mixing models of Project Open RaDOM are based off of a system of linear equations from Lecher et al. [26] originally based on the analytical model from Moore [17]. The simplest version of the model calculates the fraction of two different water sources (e.g., two different groundwaters or a groundwater and offshore water) to a sampling location, using two radium isotopes:
R a S 1 1 F S 1 + R a S 2 1 F S 2 = R a L 1 ,
R a S 1 2 F S 1 + R a S 2 2 F S 2 = R a L 2 ,
F S 1 + F S 2 = 1   ,
where Ra1 and Ra2 are the radium isotopes utilized, FS1 and FS2 are the fractions of each source to the sampling location, and RaS1, RaS2, and RaL are the radium activities in each source at the sampling location.
The model script outputs the fraction of each water source to the sampling location and a stacked bar graph where each bar represents a sampling location. The radium activities are input by the user. The equations can be modified to add additional radium isotopes (or other conservative tracers) by adding more rows of equations or to add more sources (such as a river or a second groundwater) by adding more terms to the left-hand side of the equation. In the Project Open RaDOM script collection, there are scripts for:
  • Two sources, two radium isotopes;
  • Three sources, two radium isotopes;
  • Two sources, three radium isotopes;
  • Three sources, three radium isotopes;
  • Three sources, four radium isotopes.
The advantages to this approach for calculating fractions of SGD contribution with radium activity data over the standard analytical approach is that the model can be over-constrained, whereby more tracers than unknown fractions can be used to ascertain more accurate source fractions; additionally, short-lived radium isotopes can be more easily utilized than with the analytical solution (more on this in the following section). Alternative conservative tracers to radium can be used in place of radium in this model without alteration as this model does not account for radioactive decay.

3. Examples from the Literature

3.1. SGD Flux Models

The approach of using a system of linear equations and two or more radium isotopes to simultaneously solve for SGD flux and residence time was pioneered by Hwang et al. [23,27]. Both studies used three tracers and a system of three linear equations to solve for SGD flux and residence time. Hwang et al. [23] employed silica, Ra226, and Ra224 whereas Hwang et al. [27] used silica, Ra226, and Ra223. In both cases, the equations were solved analytically and as there were more tracers than unknowns, the equations were iterative solved for pairs of tracers (e.g., silica and Ra226, then silica and Ra224, then Ra226 and Ra224) resulting in three different estimates of SGD flux and residence time for each study. For example, Hwang et al. [27] found SGD flow velocities of 72, 116, and 73 m/y and residence times of 8, 5, and 7 d for the bay they studied in Korea. Multiple flux and residence time estimates leaves a conundrum for the researcher. Which flux estimate should they use to determine the associated flux of SGD solutes like nitrogen or phosphorous (both bioactive nutrients), or to determine biological effects that depend on a residence time estimate?
Lecher et al. [22,26] built upon the equations from these previous studies to develop the method illustrated in Project Open RaDOM’s scripts. These papers quantified the SGD flux and residence time for locations in the Gulf of Alaska, Arctic Ocean, and coastal California using Ra223, Ra224, and Ra228. In the case of coastal California and the Gulf of Alaska, the models also accounted for river flux to the box in addition to SGD flux (the two source, three tracer version) versus the Arctic Ocean, where only SGD flux was calculated (the one source, three tracer version). The main advantage of the numerical approach of Lecher et al. [22,26] versus the analytical approach of Hwang et al. [23,27] is that the latter studies were able to cohesively solve the equations to determine a single SGD flux and residence time for each location.
There are some modifications to the Project Open RaDOM scripts from these versions. The scripts were converted from MATLAB (which requires a licensing fee) to R (which is free). Additionally, the original scripts included an algorithm to determine error based on the distribution of data at each site. Since the distribution of radium in the end-members must be known (e.g., normal, log-normal, uniform, etc.) for that algorithm, and the distributions may vary widely among users, the algorithm to determine error was removed to make the scripts more versatile.

3.2. Mixing Models

The approach of a mixing model based on radium isotopes was pioneered in Moore [17] and Moore [28]. Both of these studies used a mixing model with Ra226 and Ra228 to determine the contribution of two different groundwater sources and offshore ocean water to various locations on the continental shelf of Florida’s Gulf of Mexico and Sicily’s continental shelf, respectively. The use of long-lived radium isotopes was appropriate for these locations as they are more easily detected further from shore than the shorter-lived radium isotopes. Furthermore, due to the low decay rate, these isotopes can be considered conservative over longer distances and lengths of time.
Tamborski et al. [29] later used the same set of equations and tracers as in these two Moore papers to determine the contribution of two different groundwaters and sea water to Long Island Sound. This study went a step further by multiplying the fractions (percentages) output by the model by the volume of water in Long Island Sound and dividing that by a residence time of the sound to calculate an SGD flux, which was comparable to their SGD flux estimates based on mass balances.
At the same time, Moore et al. [30] modified this mixing model to work in the Okatee Estuary of South Carolina. In this case, Ra228 was paired with salinity to determine the contributions of groundwater, bay water (Port Royal Sound), and river water to the estuary. Due to low salinity being an excellent indicator of river water and high Ra228 activities being an excellent indicator of groundwater, this model was robust enough to show how the contributions of these end-members changed throughout the year, as the model was replicated multiple times yearly from 2001 to 2003. Wang et al. [31] applied the same method to calculate the relative contributions of river water, groundwater, and offshore ocean water to a coastal stretch of the Guangdong-Hong Kong-Macau Greater Bay Area.
The original Moore [17] mixing model was modified yet again by Young et al. [32] to work in a lagoon on the Yucatán Peninsula. Young et al. [32] used chlorine concentrations, Ra228, and Ra226 as tracers of two different groundwater sources and ocean water within the lagoon. Applied and solved analytically, as all the papers previously discussed, the mixing model had to be solved iteratively, i.e., chlorine and Ra228, then chlorine and Ra226, and finally Ra228 and Ra226, yielding three different results. For example, at one location in the lagoon, sea water was determined to contribute between ~25% and ~50% of the water to that location, depending on the pair of tracers used. At a second location the contribution of sea water varied between ~30% and ~65% depending on the model. This example highlights one of the biggest issues with the analytical approach: an inability to easily integrate more tracers than unknowns, resulting in wildly different results (up to 50%, as shown in these two examples). However, this also demonstrates the need for integrating as many tracers as possible into one model, so that results are not skewed based on the tracers the researchers just happen to pick. In addition, the researchers of this paper noted a difficulty in using the short-lived radium isotopes in this model, as when used alone the high activities and decay rate of those isotopes cause exceptionally low inputs of sources that were not in concert with estimates from other methods. Furthermore, this paper shows how an additional tracer (chlorine) can also be integrated into the mixing model.
The first successful attempt to over-constrain a radium-based mixing model using numerical methods was presented in Lecher et al. [26]. Lecher used the tracers of Ra224, Ra223, Ra228, nitrate, silicate, and phosphate to calculate the contributions of groundwater, river water, and deep ocean water to a bay in California. In this case, the nutrients and long-lived Ra228 isotope served as excellent tracers of deep water that was depleted in the shorter-lived radium isotopes. By integrating so many different tracers into the mixing model and therein over-constraining the model, the researchers were able to replicate the seasonal and spatial influence of upwelling and riverine flux to the bay in addition to the groundwater flux. These fluxes corresponded well to wind data, as an analogue to the likelihood of upwelling, and river discharge data from the region, despite the non-conservative nature of the shorter-lived isotopes and nutrients employed in the model. This example shows the benefits of a numerical approach that integrates more tracers than sources.
There are two main differences between the Project Open RaDOM scripts and the scripts employed in Lecher et al. [26]. First, the model was transcribed from MATLAB to R. Second, the Project Open RaDOM scripts only allow for up to four tracers instead of six. The scripts could be easily modified to integrate more tracers if needed or desired.

4. Case Studies

4.1. SGD Flux Models

To test the SGD Flux scripts, examples of radium-based box models were collected from the literature and the results of those box models were compared against the results from a Project Open RaDOM script. The results of these case studies are summarized in Table 1. The case studies were limited to articles where enough information was provided so that the end-member activities and volume of the box could be discerned. Even so, for some of the examples, the end-members had to be guessed or back-calculated, as those values were not explicitly stated. So, these comparisons must be viewed with that issue in mind. Hwang et al. [23,27] and Luo et al. [33] provided additional contributions (e.g., diffusive flux from sediments, rivers, etc.) that were added into the scripts since they were provided in the original papers. The script for each model is provided in the supplementary material.
Generally, the SGD values tend to be more conservative when using the Project Open RaDOM scripts than the original box models. The estimates from Project Open RaDOM scripts ranged from being on the same order of magnitude to three orders of magnitude less than the original models. The larger deviations occurred in calculations of an entire bay rather than normalized to area, indicating that smaller deviations may have scaled up to larger differences. The smaller estimates may be an advantage of the Project Open RaDOM models as it is very easy to over-estimate SGD based on radium activities, especially with radium-based residence times [34,35]. Residence times calculated from the Project RaDOM scripts tend to be more variable, tending towards being longer than the other methods employed.
Hwang et al. [23,27] were the most similar case studies, as they also used a combined approach of two radium isotopes, although solved analytically. In both cases the SGD estimate from the Project Open RaDOM models is at least 4× less than the original estimate and the residence time is at least 5 times longer. It is hard to say without independent verification which of the calculations is closer to the true value. At least for Yeoja Bay, an independent residence time from hydro-dynamic modeling exists of 16 days [36]. This is double the Hwang et al. [27] estimate and is about a third of the Project Open RaDOM estimate. As for the seepage rate, Hwang et al. [27] admitted in the original paper that their calculated rates are much higher than typical for coastal areas, which may indicate that the Project Open RaDOM calculated value is more representative. The residence times estimates from the data of Hwang et al. [23] are more similar. However, the SGD estimate of the Project Open RaDOM script is substantially lower. Seepage meter estimates from the shoreline indicate that SGD in those areas can range between 14 and 82 cm/day [37]. Since discharge is typically enhanced at the shoreline and the Project Open RaDOM estimate was obtained by dividing the total flux to the bay by the area of the entire bay (including low or no discharge areas) it may align with those estimates if only the shoreline areas were used for the conversion. For both of these studies the offshore and average bay radium activities was estimated from tables or back-calculated from other values, which if incorrect would obviously skew the model results, causing the different model outputs observed.
For the Shellenbarger et al. [13] comparison, the SGD fluxes were more similar, but again the Project Open RaDOM estimate was the lowest. The Project Open RaDOM script was unable to converge on a solution for the residence time of the box, which in this case was a small section of shoreline in the Red Sea, essentially just the surf zone. The residence time estimates from the original paper came from a heat model, as the radium activities did not decrease enough from the shoreline to calculate a water mass age based on radium isotope ratios. Clearly radium isotopes alone are not the best choice for calculating residence times when the residence times are small, as has been previously discussed, and the volume of the box may have just been too small for the Project Open RaDOM script to converge on a solution.
The case study of Luo et al. [33] shows the only example where the Project Open RaDOM script determined a larger SGD estimate. This is most likely due to a need to balance the much smaller residence time estimate, resulting in faster modeled removal of radium. A residence time of 2.1 h is unrealistic for such a large harbor. The average ocean radium activity was estimated from back-calculations for this system, which may have introduced error, causing the unrealistic residence time. Without knowing more about the details of the study site itself and the correct values for the radium activity end-members, it is hard to discern why the residence time estimate for this case study was so low.
The last case study explored with this model was a section of a Mediterranean karstic coastline, using only the October data from the original study. Garcia-Solonsa et al. [38] used short-lived radium isotope ratios to estimate the residence time of 2.7 days, but as previously discussed this method is inaccurate for calculated residence times of less than 5 days, throwing this residence time into question. It is also of note that the method used to compute residence time actually calculates water mass age, and although this is often used interchangeably in the SGD community with residence time, is not always valid [39]. This validity of this approach was questioned and discussed in-depth in the open access reviews published for that paper [40]. Conversely, the Project Open RaDOM model, which integrated three different isotopes output a much longer residence time of 59 days. The area of coastline this represents is 1.3 km2, and unfortunately there is no independent calculation of residence time to determine which of these two estimates are closer to the true value. As for the SGD calculations, the Project Open RaDOM script again produced a conservative estimate one order of magnitude less than the original, based off of three isotopes used. In the original estimate, the box model was calculated for each isotope individually and then averaged. Again, no independent estimates exist for SGD fluxes in this area, but both the original and Project Open RaDOM calculated fluxes are well within calculated fluxes SGD from other locations.
To summarize the take home points of these case studies: (1) Project Open RaDOM SGD estimates tend to be but are not always conservative when compared to other radium-based box models. This may be an advantage in some situations where SGD may be easily overestimated, but can pose a problem if the SGD flux is unrealistic. At the very least, Project Open RaDOM scripts could be used as a complementary method of calculating SGD to more traditional analytical box models. (2) Knowledge of the study area, appropriateness of this model, and judgement in determining end-members to input to the model is key. The estimated end-members input into the scripts may account for the differences in SGD and residence time estimates between the original studies and the Project Open RaDOM scripts. Orders of magnitude differences in SGD fluxes from the same site can result from different chosen radium activity end-members [39]. In all but two cases the model was successful in calculating realistic-looking values for residence times and SGD. Whether the values output by the models are valid, is up to the discernment of the researchers.

4.2. Mixing Models

To test the Project Open RaDOM mixing model it was compared against the original Moore [17] analytical solution from the Gulf of Mexico and Su et al. [41] from Little Lagoon, Alabama. These papers had the most clearly defined values for the end-members and the output of the mixing models in numerical form, making an accurate reproduction possible. The scripts are available in the Supplementary Materials as Moore2003.R and Su2014.R. The output of the analytical and Project Open RaDOM models for the Moore [17] data is identical (Table 2). The Project Open RaDOM model output was listed to more decimal places, but when rounded it is the same as the Moore [17] results. Most of the results from the analytical and Project Open RaDOM models for the Su et al. [41] results were identical (Table 3). The deviations may come from differences in the end-member for the deep groundwater activity of Ra226. A range was provided for this end-member in the original paper, and the midpoint of that range was used for the deep groundwater end-member in the Project Open RaDOM script. This may not have been the same value used in the original paper, resulting in slightly different results. Generally, the Project Open RaDOM mixing model scripts were able to quickly recreate the results of the analytical models for these two case studies. The utility of the scripts allowing for more end-members and more tracers is yet to be seen as they are used in the SGD community.

5. Limitations and Assumptions

5.1. SGD Flux Models

The most important assumptions for the SGD flux models are that there is little or no diffusion from the sediment. This may not work for all systems, in which sediment diffusion is important, but the model could be adapted to integrate a sediment diffusion term. If a multi-source model is utilized, then a porewater flux term could replace one of the sources to represent a diffusion flux. Since the model utilizes radium isotopes, it captures the brackish component of SGD discharge (fresh groundwater + recirculated seawater) and cannot parse out fresh groundwater flux only. The model will also give inaccurate results if a poorly-scaled matrix is input, i.e., the radium activities are orders of magnitude different. However, this can be assessed using the residuals output from the model.

5.2. Mixing Models

The biggest limitation of the mixing model is that is does not account for radioactive decay or other processes that might serve as sinks for tracers (such as uptake by autotrophs in the case of bioactive elements). However, as shown in Lecher et al. [26], this issue may be overcome by over-constraining the model. The accuracy of this model also suffers if a poorly scaled matrix is employed, but again this may be overcome by changing the units of measurement if different types of tracers are used. Again, the residuals of the fit output by the model can be utilized to assess the accuracy of the model.

6. Conclusions

This paper introduces open-source R scripts for calculating the flux and relative contribution of submarine groundwater discharge to marine systems utilizing the tracers of radium isotopes as Project Open RaDOM. The approach of these scripts builds upon previous methods as they utilize numerical methods that can force solutions based in reality, can over-constrain systems with tracers, and provide a free and user-friendly platform in R. The models can be run using the scripts attached as supplementary material or run in the online workspace. As open-source code, users are free to modify and republish the scripts, encouraging advancement of the field and integration of new techniques.

Supplementary Materials

The following supporting information can be downloaded at: https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/hydrology9060106/s1, Document S1: Project Open RaDOM Visual Instructions; Script S1: BoxModel_2Tracer_1Source_v2022.R; Script S2: BoxModel_2Tracer_1Source_WithExampleData_v2022.R; Script S3: Box_Model_3Tracer_1Source_v2022.R; Script S4: BoxModel_3Tracer_2Source_v2022.R; Script S5: MixModel_2Tracer_2Source_v2022.R; Script S6: MixModel_3Tracer_2Source_v2022.R; Script S7: MixModel_3Tracer_3Source_v2022.R; Script S7: MixModel_4Tracer_3Source_v2022.R, Script S8: MixModel_2Tracer_3Sourcev2022.R; Script S9: Garcia-Solsona2010.R; Script S10: Hwang2005a.R; Script S11: Hwang2005b.R; Script S12: Luo2014.R; Script S13: Shellenbarger2006.R; Script S14: Moore2003.R; Script S15: Su2014.R.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Moore, W.S. The role of submarine groundwater discharge in coastal biogeochemistry. J. Geochemical Explor. 2006, 88, 389–393. [Google Scholar] [CrossRef]
  2. Lecher, A.L.; Mackey, K.R.M. Synthesizing the Effects of Submarine Groundwater Discharge on Marine Biota. Hydrology 2018, 5, 60. [Google Scholar] [CrossRef] [Green Version]
  3. Moore, W.S. The effect of submarine groundwater discharge on the ocean. Ann. Rev. Mar. Sci. 2010, 2, 59–88. [Google Scholar] [CrossRef] [Green Version]
  4. Lecher, A.L. Groundwater Discharge in the Arctic: A Review of Studies and Implications for Biogeochemistry. Hydrology 2017, 4, 41. [Google Scholar] [CrossRef] [Green Version]
  5. Swarzenski, P.W.; Reich, C.; Kroeger, K.D.; Baskaran, M. Ra and Rn isotopes as natural tracers of submarine groundwater discharge in Tampa Bay, Florida. Mar. Chem. 2007, 104, 69–84. [Google Scholar] [CrossRef]
  6. De Sieyes, N.R.; Yamahara, K.M.; Layton, B.A.; Joyce, E.H.; Boehm, A.B. Submarine discharge of nutrient-enriched fresh groundwater at Stinson Beach, California is enhanced during neap tides. Limnol. Ocean. 2008, 53, 1434–1445. [Google Scholar] [CrossRef] [Green Version]
  7. Varma, S.; Turner, J.; Underschultz, J. Estimation of submarine groundwater discharge into Geographe Bay, Bunbury, Western Australia. J. Geochemical Explor. 2010, 106, 197–210. [Google Scholar] [CrossRef]
  8. Mejias, M.; Ballesteros, B.J.; Anton-Pacheco, C.; Dominguez, J.A.; Garcia-Orellana, J.; Garcia-Solsona, E.; Masque, P. Methodological study of submarine groundwater discharge from a karstic aquifer in the Western Mediterranean Sea. J. Hydrol. 2012, 464–465, 27–40. [Google Scholar] [CrossRef]
  9. Taniguchi, M.; Burnett, W.C.; Cable, J.E.; Turner, J.V. Investigation of submarine groundwater discharge. Hydrol. Process. 2002, 16, 2115–2129. [Google Scholar] [CrossRef]
  10. Rama; Moore, W.S. Using the radium quartet for evaluating groundwater input and water exchange in salt marshes. Geochim. Cosmochim. Acta 1996, 60, 4645–4652. [Google Scholar] [CrossRef]
  11. Moore, W.S.; Krest, J. Distribution of 223Ra and 224Ra in the plumes of the Mississippi and Atchafalaya Rivers and the Gulf of Mexico. Mar. Chem. 2004, 86, 105–119. [Google Scholar] [CrossRef]
  12. Crotwell, A.M.; Moore, W.S. Nutrient and Radium Fluxes from Submarine Groundwater Discharge to Port Royal Sound, South Carolina. Aquat. Geochem. 2003, 9, 191–208. [Google Scholar] [CrossRef]
  13. Shellenbarger, G.G.; Monismith, S.G.; Genin, A.; Paytan, A. The importance of submarine groundwater discharge to the nearshore nutrient supply in the Gulf of Aqaba (Israel). Limnol. Oceanogr. 2006, 51, 1876–1886. [Google Scholar] [CrossRef] [Green Version]
  14. Beck, A.J.; Rapaglia, J.P.; Cochran, J.K.; Bokuniewicz, H.J.; Yang, S. Submarine groundwater discharge to Great South Bay, NY, estimated using Ra isotopes. Mar. Chem. 2008, 109, 279–291. [Google Scholar] [CrossRef] [Green Version]
  15. Knee, K.; Street, J.H.; Grossman, E.G.; Paytan, A. Nutrient inputs to the coastal ocean from submarine groundwater discharge in a groundwater-dominated system: Relation to land use (Kona coast, Hawaii, U.S.A.). Limnol. Oceanogr. 2010, 55, 1105–1122. [Google Scholar] [CrossRef] [Green Version]
  16. Knee, K.L.; Garcia-solsona, E.; Garcia-orellana, J.; Boehm, A.B.; Paytan, A. Using radium isotopes to characterize water ages and coastal mixing rates: A sensitivity analysis. Limnol. Oceanogr. Methods 2011, 9, 380–395. [Google Scholar] [CrossRef] [Green Version]
  17. Moore, W.S. Sources and fluxes of submarine groundwater discharge delineated by radium isotopes. Biogeochemistry 2003, 66, 75–93. [Google Scholar] [CrossRef]
  18. Barnes, N. Publish your computer code: It is good enough. Nature 2010, 467, 753. [Google Scholar] [CrossRef]
  19. Lawson, C.L.; Hanson, R.J. Solving Least Squares Problems; SIAM: Philadelphia, PA, USA, 1995; Volume 15, ISBN 0898713560. [Google Scholar]
  20. Mullen, K.M.; van Stokkum, I.H.M. The Lawson-Hanson Algorithm for Non-Negative Least Squares (NNLS). R Package Version 1.4 2015. Available online: https://cran.r-project.org/web/packages/nnls/nnls.pdf (accessed on 10 May 2022).
  21. Lecher, A.L. Project Open RaDOM. Available online: https://rstudio.cloud/content/3560079 (accessed on 10 May 2022).
  22. Lecher, A.L.; Kessler, J.; Sparrow, K.; Garcia-Tigreros Kodovska, F.; Dimova, N.; Murray, J.; Tulaczyk, S.; Paytan, A. Methane transport through submarine groundwater discharge to the North Pacific and Arctic Ocean at two Alaskan sites. Limnol. Oceanogr. 2015, 61, S344–S355. [Google Scholar] [CrossRef]
  23. Hwang, D.W.; Lee, Y.W.; Kim, G. Large submarine groundwater discharge and benthic eutrophication on Bangdu Bay on volcanic Jeju Island, Korea. Limnol. Oceanogr. 2005, 50, 1393–1403. [Google Scholar] [CrossRef]
  24. Oehler, T.; Tamborski, J.; Rahman, S.; Moosdorf, N.; Ahren, J.; Mori, C.; Neuholz, R.; Schnetger, B.; Beck, M. DSi as a Tracer for Submarine Groundwater Discharge. Front. Mar. Sci. 2019, 563. [Google Scholar] [CrossRef]
  25. Lin, I.-T.; Wang, C.-H.; You, C.-F.; Lin, S.; Huang, K.-F.; Chen, Y.-G. Deep submarine groundwater discharge indicated by tracers of oxygen, strontium isotopes and barium content in the Pingtung coastal zone, southern Taiwan. Mar. Chem. 2010, 122, 51–58. [Google Scholar] [CrossRef]
  26. Lecher, A.L.; Fisher, A.T.; Paytan, A. Submarine groundwater discharge in Northern Monterey Bay, California: Evaluation by mixing and mass balance models. Mar. Chem. 2016, 179, 44–55. [Google Scholar] [CrossRef] [Green Version]
  27. Hwang, D.-W.; Kim, G.; Lee, Y.-W.; Yang, H.-S. Estimating submarine inputs of groundwater and nutrients to a coastal bay using radium isotopes. Mar. Chem. 2005, 96, 61–71. [Google Scholar] [CrossRef]
  28. Moore, W.S. Radium isotopes as tracers of submarine groundwater discharge in Sicily. Cont. Shelf Res. 2006, 26, 852–861. [Google Scholar] [CrossRef]
  29. Tamborski, J.; Cochran, J.K.; Bokuniewicz, H.; Heilbrun, C.; Garcia-Orellana, J.; Rodellas, V.; Wilson, R. Radium Mass Balance Sensitify Analysis for Submarine Groudwater Discharge Estimation in Semi-Enclosed Basins: The Case Study of Long Island Sound. Front. Environ. Sci. 2020, 8, 108. [Google Scholar] [CrossRef]
  30. Moore, W.S.; Blanton, J.O.; Joye, S.B. Estimates of flushing times, submarine groundwater discharge, and nutrient fluxes to Okatee Estuary, South Carolina. J. Geophys. Res. Ocean. 2006, 111. [Google Scholar] [CrossRef]
  31. Wang, Q.; Wang, X.; Xiao, K.; Zhang, Y.; Luo, M.; Zheng, C.; Hailong, L. Submarine groundwater discharge and associated nutrient fluxes in teh Greater Bay Area, China revealed by radium and stable isotopes. Geosci. Front. 2021, 12, 101223. [Google Scholar] [CrossRef]
  32. Young, M.B.; Gonneea, M.E.; Fong, D.A.; Moore, W.S.; Herrera-Silveira, J.; Paytan, A. Characterizing sources of groundwater to a tropical coastal lagoon in a karstic area using radium isotopes and water chemistry. Mar. Chem. 2008, 109, 377–394. [Google Scholar] [CrossRef]
  33. Luo, X.; Jiao, J.J.; Moore, W.S.; Ming Lee, C. Submarine groundwater discharge estimation in an urbanized embayment in Hong Kong via short-lived radium isotopes and its implication of nutrient loadings and primary production. Mar. Pollut. Bull. 2014, 82, 144–154. [Google Scholar] [CrossRef]
  34. Kim, G.; Ryu, J.-W.; Hwang, D.-W. Radium tracing of submarine groundwater discharge (SGD) and associated nutrient fluxes in a highly-permeable bed coastal zone, Korea. Mar. Chem. 2008, 109, 307–317. [Google Scholar] [CrossRef]
  35. Wang, X.; Zhang, Y.; Luo, M.; Xiao, K.; Wang, Q.; Tian, Y.; Qiu, W.; Xiong, Y.; Zheng, C.; Li, H. Radium and nitrogen isotopes tracing fluxes and sources of submarine groundwater discharge driven nitrate in an urbanized coastal area. Sci. Total Environ. 2021, 763, 144616. [Google Scholar] [CrossRef] [PubMed]
  36. Lee, D.-I.; Choi, J.-M.; Lee, Y.-G.; Lee, M.-O.; Lee, W.-C.; Kim, J.-K. Coastal environmental assessment and management by ecological simulation in Yeoja Bay, Korea. Estuar. Coast. Shelf Sci. 2008, 80, 495–508. [Google Scholar] [CrossRef]
  37. Kim, K.; Lee, K.; Park, K.S.; Hwang, D.W.; Yang, H.S. Large submarine groundwater discharge (SGD) from a volcanic island. Geophys. Res. Lett. 2003, 30. [Google Scholar] [CrossRef]
  38. Garcia-Solsona, E.; Garcia-Orellana, J.; Masque, P.; Rodellas, V.; Mejias, M.; Ballesteros, B.; Dominguez, J.A. Groundwater and nutrient discharge through karstic coastal springs (Castelló, Spain). Biogeosciences 2010, 7, 2625–2638. [Google Scholar] [CrossRef] [Green Version]
  39. Garcia-Orellana, J.; Rodellas, V.; Tamborski, J.; Diego-Feliu, M.; van Beek, P.; Weinstein, Y.; Charette, M.; Alorda-Kleinglass, A.; Michael, H.A.; Stieglitz, T.; et al. Radium isotopes as submarine groundwater discharge (SGD) tracers: Review and recommendations. Earth Sci. Rev. 2021, 220, 103681. [Google Scholar] [CrossRef]
  40. Garcia-Solsona, E.; Garcia-Orellana, J.; Masqué, P.; Rodellas, V.; Mejías, M.; Ballesteros, B.; Domínguez, J.A. Interactive comment on “Groundwater and nutrient discharge through karstic coastal springs (Castelló, Spain)” by E. Garcia-Solsona et al. Biogeosciences Discuss. 2010, 7, C211–C215. [Google Scholar]
  41. Su, N.; Burnett, W.C.; MacIntyre, H.L.; Liefer, J.D.; Peterson, R.N.; Viso, R. Natural Radon and Radium Isotopes for Assessing Groundwater Discharge into Little Lagoon, AL: Implications for Harmful Algal Blooms. Estuaries Coasts 2014, 37, 893–910. [Google Scholar] [CrossRef]
Table 1. Comparison of results of original box models and Project Open RaDOM models using the same data. Project RaDOM units were converted to the units of the original study.
Table 1. Comparison of results of original box models and Project Open RaDOM models using the same data. Project RaDOM units were converted to the units of the original study.
OriginalProject RaDOM
StudyIsotopesResidence TimeSGDIsotopesResidence TimeSGDSI Script
Hwang et al. [23]Ra226 and Ra223
Combined
2.0 d45–48 cm/dRa226 and Ra2239.8 d1.7 cm/dHwang 2005a.R
Hwang et al. [27]Ra226 and Ra223
Combined
8 d72 m/yRa226 and Ra22355.9 d18.3 m/yHwang 2005b.R
Shellenbarger et al. [13] (winter only)Ra223 and Ra224
Separate
3–6 h19 and 27 L/sRa223 and Ra224Un 115.3 L/sShellenbarger 2006.R
Luo et al. [33]
(upper layer only)
Ra223 and Ra224
Separate
7.2 d4.6 × 105 and 6.4 × 105 m3/dRa223 and Ra2242.1 h1.0 × 107 m3/dLuo 2014.R
Garcia-Solonsa 2010Ra224, Ra223, Ra226, and Ra228 Separate then
Averaged
2.7 d71,500 ± 11,200 m3/dRa223, Ra226, and Ra22864 d2750 m3/dGarcia-
Solosna 2010.R
1 Model was unable to converge on a solution.
Table 2. Comparison of results of original mixing model from Moore [17] and Project RaDOM models using the same data. fns is the fraction from nearshore groundwater, fos is the fraction from offshore groundwater, and fo is the fraction of ocean water. Reprinted/adapted with permission from Ref. [17]. 2003, Springer Nature.
Table 2. Comparison of results of original mixing model from Moore [17] and Project RaDOM models using the same data. fns is the fraction from nearshore groundwater, fos is the fraction from offshore groundwater, and fo is the fraction of ocean water. Reprinted/adapted with permission from Ref. [17]. 2003, Springer Nature.
Moore [17]Project RaDOM
Eastern Transect Stationfnsfosf0fnsfosfo
S40.070.020.910.0720.0170.91
S50.030.000.970.0300.0040.965
S60.020.000.980.02300.965
S70.020.000.980.0160.0040.978
S80.030.000.970.0270.0020.969
S80.010.010.980.0130.0070.978
S80.030.000.970.0300.0040.965
Table 3. Comparison of results of original mixing model from Su et al. [41] and Project RaDOM models using the same data. fSGW is the fraction from shallow groundwater, fDGW is the fraction from deep groundwater, and fSW is the fraction of sea water. Reprinted/adapted with permission from Ref. [41]. 2013, Springer Nature.
Table 3. Comparison of results of original mixing model from Su et al. [41] and Project RaDOM models using the same data. fSGW is the fraction from shallow groundwater, fDGW is the fraction from deep groundwater, and fSW is the fraction of sea water. Reprinted/adapted with permission from Ref. [41]. 2013, Springer Nature.
Su et al. [41]Project RaDOM
Surface Water SamplingfSGWfDGWfSWfSGWfDGWfSW
5 April 20100.070.050.880.0640.0530.882
15 November 20100.070.100.830.0460.0950.858
9 March 20110.110.110.780.1100.1100.778
29 November 20110.020.150.830.0080.1440.847
24 April 20120.220.320.460.2120.3220.465
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lecher, A.L. Open-Source Code for Radium-Derived Ocean-Groundwater Modeling: Project Open RaDOM. Hydrology 2022, 9, 106. https://0-doi-org.brum.beds.ac.uk/10.3390/hydrology9060106

AMA Style

Lecher AL. Open-Source Code for Radium-Derived Ocean-Groundwater Modeling: Project Open RaDOM. Hydrology. 2022; 9(6):106. https://0-doi-org.brum.beds.ac.uk/10.3390/hydrology9060106

Chicago/Turabian Style

Lecher, Alanna L. 2022. "Open-Source Code for Radium-Derived Ocean-Groundwater Modeling: Project Open RaDOM" Hydrology 9, no. 6: 106. https://0-doi-org.brum.beds.ac.uk/10.3390/hydrology9060106

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop