Conceptual uncertainties in groundwater and porewater fluxes estimated by radon and radium mass balances

Radium isotopes and radon are routinely used as tracers to quantify groundwater and porewater fluxes into coastal and freshwater systems. However, uncertainties associated with the determination of the tracer flux are often poorly addressed and often neglect all the potential errors associated with the conceptualization of the system (i.e., conceptual uncertainties). In this study, we assess the magnitude of some of the key uncertainties related to the determination of the radium and radon inputs supplied by groundwater and porewater fluxes into a waterbody (La Palme Lagoon, France). This uncertainty assessment is addressed through a single model ensemble approach, where a tracer mass balance is run multiple times with variable sets of assumptions and approaches for the key parameters determined through a sensitivity test. In particular, conceptual uncertainties linked to tracer concentration, diffusive fluxes, radon evasion to the atmosphere, and change of tracer inventory over time were considered. The magnitude of porewater fluxes is further constrained using a comparison of independent methods: (1) 224Ra and (2) 222Rn mass balances in overlying waters, (3) a model of 222Rn deficit in sediments, and (4) a fluid‐salt numerical transport model. We demonstrate that conceptual uncertainties are commonly a major source of uncertainty on the estimation of groundwater or porewater fluxes and they need to be taken into account when using tracer mass balances. In the absence of a general framework for assessing these uncertainties, this study provides a practical approach to evaluate key uncertainties associated to radon and radium mass balances.

outflows from various pathways and thus smooth out spatial heterogeneities (Charette et al. 2008). The tracer flux supplied by porewater can be estimated by measuring the tracer concentration in the studied system and accounting for all the other potential tracer sources (e.g., diffusion from sediments, sediment resuspension, riverine inputs) and sinks (e.g., radioactive decay, export offshore, evasion to the atmosphere). This tracer flux is then converted to a porewater flux by dividing the tracer flux by the tracer concentration in the discharging fluids (i.e., the porewater endmember) (Charette et al. 2008;Cook et al. 2018b).
The characterization of both the tracer flux and the endmember tracer concentration are thus critical components of tracer-derived estimates of porewater fluxes. While potential errors related to the determination of tracer endmembers have received significant attention in the literature (Michael et al. 2011;Cho and Kim 2016;Cerdà-Domènech et al. 2017;Duque et al. 2019), uncertainties associated with the determination of the tracer flux are still poorly addressed (Burnett et al. 2007). Indeed, most of the uncertainty assessments attribute all sources of uncertainty to analytical errors or spatial and temporal variability of model parameters (e.g., tracer concentration, wind speed, water volume) (e.g., Alorda-Kleinglass et al. 2019;Webb et al. 2019). However, as an abstract representation of natural systems and processes, the conceptualization of a model in itself is a considerable source of uncertainty, the so called conceptual uncertainty (Enemark et al. 2019). In the case of radon and radium mass balances, these conceptual uncertainties are often linked to (1) intrinsic assumptions of the model used (e.g., the system is at steady state, all the outputs occur at an average concentration) or (2) the choice of the equation/method used to quantify the source and sink terms (e.g., multiple possible approaches to determine diffusive fluxes or mixing losses, several winddependent empirical equations to estimate radon evasion to the atmosphere). Whereas these conceptual uncertainties are sometimes acknowledged as the main source of uncertainty on tracer flux estimates, they are rarely considered in the quantitative uncertainty assessment. Neglecting conceptual uncertainties may introduce a significant bias to the mass balance results and lead to overconfidence in mass balancederived flux estimates.
In this study, we quantify some of the key uncertainties related to the determination of the radium and radon inputs supplied by porewater fluxes into a coastal lagoon, including some of the often-overlooked conceptual uncertainties. These uncertainties are addressed here through single ensemble modeling, whereby a single model is run multiple times with variable sets of initial parameters (Uusitalo et al. 2015). We estimate the most relevant parameters of the tracer mass balance by combining multiple approaches and varying assumptions and we evaluate their influence on porewater flux estimates. We further constrain the magnitude of porewater flux estimates and their uncertainties by a multimodel approach that compares four independent approaches: (1) 224 Ra and (2) 222 Rn mass balances in overlying waters, (3) a model of 222 Rn deficit in sediments, and (4) a fluid-salt numerical transport model.

Study site: La Palme Lagoon, France
La Palme Lagoon (NW Mediterranean, France) (Fig. 1) is a small (500 ha surface area), shallow coastal lagoon, with mean and maximum water depths of $ 0.7 m and $ 2 m, respectively and is unaffected by tidal forcing. Most of the lagoon is covered by fine-to-coarse grained sands (100-500 μm) except the northern part of the northern basin that is dominated by fine-grained sediments ($ 50 μm). This study focused on the northern basin of the lagoon, which has restricted flushing from the southern basin due to dikes and represents the main waterbody ($ 85% of the lagoon area and > 95% of the total water volume; Stieglitz et al. 2013). This limited exchange, together with high evaporation rates in summer, often leads to hyper-saline conditions in the lagoon (> 60; Rodellas et al. 2018).
La Palme Lagoon receives continuous brackish groundwater inputs from a regional karst aquifer (3× 10 3 to 25 × 10 3 m 3 d −1 ; Stieglitz et al. 2013;Bejannin et al. 2017;Rodellas et al. 2018), chiefly from a karstic spring connected to the lagoon via a small stream in the northwestern part of the lagoon (Wilke and Boutière 2000;Monnin et al. 2019). For the purpose of this study, the input of low-salinity groundwater through this spring-fed stream is treated as a surficial input and thus it is not included within porewater fluxes. Some small ephemeral streams also supply freshwater to the lagoon, but they cease to flow soon after rainfall events and are typically inactive during summer months (the period of this study).
Previous studies based on radon and radium tracers have been conducted in the lagoon to estimate the magnitude of porewater fluxes (or lagoon water circulation) (Stieglitz et al. 2013;Bejannin et al. 2017;Rodellas et al. 2018;Tamborski et al. 2018). Estimated porewater fluxes vary seasonally and range from 20 × 10 3 to 90 × 10 3 m 3 d −1 , largely depending on their temporally variable driving forces, chiefly lagoon water levels and locally generated wind waves (Cook et al. 2018a;Rodellas et al. 2020).

