Assessment of turbulence effects on effective solute diffusivity close to a sediment-free fluid interface

Our work is focused on the analysis of solute mixing under the influence of turbulent flow propagating in a porous system across the interface with a free fluid. Such a scenario is representative of solute transport and chemical mixing in the hyporheic zone. The study is motivated by recent experimental results (Chandler et al. Water Res Res 52(5):3493–3509, 2016) which suggested that the effective diffusion parameter is characterized by an exponentially decreasing trend with depth below the sediment-water interface. This result has been recently employed to model numerically downstream solute transport and mixing in streams. Our study provides a quantification of the uncertainty associated with the interpretation of the available experimental data. Our probabilistic analysis relies on a Bayesian inverse modeling approach implemented through an acceptance/rejection algorithm. The stochastic inversion workflow yields depth-resolved posterior (i.e., conditional on solute breakthrough data) probability distributions of the effective diffusion coefficient and enables one to assess the impact on these of (a) the characteristic grain size of the solid matrix associated with the porous medium and (b) the turbulence level at the water-sediment interface. Our results provide quantitative estimates of the uncertainty associated with spatially variable diffusion coefficients. Finally, we discuss possible limitations about the generality of the conclusions one can draw from the considered dataset.


Introduction
The hyporheic region is a major component driving functioning of river ecosystems through regulation of processes associated with, e.g., attenuation of pollutants by biodegradation or adsorption and mixing. All of these processes are critically affected by chemical residence time and mixing rate within the hyporheic zone, these being in turn controlled by turbulent flow patterns which are documented to propagate across the interface between river flow and the hyporheic region (Buss et al. 2009).
Major research findings on stream-subsurface solute exchanges and feedbacks are summarized in Boano et al. (2014) and Rode et al. (2015). Dissolved chemicals transported through stream flow are typically subject to a temporary storage within the hyporheic zone beneath and/or alongside the stream bed followed by subsequent release into the water column. Our current framework of conceptual understanding and depiction of hyporheic exchange processes is grounded on a set of empirical observations and suggests that the characteristics of the stream flow and sediments, along with the stream bottom topography (Elliott and Brooks 1997), are the main drivers to flow and mass exchange across the hyporheic zone (Lautz and Siegel 2006). The hyporheic residence time, i.e., the time of retention of the solutes in the hyporheic zone, is then chiefly controlled by transport characteristics such as pore water velocity and dispersion as well as by the chemical reactivity of the investigated compounds (Schaper et al. 2019). In this context, our study targets the characterization of the uncertainty affecting the parameters describing solute transport, in the hyporheic zone, in the absence of chemical and biological reactive processes. Under these conditions, solute spreading and mixing is driven by molecular diffusion and by the local spatial and temporal fluctuations of the fluid velocity. A description of groundwater flow based on a classical continuum-scale approach (often relying, e.g., on the Darcy equation) is typically employed to study field-scale surface-groundwater interactions. While such an approach allows quantifying water fluxes between subsurface and surface water bodies, it neglects the effects of several physical processes that are known to influence solute residence time and mixing in the hyporheic zone (Woessner 2000). In particular, complexities of hyporheic exchanges are mainly tied to molecular diffusion, shear-driven flow, advective phenomena arising in the presence of river bed forms (also termed as advective pumping) and propagation of turbulence from the river stream into the hyporheic region (O'Connor and Harvey 2008;Chandesris et al. 2013). In practice, approaches employed to estimate hyporheic solute exchanges rely on effective models, somehow embedding the effects of all these processes (Lautz and Siegel 2006). Widely known examples of these effective approaches are the transient storage model (TSM) (Bencala 1983;Hart 1995;Triska et al. 1989;Wörman 2000) the pumping model (PM) (Elliott and Brooks 1997;Packman and Brooks 2001;Packman et al. 2000), and the slip flow model (SFM) (Fries 2007). Ultimately, the aim of all these models is to provide predictions of downstream solute transport, this objective being achieved by embedding specific parameters representing hyporheic solute exchange. In TSM the stream-bed exchange is evaluated considering (a) a storage zone, i.e., the upper region of the bed sediments, which is typically assumed to be characterized by a constant depth, and (b) a mass transfer process, which is formulated through an effective parameter. The exchange is assumed to be proportional to the difference of solute concentrations in the main channel and the storage zone (Bencala 1983). Experimental data such as those reported in Chandler et al. (2016) can be employed to constrain the parameters of TSM, as the latter approach relies on a conceptualization of hyporheic exchange in terms of solute transfer between two physical domains, i.e. main stream and bed sediments. Several mathematical formulations of such a conceptual model have been proposed in Rutherford et al. (1995), Tonina and Buffington (2007), Wörman (1998). In PM and SFM solute transport is assumed to be advection-dominated, flow across the sediment bed being governed by Darcy's Law. The exchange at the sediment-water interface is driven by pressure (i.e., head distribution) or shear velocity (velocity distribution) in PM or SFM, respectively. Molecular diffusion and turbulence effects are neglected in both models. The investigation of the influence of physical quantities such as grain size and bed shear velocity on the diffusive transport across the sedimentwater interface is a significant part of our study. PM describes the solute exchange process as a function of quantities characterizing the stream and sediment attributes, such as sediment permeability, channel geometry and stream velocity (Elliott and Brooks 1997;Marion et al. 2003). Particle-based continuous time random walk (CTRW) approaches have also emerged in the past decade as alternatives to the above-mentioned continuum-based effective approaches (Boano et al. 2007;Sherman et al. 2019).
Our work tackles the characterization of the spatial distribution of mixing and diffusion parameters in the hyporheic zone which is a key element for all of the above mentioned methodologies. To improve the quality of this characterization, propagation of turbulence from the river flow into the hyporheic region has been quantified through experiments (Chandler et al. 2016;Higashino et al. 2009;de Lemos 2005;Roche et al. 2019) and numerical simulations (Bottacin-Busolin 2019; Breugem et al. 2006;Sherman et al. 2019;Chandesris et al. 2013). Notably, in this context we refer to the high resolution experimental data collected by Chandler et al. (2016). Previous studies show the possibility of describing vertical mixing through an effective diffusion model. In these works, the authors quantify the variation of the effective diffusion coefficient with depth below the interface for various bed shear velocity and grain size combinations in the presence of turbulent flow conditions taking place above a flat sediment-water interface. The experimental evidences are interpreted upon relying on the analytical solution proposed by Nagaoka and Ohgaki (1990) and the authors conclude that the diffusion coefficient can be described by an exponential reduction with depth. In this framework, it is noteworthy to observe that also the periodic advection of turbulent eddies across the sediment-water interface exhibits an exponential decrease with depth underneath the stream bed (Bottacin-Busolin and Marion 2010). These results are at the basis of current practices that rely on embedding an exponential reduction (along the vertical) of effective diffusion effects in modeling frameworks (Bottacin-Busolin 2019; Roche et al. 2019). These studies have shown that assuming an exponential reduction of the diffusion coefficient (a) has relevant implications on the modeling of solute mixing in streams, (b) is consistent with model-based interpretations of observed reach-scale transport experiments. Assumptions on the spatial distribution of the diffusion coefficient (e.g., either constant or piece-wise uniform along the sediment bed) can markedly impact the strength of solute exchange and mixing across the surface water interface and thereby influence solute breakthrough times at downstream locations. While the experimental findings of Chandler et al. (2016) have been extensively used within transport model formulations, the uncertainty associated with the interpretation of the dataset has not yet been addressed and quantified.
Our main objective is to analyze the variability of the effective diffusion coefficient within a probabilistic perspective upon relying on a stochastic inverse modeling approach. We recall that stochastic calibration enables one to evaluate the probability density function (pdf) of each (unknown) model parameter conditional to the available dataset, i.e., the posterior parameter density. Otherwise, a deterministic approach focuses on identifying a unique set of model parameters that minimizes a target objective function. Chandler et al. (2016) obtain estimates of the diffusion coefficient below the sediment-water interface by following the procedure outlined in Nagaoka and Ohgaki (1990) and considering a simple mean square error as a metric of model performance, following a typical deterministic optimization procedure.
In our study we consider solute concentration data collected during the experiments of Chandler (2012) (see Sect. 2.1) and ground our stochastic analysis on the acceptance/ rejection methodology (Gelman et al. 2013;Tarantola 2005) to estimate the posterior (i.e., conditional on available observations) probability distribution of the diffusion coefficient which is then propagated to provide a probabilistic depiction of the ensuing chemical concentrations at diverse depths in the system. Our study considers various combinations of bed shear velocity and grain size, as embedded in the dataset analyzed. As such, the key objective of our work is the assessment of the robustness of the typically employed exponential decay relationship between the diffusion coefficient and depth (Chandler 2012;Chandler et al. 2016) and the quantification of the related uncertainty. While, as mentioned above, the assumption of the exponential decrease of effective diffusion (and ensuing mixing) below the stream bed is increasingly adopted in modeling studies, to the best of our knowledge the only experimental evidence directly supporting such a behavior is the one reported in Chandler et al. (2016). The results stemming from our contribution may then be useful to propagate parametric uncertainty to reach-scale models, as recently advocated in Tu et al. (2019).
The structure of the work is described in the following. Section 2 is devoted to the illustration of the methodology employed for stochastic model calibration through the acceptance/rejection method. The reference experimental setup of Chandler et al. (2016) and the analytical solution employed in the calibration workflow are presented in Sect. 2.1. Sect. 3 exposes the key results of the study. Finally, concluding remarks are provided in Sect. 4.

Methodology and problem setup
In the following we describe the setup and experimental data we consider in this study as well as the stochastic model calibration procedure we employ.

Reference experimental setup
We provide here a brief description of the reference experimental tests performed by Chandler (2012) which form the basis of the calibration dataset we consider. The experimental setup relies on an erosimenter of the kind depicted in Fig. 1, which is employed to characterize transport of dissolved chemicals close to the interface between a free fluid system and a porous medium and is typically considered as a proxy for the assessment of mixing processes taking place within the hyporheic zone.
The water column encompassing the main segment of the system is characterized by a height of 300 mm with an internal diameter of 96.2 mm, a 200 mm high column of porous medium being placed at the bottom of the system. Temperature is monitored during the test through a sensor placed at the top of the water column. Fibre-optic fluorometers are used to measure the temporal evolution of solute concentration in the main segment of the device as well as in the porous medium. The latter is fully saturated, fluorometers being positioned at 0.015 m, 0.049 m, 0.083 Turbulence is generated at the top of the base segment through a tri-bladed propeller positioned at a distance of 40 mm from the water-soil interface. The propeller speed is tuned to produce a bed shear velocity favoring a setting corresponding to the onset of sediment motion Chandler et al. (2016). The porous bed is initially saturated with a tracer (fluorescein disodium salt), whose concentration in the water column is null. Various combinations of bed shear velocity (u) and characteristic grain size (d g ) of the porous medium are analyzed in the experiments. Values of the experimental parameters and corresponding estimates of diffusion coefficients are listed in Table 1.

Analytical solution
Assuming that mass transfer can be described as a onedimensional (along the vertical direction) process, solute transport within the porous system is modeled through: Here, C(t, y) [ML À3 ] is solute mass concentration, y [L] is the distance from the exchange interface, D [L 2 T À1 ] is an effective diffusion coefficient, and t is time. We recall that D is considered as an effective model parameter which embeds a continuum-scale description of pore-scale processes. We follow the methodology proposed in Chandler et al. (2016) and employ the analytical solution for Eq. (1) proposed by Nagaoka and Ohgaki (1990) to analyze the depth-dependent variation of the diffusion coefficient. These authors consider the domain to be composed by N L layers. An analytical formulation is then derived for a basic system comprising two layers separated by an interface at y ¼ ÀL and characterized by effective diffusion coefficients D ¼ D i for 0\y\ À L and D ¼ D iþ1 for ðÀL\y\ À 1Þ, respectively (see the sketch in Fig. 2). Boundary and initial conditions for Eq. (1) are set as: Cðt; 0Þ ¼ f ðtÞ ð 3Þ Considering a system of N L layers, solute concentration at the interface between two layers (corresponding, i.e., at y ¼ ÀL in the setting of Fig. 2) can be evaluated analytically as Nagaoka and Ohgaki (1990): Here is the (vertical) distance between two layers, f(t) corresponds to concentration observed at y iÀ1 and t is total duration. Relying on Eq. (7) enables one to estimate the diffusion coefficient of the upper layer (D i ) once the diffusion coefficient at the lower layer (D iþ1 ) and the temporal dynamics of concentration at the top of the upper layer (f(t)) are known. For the lowest layer (i.e., i ¼ N L ), we follow Nagaoka and Ohgaki (1990) and assume that the corresponding diffusion coefficient coincides within the one associated with the layer above it, thus obtaining Eqs. (7), (8) have been employed by Chandler et al. (2016) to estimate values of effective diffusion coefficients D i in the experimental column by replacing f(t) with observations. In Nagaoka and Ohgaki (1990) model calibration is structured according to the following steps: (a) estimation of the value of the diffusion coefficient in the lowest layer upon constraining Eq. (8) through concentration values sampled at the lowest available measurement location and (b) estimation of coefficients D i (once D iþ1 is estimated) through Eq. (7). A schematic depiction of the procedure is offered in Fig. 3. Hereafter we denote the coefficients estimated according to this procedure asD i , to distinguish them from the result of our stochastic model calibration approach.

Acceptance/rejection method
We ground our stochastic model calibration analyses on the acceptance/rejection sampling (ARS) approach (see, e.g., Gelman et al. (2013), Tarantola (2005)). In this framework, our objective is to assess the posterior pdf of the effective diffusion parameter at a given depth y i (See Sect. 2.1), i.e the location of each interface, starting from an assumed prior distribution. In the context of ARS, one aims at obtaining multiple independent Monte Carlo realizations of the model output (in our case, solute concentrations at locations corresponding to monitoring ports) by sampling from the parameter distribution conditional to observations. At each iteration j model parameter values are independently sampled across the support within which the corresponding prior distribution is defined, the analytical solution (Eqs. (7), (8)) is evaluated, and the candidate parameter set is accepted or rejected upon relying on threshold criteria based on the likelihood function. Key details of the ARS are provided in the following. Let C Ã i be a vector whose entries correspond to observed values of concentrations at depth y i for a collection of N Ã discrete time levels and C i;j the corresponding values of concentration obtained by applying Eq. (7) or (8) at depth y ¼ y i and iteration j for the same time levels at which data are available. We randomly sample the assumed prior distribution of diffusion coefficients. The likelihood a j is defined as: where r 2 C y is the variance associated with observation errors (which are assumed to be zero-mean Gaussian).
The workflow depicted in Fig. 4 is employed for the implementation of ARS to all layers except for Layer 1, because no concentration data are available at the interface between the water column and the porous medium.
The procedure is repeated until a collection of N R accepted realizations is obtained or a maximum number of iterations is reached. Accepted values are used to assess the posterior probability distributions of the diffusion coefficients. Hence, as a result of the workflow we obtain a bivariate sample of accepted parameters values, i.e., S i ¼ ½D i ; D iþ1 for each pair of layers (with the exception of the univariate sample D i for i ¼ N L ), from which a sample/ empirical posterior probability distribution can be evaluated. The value of the maximum a posteriori (MAP) can then be approximated on the basis of the mode (i.e., the maximum value) of the posterior. Relying on the MAP is tantamount to considering the most likely value of the investigated parameter within each layer (Murphy 2012).
We point out that reliance on the analytical Eqs. (7), (8) implies that ARS is applied to two adjacent layers, the (1) Fig. 3 Schematic representation of the procedure employed to estimate the effective diffusion coefficient of each layer in the considered experimental setup (see also Nagaoka and Ohgaki (1990)) Stochastic Environmental Research and Risk Assessment (2020) 34:2211-2228 2215 Fig. 4 Workflow of the acceptance/rejection method Tarantola (2005) complete set of results being obtained through application of ARS to all pairs of adjacent layers in the system. It should be noted that we obtain two diverse posterior distributions for parameter D i associated with locations y i (with i [ 2 and i\N L ) because the effective diffusion coefficients are assessed in pairs in the workflow described above. This is consistent with analytical Eqs. (7), (8) so that each iteration of the ARS associates two values for the diffusion coefficients D i with a given layer. These correspond to (a) the one obtained from considering Layers i and i þ 1 and (b) the one stemming from the analysis of Layer i À 1 and Layer i, respectively. To exemplify and clarify this point, we focus here on three consecutive layers (Layer 3, 4, and 5). Once the model has been calibrated for each pair of layers separately, we obtain D 4 1 and D 5 from model calibration on the solute breakthrough curve at the interface between Layer 4 and 5 (i.e., considering C 4 Ã ðtÞ), D 3 and D 4 2 being assessed in a corresponding way upon considering C 3 Ã ðtÞ. We present our results in terms of both MAP values of the diffusion coefficient within each layer. This information can be useful to appraise conditions under which the thickness of a given layer can be considered as large enough to (a) allow process interpretations relying on solutions for unbounded domains or (b) assess the influence of boundary conditions on the system behavior. Criteria for the selection of the representative diffusion coefficient at each layer are clarified in Sect. 3.
The support of the prior distribution of the diffusion coefficients is here centered around the estimated values provided in Chandler et al. (2016) and is taken to span two orders of magnitude. For example, ifD i is of the order 10 À9 m 2 =s, the maximum and minimum value of the support of the prior pdf are taken as 10 À8 m 2 =s and 10 À10 m 2 =s, respectively. Whenever the resulting value D i ðMAPÞ is too close to (or coincides with) either of its limits the support of the prior distribution is widened and the ARS algorithm is restarted. We note that some values of the effective diffusion coefficient are not listed in Chandler et al. (2016), possibly due to some difficulties encountered during model calibration. For those cases we set the width of the support of the prior around the value corresponding to a best fit (as evaluated according to a standard least-square regression) between the analytical solution and the experimental data. Here, we rely on stochastic procedure to obtain posterior distributions of diffusion coefficients also for such layers.
The available experimental data (Chandler et al. 2016) correspond to measurements of the solute concentration at various depths under the sediment-water interface for various combination of grain size and bed shear velocity. For each scenario (sediment diameter/bed shear velocity), the test is repeated twice and the average of the two series of laboratory data is used as input in the acceptance/rejection procedure. Our analysis encompasses the conditions corresponding to the various experimental settings listed in Table 1, as expressed in terms of combinations of grain size (d g ) and bed shear velocity (u).
The experimental data are selected to ensure an accurate representation of the experimental trend. Experimental data show that equilibrium of solute concentrations at locations close to the sediment-water interface is attained at early stage (i.e., t $ 5000s for Test 1 and 2 and t $ 200000s for Test 3 and 4). Hence, model calibration at the interface Layer 2-Layer 3 is performed considering the experimental solute breakthrough curves prior to steady-state to allow for a reduced computational time. Otherwise, solute concentrations at deeper layers exhibit a slow reduction across time, a feature which is mainly seen in the scenario characterized by small sediment size. In order to reduce the computational load without losing accuracy, the experimental data at these locations are selected using a time resolution coarser than the one employed for the collection of the laboratory measurements. No smoothing filters are applied on the experimental data.

Results
We start the illustration of our results by noting that all concentration data are reported in [%], a value of 100 corresponding to the initial concentration in the porous medium. Since no details on measurement uncertainty is reported by Chandler et al. (2016), the measurement error, r C y , is here fixed to 25% with the exception of the calibration scenario associated with the pair Layer 2-Layer 3 for which we consider r C y = 45%. These values have been selected after a series of preliminary analyses to ensure an acceptance rate at least equal to 0.1%. Figures 5 and 6 document the results obtained in terms of the probability distributions for the effective diffusion coefficients. Figure 5 depicts empirical bivariate distributions for the diffusion coefficients associated with the observations pertaining to Test 1 (including the sample posterior distribution for D 5 obtained when considering Eq. (8) to interpret only data from Layer 5), Fig. 6 depicting the resulting concentration histories at various depths. Corresponding results obtained for Tests 2-4 are included in Appendix A. The solute breakthrough curves in Fig. 6 are evaluated from Eqs. (7), (8) where values of diffusion coefficients estimated for each pair of layers are employed. For example, the concentration C 4 at the interface between Layer 4 and Layer 5 is evaluated using the MAP value of the diffusion coefficients D 4 and D 5 resulting from model calibration for the pair Layer 4-Layer 5. ARS. An analogous behavior is documented for the remaining experimental tests for which estimatesD i are available (Fig. 10 in Appendix A). We note that the results depicted in Fig. 5 are obtained upon considering r C y ¼ 45% for the pair Layer 2-Layer 3, as opposed to the value of 25% employed for all of the remaining layers. We note that as a first attempt we set r C y ¼ 25% for all layers. Although this choice resulted in a generally reasonable compromise between a good acceptance rate and the loss of the quality of the data, a negligible acceptance rate was still noted when considering the shallowest layer (see Fig. 7). While we recognize that the magnitude of measurement errors is only assumed in our study, as no precise information on this aspect is available, we also note that such low acceptance rate is consistent with the observations that (i) one can argue that the turbulent behavior of flow at locations close to the sediment-water interface may influence the accuracy of the data at shallow depths more markedly than at larger depths, and (ii) data interpretation rests on a simplified diffusive model whose skill to represent the process may decrease close to the sediment-water interface. As an additional element, which might suggest that turbulence is related to highest measurement uncertainties at the shallowest layer, we observe that the need to assume the largest measurement errors at the shallowest layer was linked to all tests, independent from the sediment size. Therefore, the value of r C y was progressively increased until a reasonable acceptance rate was attained. A value r C y ¼ 45% is employed for the shallowest layers. While considering r C y \45% would provide a sufficient number of accepted diffusion coefficient values in Test 1 (see Fig. 7), for the sake of uniformity we decided to employ the same values of measurement error in all tests.
Non-zero values of posterior distributions associated with the deepest layers (Layers 4 and 5) and the smallest grain diameters tested (Tests 3 and 4) are found across a Fig. 5 Sample probability density and joint-relative frequency resulting from the stochastic inverse modeling approach for Test 1: a D 3 -D 2 , b D 4 -D 3 , c D 5 -D 4 d pdf D 5 . The red and magenta dots represent the values reported by Chandler et al. (2016) for the two replicates of this Test. We consider r Cy ¼ 25% for the calibration of D 5 , D 5 À D 4 , D 4 À D 3 , while r Cy ¼ 45% for the pair D 3 À D 2 significant portion of the support of the corresponding priors. This result could be due to the possibility that the combination of increased distance from the sediment-water interface and a small sediment size yields a reduced rate of the diffusion process. The occurrence of this phenomenon is consistent with (i) the modest reduction in time of solute concentration and (ii) the overlap observed by the experimental breakthrough curves at the lowest interfaces (Layer 3-Layer 4, Layer 4-Layer 5). Thus, diffusion coefficients at the deepest layers vary across ranges of similar magnitude. This, in turn, leads to an increase of the uncertainty associated with the estimated diffusion coefficients (as indicated in Fig. 10d by the large spread of accepted values across the support). Along the same lines, differences in the shape of the pdf of D 5 , as well as in terms of the number of accepted values associated with estimates of the pairs D iþ1 À D i , are observed between the various experimental  tests (see Figs. 5 and 10). For example, in Test 4, where the experimental breakthrough curves C 5 Ã ,C 4 Ã are essentially overlapped, the presence of multiple peaks in the pfd of D 5 (see Fig. 10 ((iii)-d)) can introduce elements of ambiguity in the interpretation, as opposed to the remaining. We further note that decreased concentration values are observed at shallow locations at early times, where the diffusion model employed may be inaccurate due to effects of marked spatial concentration gradients. Additionally as these points are located close to the sediment-water interface, transport involves small space-time scales, and may exhibit a non-Fickian pre-asymptotic behaviour. We observe a good agreement between the C i ðMAPÞ andC i in Test 1, the ensuing results being virtually indistinguishable for Layer 5 and displaying only minor differences for Layers 3 and 4. The trends displayed by the experimental observations are reasonably well reproduced by the modeling results, even as local oscillations in the concentration temporal gradients are not represented through the diffusive model. With reference to the latter point, the observed concentration C Ã 4 appears to display a steeper decrease than the corresponding model results, which display a reduced rate of convergence to the final (equilibrium) concentration (that is approximately equal to 20%). Experimentally observed concentrations within Layer 5 (i.e., C Ã 5 ) initially decrease faster than their model-based counterparts. Otherwise, we observe that the values C 5 ðMAPÞ are close to the experimental observations for long times (t [ 0:75 Â 10 5 s).
We recall that, considering a given pair of layers in the model calibration procedure of Nagaoka and Ohgaki (1990), the pattern of solute concentration at the top of the upper layer plays a significant role on the determination of the breakthrough curve at the interface between the layers. In other words, and referring to the bottom layer, the behavior displayed by C 5 ðMAPÞ calculated from Eq. (8) using as input values the diffusion coefficient D 4 ðMAPÞ, D 5 ðMAPÞ, and the experimental concentration C 4 Ã , is considerably affected by the temporal variation of the latter. The strong dependence between C 5 ðMAPÞ and C 4 Ã may contribute to the discrepancy observed between C 5 Ã and C 5 ðMAPÞ.
Results of corresponding quality are obtained for Test 2 and Test 3 (Fig. 6b,c). Estimated C 5 ðMAPÞ and C 3 ðMAPÞ values in Test 3 (Fig. 6c) display a temporal pattern which is in close agreement with their experimental counterpart (see in C Ã 5 Test 3 1 and in C Ã 3 Test 3 2, respectively). Otherwise, an increased discrepancy can be observed between the experimental data and their model-based counterparts for C 4 , the latter (i.e., C 4 ðMAPÞ andC 4 ) decreasing faster than the experimentally-based results and attaining a final concentration which is 20% lower than the experimental one.
With reference to Test 4 (Fig. 6d), we note that only C i ðMAPÞ values can be compared against the experimental data set, as estimated values for the effective diffusion coefficient are not documented by Chandler (2012), Chandler et al. (2016). The analytical curve C 5 ðMAPÞ exhibits a behavior which is similar to the experimental one, i.e., C Ã 5 Test4 2. A similar result is observed for the curve C 3 ðMAPÞ for t [ 10 5 s. On the other hand, model results corresponding to C 4 ðMAPÞ decrease following a trend significantly different from the one displayed by the experimental data. We observe that Test 4 displays similar late-time slopes for the experimental curves associated with C 3 Ã , C 4 Ã ,C 5 Ã , differing slopes being observed in the remaining tests. As mentioned above, estimation of the pair of coefficients D 2 À D 3 using concentration data C Ã 2 is characterized by a significantly reduced acceptance rate as compared to the remaining layers. As such, it has been subject to a specific analysis. Figure 7 depicts the results of the investigations corresponding to increasing values of r C y , to investigate the impact of measurement errors on the acceptance rate, the latter being quantified as the ratio between the number of accepted pairs D 2 À D 3 and the total number of realizations tested. As stated earlier, the rationale underlying this analysis corresponds to the observation that measurement (and model) error may increase as sampling locations are closer to the interface between the porous medium and the water. As shown in Fig. 7, the highest acceptance rate is associated with the experiments characterized by the largest grain size and bed shear velocity (i.e., Test 1), while the lowest acceptance rate corresponds to the scenario with the smallest grain size and lowest velocity (i.e., Test 4), the remaining two tests being positioned between these two extremes. This analysis suggests that the data quality is not uniform across the tests: parameter estimation in Test 1 appears to be associated with lower measurement errors than in the remaining settings, while data of Test 4 tend to be characterized by the lowest reliability (i.e., the largest estimated measurement error). This result is consistent with the observation that Chandler (2012), Chandler et al. (2016) did not report results associated with effective diffusion estimates in Test 4.
Our results indicate that acceptance rates increase with the grain diameter and with bed shear velocity, thus suggesting that the quality of the data and the ability of the considered diffusion model to interpret these are directly influenced by these two physical parameters. This result may also be related to the specific equipment employed during the tests. For example, the sensors used to measure concentrations are sensitive to temperature variations. As experiments characterized by small grain diameters are associated with a longer duration, the related measurements may display a time dependent error which is hard to model explicitly in the absence of further information. We observe that the overall performance of the model considered (i.e., Eqs. (7), (8)) is less satisfactory for the shallowest layers, a significant discrepancy between the temporal pattern of C 2 ðMAPÞ and its experimental counterpart (C Ã 2 ) being observed at Layer 2 (see Fig. 11 in Appendix B) in all tests.
We then investigate the influence of the grain size and the bed shear velocity on the diffusion process and assess if the signature of the exponential reduction of the diffusion coefficient with depth reported in the literature Chandler et al. (2016) is observed also through our results. We do so by normalizing the vertical coordinate below the watersediment interface through the characteristic sediment size (d g ), the diffusion coefficient being normalized by ðud g Þ. The results of the stochastic calibration process are then summarized in Fig. 8, where one can also compare the MAP parameter values resulting from the acceptance/rejection procedure against the corresponding estimates obtained by Chandler et al. (2016). We recall that examining two consecutive pairs of layers yields two distinct values for the diffusion coefficient of a given layer. As such, two values of the MAP (red squares) are displayed at each layer in Fig. 8. Fig. 8 Summary of the acceptance/rejection stochastic calibration results for all experimental tests: a Test 1, b Test 2, c Test 3, d Test 4, and the corresponding regression fits A graphical depiction of the acceptance range corresponding to the 25th and 75th percentile is included in Fig. 8 for each MAP value, thus providing quantitative information on the uncertainty about the estimated diffusion coefficient at each layer. Furthermore, our results allow quantifying the range of variability of accepted values for a given layer which stems from analyzing the results of two consecutive pairs of layers, i.e., an estimate of D i is obtained through the concentration histories C Ã i and C Ã iÀ1 . These two estimates are here obtained independently and results in Fig. 8 show that they identify distinct ranges for each D i . This result is possibly related to the different role played by the same coefficients in the two consecutive steps and suggests that the workflow proposed in Nagaoka and Ohgaki (1990) and depicted in Fig. 3 may not lead to optimal estimation results. In the following we consider for each coefficient D i only the probability distribution obtained using the data C Ã i , consistent with the interpretation given in Nagaoka and Ohgaki (1990).
We consider next the following relationship suggested by Chandler (2012) to characterize the dependence between (dimensionless) effective diffusion and depth Here, A and B are model parameters which we estimate by linear regression, upon substituting the values ofD i or D i ðMAPÞ in (10). The results of such an analysis are listed in Table 2. As shown in Fig. 8, our result display an exponential reduction of the diffusion coefficient with depth for Tests 1 and 2, i.e., for a grain size of d g ¼ 5 mm. Otherwise, the occurrence of a decreasing pattern is not consistent with the results obtained for Tests 3 and 4, i.e., for d g ¼ 0:625 mm. In Test 3, we observe a decay of the diffusion coefficient up to Layer 4, results obtained for Layer 5 being inconsistent with the trend found by Chandler et al. (2016). We observe that the effective diffusion coefficient displays an increasing trend with depth for y=d g \ À 100 in Test 4. This behavior is not supported by a well-defined physical interpretation. Note that, as mentioned above, concentration values observed at different depths for Test 4 display minimal differences (see Fig. 6). Similarly, concentrations observed in layers 4 and 5 in Test 3 display no appreciable decreasing trend (see Fig. 6d). These observations suggest that the available time series are too limited to enable a detailed assessment and quantification of the transport process taking place under these conditions and for y=d g \ À 100. Figure 9 juxtaposes the collection of values D i ðMAPÞ obtained for each experimental test analysed. We observe that the results obtained for Test 1 are similar to those of Tests 2. Otherwise, results obtained for Tests 3 and 4 follow a different trend. The results depicted in Fig. 9 suggest that recasting the problem in terms of dimensionless parameters of the kind included in Eq. (10) does not provide a unique interpretation to the available data and that the exponential trend may be representative solely of the grain diameters analyzed. While a physical explanation of the observed behavior is not clear at the current stage, the analysis of additional experimental data would be required to shed more light on this element.

Discussion and conclusions
Our study is aimed at the quantification of the uncertainty associated with effective diffusion coefficients employed to describe solute transport and mixing in the hyporheic region. We do so upon relying on a set of experiments performed at the laboratory scale considering solute transport close to a sediment-water interface (Chandler 2012;Chandler et al. 2016). These are associated with various combinations of (a) bed shear velocity, governing the intensity of turbulent fluxes propagating from the water body to the underlying porous bed, and (b) characteristic size of the grains forming the solid matrix of the porous medium. Similar to the analysis of Chandler (2012) and Chandler et al. (2016), we rely on the analytical formulation of Nagaoka and Ohgaki (1990) as our process model.
Our analysis then rests on a stochastic inverse modeling approach based on the acceptance/rejection algorithm. Our findings contribute to enhance the information content associated with the results in Chandler et al. (2016), as we provide a quantification of the uncertainty associated with estimated diffusion coefficients. Our study explicitly yields the posterior (i.e., conditional to observations) distribution Further to these elements, our procedure allows providing an appraisal of the quality of the experimental data, as quantified in terms of an experimental error considered for each test. As such, recognizing this unique feature of a stochastic approach of the kind we consider can contribute to favor the use of enhanced uncertainty propagation approaches, which are not always considered in hyporheic region studies. In the absence of more specific information, the acceptance rate can be viewed as a combined indicator measuring (i) data quality and (ii) the skill of the assumed diffusive model to interpret available observations. We observe that experiments associated with the smallest sediment diameter (i.e., d g = 0.625 mm in Tests 3 and 4) display (i) a low acceptance rate at shallow locations, and (ii) large uncertainty in the estimated diffusion coefficients far from the sediment/water interface.
Our results enable us to assess the uncertainty related to the exponentially decreasing trend of the effective diffusion with depth in the sediment bed, as embedded in Eq. (10). This type of trend has been proposed by Chandler et al. (2016) and is commonly assumed in current modeling of mixing processes within the hyporheic region. An exponential reduction of the diffusion coefficient with depth under the sediment-water interface is observed for most of the combinations of sediment size/bed shear velocity here analyzed. However, it is noted that the results obtained significantly differ between experiments performed with different sediment size. Our results show that the experiments performed with the largest sediment diameter (i.e., d g = 5 mm in Tests 1 and 2) can be characterized by diffusion coefficients that are well interpreted by Eq. (10). Contrary to what reported by Chandler et al. (2016), our analysis indicates an increase of the effective diffusion with depth in the experiments performed with the smaller grain diameter (i.e., d g = 0.625 mm in Tests 3 and 4). An increase of the diffusion coefficient is here documented at locations corresponding to depths larger than 100 grain diameters for the scenario with small-diameter sediments. This result might be related to a low accuracy of the experimental data associated with d g = 0.625 mm, as discussed above. This element is also recognized by Guymer (2020, personal communication), according to whom the data related to the tests with small sediment size might be affected by higher sensitivity to the variation of the room temperature, thus resulting in a possible reduction in the accuracy of the experimental results. Thus, it is still not possible to draw general conclusions which could be unambiguously related to the occurrence of small grain diameters.
We do recognize that the analytical formulation we employ is simple and there are opportunities for further advancements. At the same time, it is also apparent that having at our disposal enhanced datasets, eventually collected under a variety of experimental conditions, can contribute to identify processes which are only partly included, or eventually disregarded, within a given model. As such, our study shows that the currently available dataset does not allow generalizing our findings and formulating unambiguous conclusions on the contribution of quantities such as sediment size and/or bed shear velocity to solute exchange across and above the sediment-water interface. Further experimental observations, e.g. extended to longer temporal windows and with increased spatial detail, and an interpretive model characterized by less stringent assumptions than those proposed in Nagaoka and Ohgaki (1990) could be beneficial to improve the characterization of effective diffusion at larger depths. In addition to these elements, our study highlights the need for improved and rigorous quantifications of measurement uncertainties. These should then be fully considered to constrain estimates of model parameters as well as their estimation uncertainty. Acknowledgements The authors would like to thank Professor Ian Guymer (University of Sheffield) for providing experimental data and advices provided to third author on the hydro-physical process occurs in the water-riverbed interface.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.

Appendices
Calibration of the diffusion coefficient at the Layer 3, Layer 4 and Layer 5 for Test 2, Test 3 and Test 4 Sample probability density and joint-relative frequency resulting from the stochastic inverse modeling approach in Test 2, Test 3, Test 4 We display here for completeness the relative frequency distributions of the diffusion coefficients obtained for Tests 2, 3, 4. Results for Test 1 are illustrated in Fig. 5. (Fig. 10) Fig. 10 Sample probability density and joint-relative frequency resulting from the stochastic inverse modeling approach for Test 2 (i), Test 3 (ii) and Test 4 (iii): a D 3 À D 2 , b D 4 À D 3 , c D 5 À D 4 d pdf D 5 . The red and magenta dots represent the values reported by Chandler et al. (2016) for the two replicates of each Test. We consider r Cy ¼ 25% for the calibration of D 5 , D 5 À D 4 , D 4 À D 3 , while r Cy ¼ 45% for the pair D 3 À D 2 Calibration of the diffusion coefficient at the pair Layer 2-Layer 3 Range of accepted solute concentration at the interface Layer 2-Layer 3 for Test 1, Test 2, Test 3, Test 4 The Appendix shows the evolution in time of the concentration C 2 at the interface Layer 2-Layer 3 corresponding to the accepted values of the diffusion coefficient D 2 and D 3 for Test 1, 2, 3, 4. Additional results to the ones displayed in Fig. 6. (Fig. 11) Fig. 11 Evolution in time of the concentration at the interface Layer 2-Layer 3: Test 1 a, Test 2 b, Test 3 c, Test 4 d. The curves compare the analytical solution calculated with MAP and the one computed with the diffusion coefficients resulting from the calibration procedure by Chandler et al. (2016), the shaded areas indicate the envelope of the realizations obtained with accepted values of the diffusion coefficient