Radon and radium data
Most of the radon and radium data used in the mass balance is derived from two previous studies that were conducted concurrently in La Palme Lagoon in June 2016 (Rodellas et al. 2018 for 222 Rn and Tamborski et al. 2018 for 224 Ra) (Fig. 1). This data include 222 Rn and 224 Ra concentrations in (1) lagoon surface waters, (2) pore waters collected at four different locations (Pz1, Pz2, Pz3, and Pz4) and at different depths (from 5 to 50 cm below the sediment-water interface), and (3) discrete water samples from the permanent stream that flows into the lagoon and is fed by the main karstic groundwater spring (Fig. 1). Additional data was also collected to establish the short-term temporal variability of the short-lived tracers 222 Rn and 224 Ra in the lagoon and to determine the porewater 222 Rn concentration in equilibrium with production from sediments. A description of the data used in the mass balances and the raw data are available as Supporting Information (S1 and ResearchData, respectively).

Radon and radium mass balances
Identified sources and sinks of 222 Rn and 224 Ra for the northern section of La Palme Lagoon are shown in Fig. 2. Based on these box models, mass balance equations for 222 Rn (Eq. 1.1) and 224 Ra (Eq. 1.2) can be written as: where V and A are the water volume (m 3 ) and surface area (m 2 ) of the northern basin, respectively; t is time in days (d); Q str is the water inflow from streams (including the spring-fed stream in the northwestern part of the lagoon) (m 3 d −1 ) and Q out is the water outflow to the southern lagoon (m 3 d −1 ); C and C str are the mean concentration of 222 Rn and 224 Ra in lagoon waters and in streams (Bq m −3 ), respectively; C Ra is the mean concentration of dissolved 226 Ra in lagoon waters (Bq m −3 ); X sed and C sus are the amount of radium available for desorption from resuspended sediments (Bq kg −1 ) and the concentration of resuspended sediments in lagoon waters per day (kg m −3 d −1 ), respectively; F dif and F PW are the net fluxes of 222 Rn and 224 Ra per unit area from underlying lagoon sediments due to molecular diffusion and porewater fluxes (Bq m −2 d −1 ), respectively; λ is the radioactive decay constant of the specific tracer ( 222 Rn and 224 Ra) (d −1 ); and k is the gas transfer velocity (m d −1 ). The definitions of all these terms are summarized in Table 1.

Initial estimation of mass balance fluxes
An initial estimate of the fluxes in the mass balance is needed to subsequently assess the sensitivity of porewater flux estimates to different mass balance parameters. Mass balance equations (Eqs. 1.1, 1.2) are solved analytically for 222 Rn and 224 Ra, by initially assuming that the system is in steady state ( ∂CV ∂t = 0). The different parameters in Eqs. 1.1, 1.2 are determined using previously estimated values or applying common approaches in groundwater/porewater studies. Values of V, A, Q out , C Ra , X sed , and C susp for 222 Rn and 224 Ra are obtained directly from Rodellas et al. (2018) and Tamborski et al. (2018) ( Table 1). The main surface inputs to the lagoon (Q str ) occurred via the spring-fed stream, which was directly measured by Tamborski et al. (2018) with a flow meter in the 2016 dry season. The tracer concentration in the karstic spring is used as C str , assuming that there were no major tracer losses (radioactive decay and evasion to the atmosphere) or  (Rodellas et al. 2018;Tamborski et al. 2018). Dots represent sampling points used to derive the interpolation (kriging). Note that the different gradients observed for salinity, 222 Rn and 224 Ra distributions are mainly caused by the different sources and sinks that control the distribution of these parameters.
additional inputs between the spring and the lagoon over the $ 250 m length of the stream (Bejannin et al. 2017).
To estimate the weighted average 222 Rn and 224 Ra concentrations in lagoon waters (C), we determine the tracer inventories by multiplying the 222 Rn and 224 Ra concentrations at each station by the respective water depth and subsequently dividing averaged inventories by the average water depth in the lagoon. The diffusive fluxes of 222 Rn and 224 Ra from lagoon sediments are obtained directly from Rodellas et al. (2018) and Tamborski et al. (2018), who used experimental approaches frequently applied in the literature: 222 Rn diffusive fluxes were estimated from a depth-independent approach based on the radon concentration in equilibrium with sediments and 224 Ra diffusive fluxes were estimated from sediment incubations. The radon gas transfer velocity (k) is estimated from the wind-dependent empirical equation of MacIntyre et al. (1995) normalized for radon in seawater and the specific water temperature. To obtain a representative value for k, we use the average wind speeds recorded in the 48 h preceding the final sampling period (Rodellas et al. 2018). Solving the tracer mass balance equations using these parameters yields porewater-driven tracer fluxes (F PW ) of 70 Bq m −2 d −1 for 222 Rn and 1.1 Bq m −2 d −1 for 224 Ra ( Table 1). As a first approximation of the uncertainties of F PW , we propagate parameter uncertainties derived from analytical uncertainties, standard errors of the mean or from the literature. This results in a relative uncertainty associated with F PW of $ 35% for both 222 Rn and 224 Ra (Table 1).

Sensitivity analysis and single model ensemble
We conducted a sensitivity analysis to determine the parameters to which the mass balance is most sensitive to. A relative variability of ± 25% is attributed to the estimated value for each parameter in Table 1. Because lagoon water volume (V) is related to lagoon area (V = AÁh, where h is the water depth), we evaluated the sensitivity of porewater fluxes to h rather than to V. Given that the initial estimation of ∂CV ∂t was 0 (steady state), a daily change in tracer concentrations of ± 5% is attributed to this parameter. The key model parameters determined from the sensitivity analysis are subsequently quantified using different approaches that are commonly applied in the literature to determine their potential variability.
To assess the uncertainties associated to porewater fluxes estimated from 222 Rn and 224 Ra mass balances (F PW ), we applied a single model ensemble: a unique model (the tracer mass balance in Eqs. 1.1, 1.2) is run multiple times with the different quantifications of the key model parameters (previously determined in the sensitivity analysis). We conducted a Monte Carlo analysis using 10,000 stochastic simulations for each tracer ( 222 Rn and 224 Ra) where different parameter combinations are produced. In these simulations, (1) the approach used to quantify each key parameter in the mass balance is randomly selected between the different methods evaluated, and (2) model parameters are generated via normally distributed random values based on the mean and standard deviation (SD) of each parameter. The probability density function that best fits the data distribution obtained from the Monte Carlo simulations was determined using least-squares fitting (SciPy library for the Python programming language). Statistical distributions tested included Beta, Binomial, χ 2 , Noncentral χ 2 , F, Noncentral F, Gamma, Negative Binomial, Normal, Poisson, Student's t, and Noncentral t. Distribution parameters derived from the best fit to the data were used to derive the expected (mean and SD) porewater fluxes for 222 Rn and 224 Ra.

Results and discussion
Sensitivity of porewater fluxes to mass balance parameters The results of the sensitivity analysis are shown in Fig. 3. Porewater fluxes (F PW ) are not sensitive to changes in lagoon area (A) because both relevant input and output terms are scaled to A (see Eqs. 1.1, 1.2), counteracting the effects of potential changes on this parameter. Porewater fluxes are highly sensitive to water depth for 224 Ra because the main sink term (radioactive decay) depends on h, but none of the relevant sources do. For 222 Rn, porewater fluxes are less sensitive to changes in h because its main sink is gas evasion to the atmosphere, which does not depend on water depth. Both A and h (or V) can however represent an important source of uncertainty in other shallow systems with intricate and variable shapes and bathymetries, and thus they need to be accurately determined (Trigg et al. 2014). In the case of this study, A and h were determined from satellite images and a detailed bathymetry map by interpolating echosounder data accurately georeferenced with a differential GPS (Rodellas et al. 2018), and thus we assume that these parameters are well constrained and they are not further discussed here. Q str , C str , Q out , C Ra , X sed , and C sus can be sometimes difficult to accurately constrain. However, errors on their quantification for La Palme Lagoon produce relatively small changes on porewater flux estimates (a change of 25% produces errors < 5% in F PW ).
Sensitivity to average tracer concentrations in lagoon waters (C) and the diffusive flux from sediments (F dif ) is high . A relative variability of ± 25% is attributed to each parameter, with the exception of ∂C ∂t (a daily change in tracer concentrations of ± 5% is attributed to this parameter).
for both radon and radium estimates, and thus their accurate estimation is important to obtain reliable estimates for porewater fluxes. In the case of radon, the results are also highly sensitive to the gas exchange velocity (k), for example, a change of 25% in k will produce a change of $ 30% in F PW . The final estimates are also highly sensitive to the assumption of steady-state conditions. A daily change in tracer concentrations ( ∂C ∂t ) of 5% will result in relative errors in the final F PW estimates of $ 7% for 222 Rn and as much as $ 60% for 224 Ra.
Relative uncertainties associated with mass balance parameters can greatly vary depending on the parameter and the approach used to quantify it (from < 10% to > 100%). For instance, recent studies estimating porewater fluxes using 224 Ra or 222 Rn have attributed relative uncertainties of 10-90% to desorption from resuspended sediments (VÁX sed ÁC sus ), 20-50% to diffusion from sediments(F dif ÁA), 0-30% to surface water inputs (Q str ÁC str ), 20-100% to atmospheric evasion (kÁAÁC), 20-100% to export offshore (Q out ÁC), and 10-70% to radioactive decay (λÁVÁC) (Wong et al. 2013;Garcia-Orellana et al. 2014;Luek and Beck 2014;Rocha et al. 2016;Sadat-Noori et al. 2017;Chen et al. 2020). Figure 4 depicts the uncertainty introduced to final estimates of F PW assuming that total inputs (excluding F PW ) or outputs have an arbitrary relative uncertainty of 25%. A higher relative contribution of porewater fluxes as a tracer source into the studied system (i.e., the ratio F PW ÁA/total inputs) will always result in a lower relative uncertainty to the final fluxes. For example, total outputs with a relative uncertainty of 25% will always produce estimates of F PW with an associated uncertainty higher than > 100% if the relative contribution of porewater fluxes as a tracer source is lower than 25% (Fig. 4). The importance of specific input (e.g., diffusion from sediments, inputs from streams or rivers) and output fluxes (e.g., radioactive decay, export offshore, or evasion to the atmosphere) as a source of uncertainty additionally depends on their importance relative to other tracer sources and sinks in the studied system. For instance, F PW might be highly sensitive to tracer inputs from streams or rivers in areas with significant surface water inputs (e.g., estuaries) (Luek and Beck 2014), whereas it might be highly sensitive to diffusive fluxes in shallow areas covered by fine-grained sediments rich in uranium and thorium series nuclides (Rodellas et al. 2015b). Similarly, the export of tracer offshore in open systems (e.g., coastal ocean) is likely to represent the main 222 Rn and 224 Ra loss term in the mass balance and the porewater flux estimates will be highly sensitive to the estimation of this term (Burnett et al. 2007).

Conceptual uncertainties intrinsic to the mass balance
Aside from the key parameters determined from the sensitivity analysis, there are additional conceptual uncertainties derived from the numerous simplifications and assumptions that are required for the parametrization of the mass balance (e.g., what processes to include and to neglect, how to parametrize them, etc.). These errors in approximating the real system are intrinsic to each model and cannot be addressed through the evaluation of model parameters conducted in this study (Tebaldi and Knutti 2007). Constraining all these structural uncertainties is beyond the scope of this study. However, we would like to highlight here a key general assumption in the mass balance approach: all the model parameters are assumed constant with respect to the timescale of the tracer used. If transient processes occur during the time that the tracer resides in the system (e.g., temporally variable stream inputs, change of wind conditions that influence radon evasion to the atmosphere), this could lead to significant errors in porewater flux estimations and transient models should be implemented instead (Gilfedder et al. 2015). All studies based on tracer mass balances should therefore estimate the timescale of the tracer applied to assess the values that best represent the transient parameters over this timescale. The tracer residence time can be calculated by dividing the tracer inventory in surface waters by the sum of all tracer losses. Based on the 222 Rn and 224 Ra mass balances (Eqs. 1.1, 1.2), the average tracer residence time in the system (t tracer ) can thus be estimated as follows: where Q out V can be approximated by the water residence time (τ) if the role of evaporation is assumed to be negligible (description of terms in Table 1). Notice also that k only applies to 222 Rn. The average tracer residence times for 222 Rn and Ra by total inputs (excluding F PW ) or outputs as a function of the relative contribution of porewater fluxes as a tracer source (F PW ÁA/Total inputs). An arbitrary relative uncertainty of 25% is attributed to either total inputs (F in ) or outputs (F out ), which are considered as the only source of uncertainty. This figure is calculated from Eqs. 1.1 and 1.2 by assuming steadystate conditions ( ∂CV ∂t = 0) and thus recognizing that isotopes for study sites with different characteristics are shown in Fig. 5. Using the data in Table 1 and Eq. 2, the average 222 Rn and 224 Ra residence time in La Palme Lagoon is $ 1 d and $ 5 d, respectively. The different parameters in the mass balances should thus be representative of this indicative timescale. Constraining the tracer timescale is also particularly relevant to validate the assumption of steady state that will be discussed below.

Conceptual uncertainties associated with mass balance parameters
The critical parameters for 222 Rn and 224 Ra mass balances in La Palme Lagoon determined from the sensitivity analysis ( Fig. 3) are: (1) average tracer concentration in surface waters (C); (2) diffusive flux of tracer from sediments (F dif ); (3) gas exchange velocity for radon (k); and (4) change of tracer inventory over time ( ∂CV ∂t ), that is, the assumption of steadystate conditions. Each of these parameters is discussed in the following subsection and is reevaluated using different approaches and assumptions to constrain their uncertainties.

Average tracer concentration in surface waters C
Like in most published studies, the 222 Rn and 224 Ra mass balance used here is a simple zero-dimensional single-box model and thus it assumes that all the output fluxes are derived from C (i.e., radioactive decay, evasion to the atmosphere, mixing with low-concentration waters) (Fig. 2). Therefore, errors on the determination of C will introduce proportional errors in all the tracer output fluxes. This assumption is correct for radioactive decay. However, significant bias can be introduced to tracer estimates when assuming that water outflows at the "boundaries of the box" occur at mean 222 Rn and 224 Ra concentrations. The offshore tracer flux (Q out ÁC) plays a minor role as a sink of 222 Rn and 224 Ra in La Palme Lagoon, and thus this simplification has minor effects on porewater flux estimates, but this can be relevant in systems where water outflow is an important sink term. Assuming that 222 Rn evasion to the atmosphere occurs at C might also introduce some additional uncertainties to the final estimates, because gas exchange will occur at the average concentration immediately below the air-water interface, which might be different from the average concentration within the defined box (C).
Spatial variability of tracer concentrations can compromise the results of the mass balance if the sampling distribution is not adequate and if the spatial interpolation of data does not produce representative results. Rapaglia et al. (2012) applied a location-allocation model to determine the optimal distribution of samples for field collection. Some authors have also used inverse modeling techniques to obtain an average tracer concentration in large regions with low sampling densities (Kwon et al. 2014;Le Gland et al. 2017); however, these approaches are rarely applied in small-scale studies and most of the authors collect a relatively large number of samples to derive average tracer concentrations.
In La Palme Lagoon, a dense sampling grid distributed throughout the lagoon was used to characterize the radon and radium concentrations in lagoon waters (n = 59 for 222 Rn; n = 42 for 224 Ra; Fig. 1). However, we cannot exclude a potential misrepresentation of specific areas: some samples could have been collected in stations not representative of surrounding waters (e.g., station located just above a submerged spring) or the sampling strategy could have missed a relevant tracer source located in areas that were not sampled. Assuming that the sample distribution is adequate to capture local conditions and that each sample is representative of the whole water column could thus introduce additional uncertainties to C.
Here we evaluate only the uncertainties introduced by the spatial-averaging approach by comparing four different commonly applied ways to obtain C: (1) direct averages of measured tracer concentrations (Beck et al. 2007;Tamborski et al. 2017), (2) averages of tracer inventories divided by average water depth Rodellas et al. 2017), (3) averages from a deterministic interpolation method (inverse distance weighting [IDW]; using a power parameter P of 5) (Sadat-Noori et al. 2017;Correa et al. 2020); and (4) averages from a geostatistical interpolation method (kriging) (Gu et al. 2012;Baudron et al. 2015). It should be noted that results from interpolation methods can also vary depending on the parameters and analysis used (e.g., maximum C changes of $ 15% when using different power parameters for IDW). Since sampling stations were homogeneously distributed, results obtained from the three latter approaches are similar, with C ranging from 112 to 136 Bq m −3 for 222 Rn and 19 to 23 Bq m −3 for 224 Ra. However, obtaining C from the first approach (i.e., direct averages of  concentrations; 174 Bq m −3 and 28 Bq m −3 for 222 Rn and 224 Ra, respectively) leads to significantly greater values for C (Kruskal-Wallis test; p < 0.05). This is likely because direct averages of tracer concentrations will tend to over-represent conditions in shallower areas and, given that shallower areas usually concentrate higher tracer inputs per volume unit, they may often produce an overestimation of C. Excluding the results from direct concentration averages, different approaches produce C varying within 10-20% (Table 2). These conceptual uncertainties, which only consider the spatialaveraging approach and exclude other potential sources of uncertainty (e.g., inappropriate sample distribution), are comparable to the analytical uncertainties associated with the measurements of 222 Rn (e.g., 10-50% in La Palme Lagoon), but considerably higher than those associated with 224 Ra analysis (e.g., 5-10%). These uncertainties might be significantly higher in other systems, particularly when volumetrically large areas of the study site are not considered, when there are strong tracer gradients (e.g., point-sourced springs), when the system is deep and/or the water column is stratified. In this case, the interpolation approach needs to account for spatial differences both in the horizontal and vertical dimensions.
Diffusive flux of tracer from sediments F dif Diffusion of radon and radium isotopes from bottom sediments into the overlying water column may represent a major source of tracers in some systems, particularly in areas covered by fine-grained sediments or sediments that contain high concentrations of 226 Ra and 228 Th (for 222 Rn and 224 Ra, respectively) and/or when inputs of tracers from pore water and other sources are comparatively low (Lambert and Burnett 2003;Rodellas et al. 2015b). Diffusive fluxes of radon and radium are routinely assessed via different methods and approaches, including laboratory incubations of sediment cores (Beck et al. 2007;Tamborski et al. 2017), sediment equilibration experiments (Cable et al. 2004;Wang et al. 2017), tracer concentrations measured in overlying waters (Gilfedder et al. 2015), or the 226 Ra concentration in sediments Baudron et al. 2015). Other studies have directly assumed that the tracer flux from diffusion is negligible (Schmidt et al. 2010;Trezzi et al. 2016). These commonly applied approaches/assumptions are used here to constrain the magnitude of diffusive fluxes of 222 Rn and 224 Ra from sediments of La Palme Lagoon: Table 2. Quantification of the critical parameters in the 222 Rn and 224 Ra mass balances using different approaches. Relative uncertainties of 15%, 25%, 25%, and 25% are assumed for the different choices of C, F dif , k, and dC/dt, respectively (comparable to the common uncertainties associated to these approaches). (i) Incubation of sediment cores in the laboratory following the methods of (Chanyotha et al. 2016) for radon and (Beck et al. 2007) for radium isotopes. Both of these methods monitor the increase of the tracer concentration in waters overlying the sediment core over time. Diffusive radium fluxes based on this approach are directly obtained from  and measured for radon using sediment cores collected at Pz1, Pz2, Pz3, and Pz4 (see S1 section in Supporting Information). The results obtained for individual locations are area-weighted according to their respective sediment distributions (IFREMER 2003) to produce average values for the entire lagoon.

Parameters
(ii) A tracer balance in sediments, where the tracer distribution in sediments is controlled by a balance between production from its parent isotope, diffusion and radioactive decay (Martens et al. 1980;Cook et al. 2008). The tracer diffusive flux (F dif ) can be estimated from this balance (sometimes referred to as a depth-independent approach) assuming a constant concentration upper boundary condition, as follows: where θ is sediment porosity (obtained from Tamborski et al. (2018), C eq and C L are the 222 Rn or 224 Ra concentrations in equilibrium with sediments and in lagoon water (assumed negligible in comparison with C eq ), respectively, and D s is the diffusion coefficient of radon or radium in sediments (estimated from Peng et al. 1974 andGregory 1974, respectively, and corrected for tortuosity following Ullman and Aller 1982). C eq was derived from equilibration experiments for radon (obtained from Rodellas et al. 2018) and from porewater concentrations for radium. Results were obtained from individual sites (Pz1, Pz2, Pz3, and Pz4) and area-weighted to produce absolute fluxes.
(iii) A simplified tracer balance in lagoon waters from the station with the lowest measured activity, assuming that in this station the tracer concentration is only supported by diffusion from sediments and that these inputs are balanced by radioactive decay and (only for radon) degassing to the atmosphere (Gilfedder et al. 2015). This implies that tracer inputs from other sources (e.g., porewater exchange) are negligible at the station with the lowest activities, and therefore this approach estimates maximum diffusive fluxes, therein resulting in conservative porewater fluxes (Gilfedder et al. 2015). The diffusive flux can be calculated from the minimum tracer activities in the lagoon (C min ) as follows: where z is the water depth at the station with the minimum tracer activity, λ is the 222 Rn or 224 Ra decay constant, and k is the gas transfer velocity.
(iv) (Only for 222 Rn) An empirical linear relationship between 226 Ra concentration in sediments and radon diffusion from sediments derived from . Parameters other than the radium content in sediments may influence the diffusion of radon (e.g., grain-size, temperature, porosity), and so this approach is likely to be site-specific and subject to significant error. Nevertheless, it provides a useful first approximation. The equation calculates the diffusive radon flux (in Bq m −2 d −1 ) as: where C Ra226-sed is the 226 Ra concentration in bulk sediments (Bq kg −1 ), taken as the average concentrations in sediments from La Palme Lagoon from (Stieglitz et al. 2013) (C Ra226-sed = 32 Bq kg −1 ).
(v) Assuming that the flux of radon and radium from diffusion is negligible in comparison to tracer inputs from other sources.
The tracer diffusive fluxes from sediments to La Palme Lagoon varies by a factor of 2-3 for 222 Rn and 224 Ra depending only on the approach or assumption used (Table 2). This comparison demonstrates that the selection of the approach itself can add uncertainties most likely higher than 50% to the estimates of diffusive fluxes of radon and radium. There are additional uncertainties linked to the representativeness of the sediments used for determining absolute fluxes. However, the differences between diffusive fluxes obtained from different types of sediment in the lagoon (when using the same approach) are significantly smaller than the differences in fluxes estimated by different approaches. Therefore, we deduce that the uncertainty associated with spatial variability is minor in comparison to that introduced by the choice of approach used to estimate the diffusive flux.

Gas transfer velocities for radon k
Radon loss to the atmosphere constitutes one of the major sinks in the radon mass balance (Burnett et al. 2007). The determination of the gas transfer velocity (k) is thus one of the most crucial parameters in the estimation of porewater fluxes. Gas transfer velocities in radon studies are generally estimated using empirical equations that estimate k as a function of wind speed (Dulaiova et al. 2010;Dimova and Burnett 2011). Whereas these equations are relatively well parameterized for open ocean conditions (Wanninkhof 2014), using them for shallow open water surfaces at land-ocean interfaces (e.g., lakes, lagoons, estuaries, or the coastal ocean) is expected to add significant uncertainties (Kremer et al. 2003;Cockenpot et al. 2015;Weber et al. 2019). In shallow and small environments, k might be significantly influenced by parameters others than wind speed, such as fetch area, bottom roughness, turbidity, current velocity, or ecosystem size (Borges et al. 2004;Vachon and Prairie 2013).
Another important issue is how to apply these empirical relationships correctly in systems with varying wind speeds. The gas transfer velocity is frequently calculated based on punctual wind speed observations during the sampling campaign or averages of specific time periods (few hours, 1-2 d, etc.). However, errors can be introduced if prior conditions that might have affected observed radon inventories are not considered (Schubert et al. 2019). In addition, wind measurements are commonly made from nearby meteorological stations located on land, and the different roughness lengths (i.e., friction) of vegetated areas vs. open water complicate extrapolating the land-based wind measurements to open water conditions (Vachon and Prairie 2013).
To evaluate the conceptual uncertainties associated with the determination of gas transfer velocities (k) for La Palme Lagoon, we compare how different parameterizations of gas transfer velocities influences the results. First, different ways of calculating wind speeds are used as input parameters for wind-k relationships, using the empirical relationships of MacIntyre et al. (1995) as a base equation. Different empirical and experimental approaches are then compared.
Different wind speeds to estimate k. Hourly data on wind direction and speed at the nearby (< 5 km) meteorological station "Leucate" was extracted from the database of the French meteorological service (Météo France). This assumes that this meteorological data is representative for La Palme Lagoon, introducing another source of error that is difficult to quantify and not explicitly considered here. Representative average wind speeds (scaled to 10 m above the surface) are estimated in four different ways and then subsequently used to derive gas transfer velocities from the wind-k equation of MacIntyre et al. (1995) (considering the indicative 222 Rn residence time in the system of $ 1 d; Eq. 2): (1) average wind speed at La Palme Lagoon during the 6-h sampling period (u = 2.2 m s −1 ); (2) average wind speed during the 24 h before the end of the sampling (u = 2.0 m s −1 ); (3) average wind speed during the 48 h before the end of the sampling (u = 3.5 m s −1 ); (4) weighting the influence of degassing on radon concentration depending on their proximity to the sampling time (i.e., events closer to the sampling campaign are more influential). A weighting factor for the importance of degassing events depending on its proximity to the sampling campaign (w t ) is adapted from (Schubert et al. 2019) based on the radon decay constant (λ) and the gas transfer velocity (k) as primary controls on radon inventories for La Palme Lagoon: where h is the average water depth (0.57 m), t are different time steps before the sampling, and k is estimated by an iterative process.
Results of k obtained applying the MacIntyre et al. (1995) equation but using wind speeds calculated from these different approaches range from 0.24 to 0.59 m d −1 ( Table 2). Given that wind speeds were the lowest on the day of sampling, using k derived from wind measurements concomitant to sampling (or the 24 h before) neglects the radon losses that occurred the days before, whereas k derived from the 48-h approach underrepresents the wind conditions during the sampling time. Differences between approaches thus depend on the specific wind conditions prior to and during a survey. We thus suggest that introducing the weighting factor w t (Eqs. 6.1, 6.2) is the best option to account for degassing events occurring at different times.
Different approaches to estimate k. To compare the estimates of k derived from different approaches, we use four commonly applied empirical degassing equations that relate k with wind speeds in different settings (i-iv), an empirical relationship that introduces the effects of system size (v) and a site-specific experimental estimation of k for La Palme Lagoon (vi). This latter experimental determination is based on a newly developed method to estimate k from continuous measurements of dissolved noble gases (Ar and Kr) (Weber et al. 2019). Weber et al. (2019) applied this method to La Palme Lagoon in a period with similar meteorological conditions to those occurring during the June 2016 sampling (Weber et al. 2019). The different approaches used are summarized in Table 3. All of the equations are corrected via the Schmidt number ratio estimated for radon at the sampling temperature and seawater conditions (Wanninkhof 2014). Using the factor w t to weight the importance of different wind events occurring at different times (Eqs. 6.1, 6.2), the different empirical equations produce gas transfer velocities for La Palme Lagoon ranging from 0.4 to 1.2 m d −1 . This range of Table 3. Approaches used to estimate the gas transfer velocity (k). k 600 and k 660 refer to normalized gas transfer velocities (cm h −1 ) for CO 2 at 20 C for freshwater and seawater, respectively; u is the wind speed (m s −1 ); and LA is the lake area (km 2 ). k is lower than the k estimated experimentally during comparable conditions (k = 1.3 ± 0.2 m d −1 , when normalized for radon) (Weber et al. 2019). The empirical parameterization of the gas transfer velocity itself can thus be in error by up to a factor of four. While these conceptual uncertainties could be reduced by estimating k from site-specific experiments during sampling time (e.g., injections of artificial tracer gases, monitoring of dissolved gases, or floating chambers), this approach is experimentally more demanding and rarely applied.

Change of tracer inventory over time ∂CV ∂t
Most studies assume steady state conditions in radon and radium mass balance models (Beck et al. 2007;Liu et al. 2017;Sadat-Noori et al. 2017). This implies that tracer inputs equal tracer outputs and that these fluxes are constant for the timescale of the tracer. Achieving steady state conditions for 222 Rn would require few days of constant meteorological conditions, as a consequence of the rapidly changing wind-dependent 222 Rn evasion to the atmosphere (Gilfedder et al. 2015). Several studies have indeed suggested that dynamic modeling is often a more accurate representation than steady state mass balances for 222 Rn -based approaches (Cook et al. 2008;Dimova and Burnett 2011;Gilfedder et al. 2015), but this still requires knowledge of the transience of all model parameters. As radium is not sensitive to the gas exchange rate, which is highly temporally variable, it might be reasonable to assume that 224 Ra concentrations in surface waters are less variable than 222 Rn concentrations. However, 224 Ra is generally more sensitive to the steady state assumption than 222 Rn (i.e., a similar variability in the concentration results in a larger error in the porewater flux: Fig. 3). Thus, the change of 224 Ra inventory over time needs to be accurately determined in dynamic water bodies with transient sources or sinks on the timescale of the tracer (e.g., tidal-influenced systems; Sadat-Noori et al. 2015).
To constrain the potential conceptual uncertainties linked to the determination of ∂CV ∂t , we compare two cases: (1) Assuming steady-state conditions ( ∂CV ∂t = 0); and (2) Using variable tracer concentration and lagoon water volume over time ( ∂CV ∂t = C ∂V ∂t + V ∂C ∂t ). The water volume change over time in La Palme Lagoon was measured in June 2016 ( ∂V ∂t = −2.6 × 10 5 ± 1.3 × 10 5 m 3 d −1 , i.e., $ 1% of the lagoon volume per day; Rodellas et al. 2018). Given that the change of tracer concentrations over time was not measured during sampling time, we are using complementary samplings conducted in periods with similar conditions to constrain ∂C ∂t . These numbers thus need to be treated with some caution. Minimum and maximum daily changes of 222 Rn concentrations ( ∂C ∂t ) are obtained from the continuous 5-d monitoring of 222 Rn concentrations (26 April 2016-01 May 2016, using the 1 st and 3 rd quartiles of all the measurements in wind conditions comparable to those found in June 2016; Fig. 6). Representative 224 Ra changes over time ( ∂C ∂t ) are obtained from the change in the weighted average radium concentrations during the two consecutive samplings conducted 12 d apart in 2017; 14.3 Bq m −3 (22 June 2017) and 16.2 Bq m −3 (04 July 2017). To constrain this rate of change, we assume that this change in concentration is occurring over the average time 224 Ra resides in the system (i.e., $ 5 d; Fig. 5). While maximum changes in 222 Rn concentrations over 1 d are $ 40% (comparable to maximum changes observed in other sites with similar characteristics; Dimova and Burnett 2011;Gilfedder et al. 2015), the rate of daily change for 224 Ra is only $ 3%. It should be however noticed that these numbers cannot be directly compared because they are obtained from different approaches (highresolution time series in one station for 222 Rn vs. two spatially distributed samplings few days apart for 224 Ra). While we expected the rate of daily concentration change to be higher for 222 Rn than for 224 Ra, the approach used to estimate ∂C ∂t for 224 Ra is likely significantly underestimating this term.
Changes in volume over time can be neglected when concentration variability is much greater than volume variability (e.g., as in the case of 222 Rn). Changes in tracer inventories over time obtained using these assumptions are reported in Table 2, which range from −2 × 10 7 to 7 × 10 7 Bq d −1 for 222 Rn and from 0 to 3 × 10 5 Bq d −1 for 224 Ra.

Uncertainties of porewater fluxes derived from single model ensemble
To highlight the importance of different parameters as a source of uncertainty to the final estimates, we first estimate porewater fluxes to La Palme Lagoon by only varying one parameter at a time (Fig. 7). Rn-based porewater flux estimates can vary by a factor of $ 2 when different methods are used to evaluate C, F dif or ∂CV ∂t , or up to $ 4 times for the gas transfer coefficient (k) (Fig. 7a). In the case of 224 Ra (Fig. 7b), changing only the approach used to evaluate C or F dif produces final porewater flux estimates ranging by a factor of $ 2, whereas variations on the methods used to estimate ∂CV ∂t produce comparatively lower changes on the final porewater fluxes.
The conceptual uncertainties associated with the final porewater flux estimates are constrained through a single model ensemble using Monte Carlo simulations (see "Methods" section). The data distributions obtained from the simulations are best described by a F distribution (d1 = 369, d2 = 18.7, loc = −89.4, scale = 161) for 222 Rn and a normal distribution (μ = 2.15, σ = 0.59) for 224 Ra (Fig. 8). Porewater fluxes to La Palme Lagoon estimated from this approach are 91 ± 68 Bq m −2 d −1 for 222 Rn and 2.2 ± 0.6 Bq m −2 d −1 for 224 Ra (mean ± SD of the simulations performed) (Fig. 8). Estimations of porewater fluxes derived from the single model ensemble span a significantly wider range than those derived from the initial direct base estimation, which is the commonly applied approach (70 ± 25 Bq m −2 d −1 for 222 Rn and 1.1 ± 0.4 Bq m −2 d −1 for 224 Ra; Table 1) (Fig. 8). Interestingly, 224 Ra porewater fluxes estimated from the direct quantification are outside and below the range of fluxes obtained with the single model ensemble. This is mainly because the direct estimation was obtained from the lowest average concentration (C obtained from averaging inventories) and the highest diffusive flux (F dif obtained from incubation experiments), which are the approaches that produce the lowest porewater fluxes. This disagreement reveals that the real uncertainties associated with porewater flux estimates are much higher than the propagated uncertainties initially estimated. Relative uncertainties associated with 222 Rn-derived fluxes ($ 70%) are significantly larger than those associated to 224 Ra estimates ($ 30%). These uncertainties should not be directly compared because different approaches have been used to evaluate the different mass balance parameters for 222 Rn and 224 Ra. However, higher uncertainties are expected for 222 Rn in La Palme Lagoon as a consequence of the large uncertainties associated with the determination of the gas transfer velocity (k), which only affects 222 Rn. This might not be the case for other systems where mixing losses due to exchange with lowconcentration offshore waters (Q out ) play a comparatively major role, because radium mass balances (particularly for the long-lived radium isotopes) are more sensitive to these outflows than radon mass balances are (Tamborski et al. 2020).
The significant uncertainties introduced to ensemble estimates of porewater fluxes are a result of the conceptualization of the approaches and of assumptions made in determining the different model parameters, which are neglected in the direct base estimation. It should be noted that we have determined model parameters using approaches commonly applied in the literature, without determining the appropriateness of each method and thus assuming they are all equally valid. Conceptual uncertainties linked to parameters and final estimates can be significantly reduced if ancillary data (observational or theoretical arguments) are explored to validate if a choice of a method is appropriate. For instance, evidence that direct averages of measured tracer concentrations can overestimate C would suggest not including this approach in the ensemble model. The single model ensemble proposed here is also limited in its ability to capture the total uncertainties of tracer-based porewater fluxes, as there are different ways to design the parametrization of the mass balance and there are many assumptions necessary to develop each approach.

Multimethod comparison: An alternative approach to constrain uncertainties
The evaluation of uncertainties associated with final porewater flux estimates can also be significantly improved by combining different independent methods to quantify these fluxes (multimodel ensemble) (Tebaldi and Knutti 2007;Uusitalo et al. 2015). Comparisons of independent methods (e.g., use of various tracers, seepage meters, hydrological flow modeling) have been applied in porewater flux studies to provide additional confidence in the obtained results (Oberdorfer 2003;Mulligan and Charette 2006;Burnett et al. 2008). However, the differences between estimates obtained using different methods are not directly related to the uncertainties of porewater fluxes because different techniques often capture different components of porewater fluxes or processes occurring over different temporal and spatial scales (Burnett et al. 2006;Taniguchi et al. 2019). A comparison of methods is addressed here to better constrain the magnitude of conceptual uncertainties using this multimodel ensemble. Porewater fluxes obtained from four different and independent models are ensembled here: (1) 222 Rn mass balance in overlying waters; (2) 224 Ra mass balance in overlying waters; (3) 222 Rn deficit in sediments; (4) fluid-salt numerical transport model. To normalize results from the different models, we estimate volumetric water fluxes (which is the end result of most studies).

Rn and 224 Ra mass balances in overlying waters
The ranges of porewater fluxes of 222 Rn and 224 Ra obtained from the single model ensemble are converted into water flows by dividing them by the tracer concentration in discharging fluids (the porewater endmember). Porewater endmembers estimated in Rodellas et al. (2018) and Tamborski et al. (2018) are used here as representative tracer concentrations of pore waters inflowing to La Palme Lagoon in June 2016 (5100 Bq m −3 and 177 Bq m −3 for 222 Rn and 224 Ra, respectively). Neglecting the potentially large uncertainties associated to endmember selection (which are not evaluated in this study), estimated volumetric porewater flows are 1.8 ± 1.3 cm d −1 and 1.2 ± 0.3 cm d −1 for 222 Rn-and 224 Ra-based estimates, respectively. These estimates are in relatively good agreement with each other, suggesting that the magnitude of porewater fluxes is relatively well constrained. However, estimated porewater flows can vary by a factor of 6 ( Fig. 9) as a consequence of the large conceptual uncertainties determined here, which would be considerably higher if other key uncertainties were included in the evaluation, including the uncertainty of the porewater endmember.

222
Rn deficit in sediments The exchange of 222 Rn between overlying waters and sediments produces a deficiency of 222 Rn in sediments relative to the equilibrium porewater concentration that can be used to estimate porewater fluxes (Martin et al. 2007;Cable and Martin 2008;Cook et al. 2018a). Using the 222 Rn profiles reported in Rodellas et al. (2018), we estimated porewater fluxes from this deficit model of 15 ± 3 Bq m −2 d −1 of 222 Rn (weighted average of the four piezometers), which can be converted to a water flow of 0.3 ± 0.1 cm d −1 (see Supporting Information).

Fluid and salt transport model
A vertical one-dimensional salt and transport model was developed to investigate temporally variable porewater fluxes in La Palme Lagoon, and solved using a finite element numerical approach (details in Rodellas et al. 2020). Variations of surface and subsurface salinities and lagoon water levels were continuously monitored at station Pz1 over almost 2 yr (2016)(2017) and coupled with the transport model to produce continuous estimates of porewater fluxes to the lagoon (which were found to be mainly driven by temporal variations of lagoon water depths). Estimated porewater fluxes to La Palme Lagoon for the 5 d previous to the sampling (comparable to 224 Ra residence time) range from 1.7 to 2.0 cm d −1 (from 0.8 to 2.7 cm d −1 if the entire month of June 2016 is considered). ity density function that best fits the data (PDF) is shown, together with the ranges of the mean ± SD of the distribution (Mean ± SD PDF) and of the direct initial base estimation of Table 1 (± SD Direct).
Ensembled porewater fluxes considering the four models range from 0.2 to 3.1 cm d −1 (Fig. 9). The relatively good agreement between individual models provides additional confidence in the obtained results. Assuming that these are four plausible models, average porewater flux estimates could be obtained by integrating all individual model predictions (e.g., 1.7 ± 1.5 cm d −1 ). More complex and adequate evaluations can be achieved by attributing weights (likelihoods) to aggregate multiple model outputs (Rojas et al. 2008), but this weighting will require information on site-specific model performance that is rarely available. However, this multimethod ensemble needs to be addressed in a qualitative way. First, there are large structural assumptions on the application of the ensembled approaches that are not considered here (e.g., tracer endmember selection, constant porosities and 222 Rn production rates, representativeness of selected stations, appropriate parametrization of the models), which could significantly increase the overall uncertainty. Second, the accuracies of the methods are not comparable; whereas the estimates derived from the 222 Rn and 224 Ra mass balances and the fluidsalt transport model are derived from comprehensive models (Rodellas et al. 2018Tamborski et al. 2018), the 222 Rn deficit approach is based on a simple and spatially coarse estimation that could contribute explaining the comparatively lower fluxes obtained from this approach. In addition, the fluid-salt transport model and the 222 Rn porewater profiles are based on upscaling point-measurements, whereas 222 Rn and 224 Ra mass balances constitute integrated whole-of-lagoon approaches. Upscaling point-measurements from unique sediment environments to the whole lagoon basin may constitute an additional source of uncertainty not presently accounted for. Finally, the different models are most likely capturing different components of porewater fluxes. For example, the 222 Rn deficit approach focuses on shallow fluxes, whereas the fluid-salt transport model targets larger scale fluxes.

Conclusions and recommendations
Some of the key uncertainties associated with the estimation of porewater fluxes to La Palme Lagoon by using radon ( 222 Rn) and radium ( 224 Ra) mass balances were evaluated in this study. While some of the results are site-dependent (it is a shallow, semi-enclosed lagoon), the general understanding and conclusions derived from this work provide guidelines to better constrain the magnitude of uncertainties of porewater fluxes estimated from radioactive tracers in coastal and freshwater ecosystems (e.g., lakes, lagoons, marshes, bays, coastal sea). This understanding can therefore be extended to other systems worldwide, as well as to other radium isotopes ( 223 Ra, 228 Ra, and 226 Ra) and tracers (e.g., dissolved silica, CH 4 , salinity).
Uncertainty assessments in radon-and radium-based porewater flux estimates commonly attribute all sources of uncertainty to analytical errors in mass balance parameters or to their natural variability, neglecting conceptual uncertainties implicit in their assumptions. Uncertainties linked to the conceptualization of the mass balance are a primary source of uncertainty of porewater flux estimates. Neglecting conceptual uncertainties may thus lead to unrepresentative estimations of tracer-derived porewater fluxes.
There is not a general framework for assessing these conceptual uncertainties, but in this study, we provide an example of an accessible approach to address uncertainties associated to the estimation of porewater fluxes. We used a single model ensemble where the most sensitive mass balance parameters were identified and subsequently evaluated using different commonly applied approaches. The uncertainties associated with porewater flux estimates can then be quantified by using a set of simulations with a single tracer mass balance but a set of different options for the relevant model parameters. The use of multiple independent methods (e.g., hydrogeological approaches, seepage meters, different tracers) can also contribute to constraining the magnitude of porewater fluxes, although this interpretation needs to be addressed in a qualitative way because the different models may capture different components of porewater fluxes (e.g., shallow vs. deep fluxes, different spatiotemporal integrations, etc.).
We also reviewed different approaches used to evaluate the most sensitive terms in the 222 Rn and 224 Ra mass balances for La Palme Lagoon (C, F dif , k, and ∂CV ∂t ), which are essential components of tracer mass balance models for most coastal systems. Obtaining high quality data for these parameters and appropriately constraining their associated uncertainties is required to produce reliable porewater flux estimates using radon and radium isotopes. Recommendations derived from this study are: • The average radon or radium concentration in surface waters (C) is always a critical component of the mass balance (C directly influences all the tracer output fluxes). An appropriate characterization of C relies on an optimal distribution of sampling stations and on their appropriate spatial averaging. Using direct averages of tracer concentrations to obtain C (a commonly applied approach) will always overrepresent tracer concentrations in shallow areas, and should be avoided. Other averaging methods (e.g., averaging inventories, interpolation methods) are likely to produce more accurate estimates of C, but there can be issues related to the representativeness of sampling stations. Uncertainties associated with interpolation methods should not be neglected. • Diffusive fluxes of radon and radium estimated using different approaches may easily vary by a factor of 2-3. The selection of the approach itself can thus produce tracer diffusive fluxes with large relative uncertainties (most likely higher than 50%). • The determination of the gas exchange rate (k) is a critical source of uncertainty to final 222 Rn-derived porewater flux estimates in shallow systems with relatively long residence times (e.g., coastal lagoons and semi-enclosed basins). The empirical parameterization of k itself can produce k estimates varying by up to an order of magnitude. Unless k can be deduced from site-specific experiments reflecting marine and atmospheric conditions during the time of sampling, the relative uncertainties associated with the determination of this parameter are likely to exceed 50%. Importantly, the determination of representative wind speeds for these parametrizations (u) can also introduce large uncertainties to the determination of k. Introducing a weighting factor to wind speeds occurring at different times is likely the best option for constraining u. • A major assumption in the application of tracer mass balances is that all the model parameters are constant during the time the tracer resides in the system (if a transient model is not used). The potential temporal variability of all the parameters should be considered relative to the residence time of the tracer in the system. This will be tracer and site dependent. • Steady-state conditions for radon are unlikely in systems where gas exchange represents a primary loss term (i.e., shallow lakes, lagoons, and wetlands). Dynamic modeling should therefore be considered for these systems. Although radium is not subject to gas exchange, estimates of porewater fluxes can also be affected by transient sources and sinks (e.g., wind-driven changes on export offshore, temporally variable stream discharges, or porewater fluxes) and this variability should also be considered.
• Radon and radium uncertainties estimated via the approach followed in this study are only capturing some of the key conceptual uncertainties of mass balance estimates. They should be neither considered as absolute uncertainties nor directly compared because they are largely dependent on the approaches considered.