3.1 Numerical model

We use the mHM-OGS coupled model, proposed by Jing et al. (2018a), to simulate terrestrial hydrological processes. This coupled model was developed for extending the predictive capability of mHM from land surface processes to the subsurface flow and transport processes. Specifically, mHM (Kumar et al., 2013; Samaniego et al., 2010) is used to partition water budget components, while the OGS porous media simulator (Kolditz et al., 2012) is used to compute groundwater flow and transport processes by using mHM-generated recharge as driving forces. For details on the mHM-OGS coupled model, please refer to Jing et al. (2018a).

The catchment water storage is conceptually partitioned into soil zone storage and deep groundwater storage; the two corresponding components are computed by mHM and OGS, respectively. The soil-zone dynamics of TTD has been well studied using mHM in a previous work (Heße et al., 2017). Hence, in this paper, we perform explicit forward modeling of the saturated-zone TTD through a 3-D OGS groundwater model by using the mHM-generated recharge as the external forcing.

In this study, we focus on the travel times in the saturated zone. Saturated groundwater flow is characterized by the continuity equation and Darcy's law:

where S is the specific storage coefficient in confined aquifers, or the specific yield in unconfined aquifers (1 / L); ψ_p is the pressure head in the porous medium (L); t is the time (T); q is the specific discharge or the Darcy velocity (LT⁻¹); q_s is the general source/sink term (T⁻¹); K_s is the saturated hydraulic conductivity tensor (LT⁻¹); and z is the vertical coordinate (L).

We use the RWPT method to track the particle movement. The RWPT method is embedded in the source codes of OGS (Kolditz et al., 2012; Park et al., 2008). Derived from stochastic physics, RWPT is under the assumption that the advection process is deterministic and the diffusion–dispersion process is stochastic. The theoretical background of the RWPT method is described in detail in Appendix A.

3.2 Numerical model setup

3.2.2 Modeling procedures

Table 1Adjustable ranges of the hydraulic parameters.

Download Print Version | Download XLSX

Figure 2Recharge realizations used in this study (unit: mm). They were sampled from a high-resolution dataset of land surface fluxes for Germany (Zink et al., 2017).

Download

The numerical experiment to explore the uncertainty of TTDs is performed through the following workflow.

1. Eight spatially distributed recharge realizations are sampled from a high-resolution dataset of land surface fluxes for Germany, in which mHM is used to simulate land surface hydrological processes. The details of the dataset and the sampling method are described in the following section.

2. For each recharge realization, a series of equally probable realizations of K_s fields are generated using the null-space Monte Carlo (NSMC) method.

   The NSMC method takes advantage of the hybrid Tikhonov–TSVD (truncated singular value decomposition) method in the PEST parameter estimation code to produce Monte Carlo realizations of parameters (Doherty and Hunt, 2010; Tonkin and Doherty, 2009). This approach is able to efficiently generate an ensemble of parameter fields that are conditioned to expert knowledge and measurements. Here, the observations of groundwater levels from 18 spatially distributed monitoring wells are used to calibrate the OGS model (the locations of monitoring wells are illustrated in Fig. 1a). Before generating parameter sets, we calibrate the model to obtain the best-fit hydraulic conductivities, as well as a covariance matrix of the parameter probability distributions. On the basis of this information, many fully distributed K_s fields are randomly generated from a uniform distribution of hydraulic conductivity values (Doherty, 2015). As shown in Table 1, the range of hydraulic conductivities is predefined based on values obtained from a geological survey (Wechsung et al., 2008). As a result, a total of 400 parameter sets conditioned on both observations and reality are generated for the uncertainty analysis.

3. In each parameter realization, a large number of particles are injected through the top surface of the groundwater model. The spatial density of particles is proportional to the spatially distributed recharge rates.

   To accurately interpret the travel time distribution, a large number of particles (e.g., approximately 80 000 particles in the case study) is released into the top surface of the groundwater model. The released particles serve as samples of water parcels for deriving their travel time distributions. In doing so, the density of particles is set proportional to the recharge at the corresponding grid cell (Fig. 3). Each particle tracer represents a volumetric recharge rate of around 700 m³ year⁻¹.

4. An ensemble of forward simulations using the RWPT method is performed over all realizations of K_s fields.

   In each realization of the ensemble parameter sets, forward simulations of particle tracking are performed. In this study, we focus on the predictive uncertainty within the convection process. Therefore, the molecular diffusion coefficients are universally set to 0 for all ensemble simulations. The porosity of the study domain is set to 0.2 universally. Through the above procedures, the flow paths and the corresponding residence times can be fully traced in the model at random times and locations, facilitating the detailed characterization of TTDs.

In parallel to this analysis, a sensitivity analysis for the spatial variability of recharge is also performed. Two different recharge scenarios are compared for this purpose: (1) the spatially distributed recharge generated by mHM and (2) the uniform recharge that is equal to the spatial average of the distributed recharge. Other parameters, including the porosity and the hydraulic conductivity, remain identical in these two recharge scenarios.

Figure 3Two different spatial distributions of particle tracers for the RWPT method. (a) The mass-weighted distribution of particles based on the recharge estimated by mHM. This is the default spatial pattern of particle tracers in this study. (b) The uniformly distributed particle tracers used in the uniform recharge scenario.

Download

3.2.5 Parameter uncertainty

Figure 4Box-plot of stochastically generated hydraulic conductivity (K_s) for each geological layer in eight recharge realizations. Note that the parameters mo2, mm2, and mu2 are not shown in this figure because the less permeable subunits of the Muschelkalk (mo2, mm2, and mu2) are tied with the respective more permeable subunits (mo1, mm1, and mu1) with a factor of 0.1.

Download

Figure 5Observed and simulated groundwater heads for each parameter and recharge realization. The results of 400 realizations (R1K1–R8K50) are categorized by recharge realization and shown in different panels.

Download

Multiple calibration-constrained K_s fields were stochastically generated for each recharge realization. Figure 4 shows the box-plot of generated hydraulic conductivities in all realizations categorized by geological unit. The hydraulic conductivity of the Lower Muschelkalk (mu) has the highest uncertainty (10⁻⁸–10⁻⁵ m s⁻¹) because the conductivity of mu is insensitive to groundwater head observations, given that it is the deepest geological layer and that no monitoring well is located in this layer (Table B1). The other parameters fluctuate moderately and are constrained within 1 order of magnitude in most of the recharge realizations. Hydraulic conductivities of several permeable layers (mo, mm, alluvium, and soil) increase from R1 to R8, which is not surprising because the hydraulic conductivity increases with increasing recharge and constant groundwater head. Moreover, the hydraulic conductivities of the above layers are roughly linearly correlated with the corresponding recharge in each recharge realization. Figure 5 shows the simulated and observed groundwater heads for all 400 realizations. All of the 400 realizations are well constrained to observations, with the root mean square error (RMSE) of groundwater level residuals being lower than 4.6 m in all of the considered recharge realizations.

3.3 Theory of analytical StorAge Selection function

The travel time is defined as the time spent by a moving element (either a water particle or a solute) in a control volume of a hydrologic system. In principle, the control volume can be defined at arbitrary spatial scales (i.e., from the molecular scale to the regional scale). Considering a hydrologic system in which the input flux (J) and the output fluxes ($Q_{\mathrm{1}},Q_{\mathrm{2}},\mathrm{\dots },Q_n$) are known, each parcel of water within the system is tagged using its current age τ. The age-ranked storage $S_T=S_T(T,t)$ is defined as the mass of water in the system with age τ<T. The backward form of the master equation (ME) for TTD in a control volume can be expressed as follows (Botter et al., 2011; Harman, 2015; Van Der Velde et al., 2012):

with boundary condition $S_T(\mathrm{0},t)=\mathrm{0}$, where T is the residence time of the oldest water parcel in storage S_T; t is the chronological time; ${\overleftarrow{P}}_{Q_j}(T,t)$ is the cdf of the backward travel time distribution of output flux Q_j; J(t) is the input flux at time t; and Q_j(t) is the output flux at time t. Specifically in this study, J is the groundwater recharge, and Q is composed of two components: the stream baseflow and the abstraction at production wells.

The SAS function describes the fraction of water parcels leaving the control volume at time t, which is selected from the age-ranked storage S_T. Following the above definition, the SAS function can be linked with the backward travel time distribution ${\overleftarrow{P}}_Q(T,t)$ (Harman, 2015)

for $S_T=S_T(T,t)$. Ω_Q is the cumulative form of the SAS function.

Three instances of SAS functions using gamma distribution are shown in Fig. 9a. In case the age distribution of each outflow is uniformly selected from all water storages with various ages, the outflux TTDs turn into a random sample of the storage residence time distribution (RTD). The random sampling (RS) is equivalent to the uniform SAS function (Fig. 9a). Many past studies have also considered the random sampling as a proper description of the sampling behavior for heterogeneous catchments (Benettin et al., 2015; Cartwright and Morgenstern, 2015; Heße et al., 2017). Equation (4) in this case has the analytical solution (Danesh-Yazdi et al., 2018; Harman, 2015)

where p_S(T,t) is the pdf of the residence time distribution and S(t) is the storage at time t. Specifically, in the case of a steady-state hydrodynamic system, Eq. (5) is further simplified into an exponential form:

Equation (6) is the analytical solution of backward TTD under the RS assumption.

In the idealized saturated groundwater aquifer, Eq. (6) is equivalent to the analytical solution derived by Haitjema (1995). Based on the Dupuit–Forcheimer assumption, Haitjema (1995) derived a formula about the frequency distribution of residence time:

provided that nH∕J is constant over the entire domain, the recharge is spatially uniform, and the aquifer is locally homogeneous, where n is the porosity, H is the saturated aquifer thickness, and $\stackrel{\mathrm{‾}}{T}$ is the weighted mean travel time in the aquifer.

3.4 Linking the SAS functions to the physically based numerical model

Danesh-Yazdi et al. (2018) developed an approach to link the analytical SAS functions to the fully distributed numerical model. Although differing in numerical model and particle tracking scheme, we apply the same approach to link the SAS functions and the numerical model in this study. Equation (3) under the steady-state condition can be further simplified as

By combining Eqs. (4) and (9), the age-ranked storage S_T can be calculated directly under the steady-state assumption:

Combining Eq. (10) with Eq. (4), the SAS function can be directly computed from the simulated backward TTD using the numerical model.

3.5 Predictive uncertainty of TTDs

The theoretical framework of predictive uncertainty in this paper is based on Doherty (2015). As indicated in Bayes' theorem, the parameters of a model retain uncertainty, given that they have been adjusted to best-fit values achieved during calibration. Nevertheless, the uncertainty of parameters is subject to constraints. One of the constraints resides in the fixed adjustable range of parameters, in which expert knowledge must be respected. Another constraint is exerted by the parameterization process.

While the computationally expensive Bayesian approach offers a complete theoretical framework for predictive uncertainty evaluation, practical modeling efforts are often based on model calibration and a subsequent analysis of error or uncertainty in post-calibration predictions (Doherty, 2015). Ideally, the best-fit parameters achieved through calibration can reduce the predictive error to a minimum, with the minimum predictive error being the inherent uncertainty. However, the best-fit parameter is always biased from the true parameter because the essential imperfection of the model may facilitate or hamper the achievement of the minimum. Therefore, the motivation of uncertainty analysis in this study is to quantify and minimize the predictive uncertainty of travel time distributions, given that the parameters are plausible and that the model can reproduce well the groundwater heads.

4.1 Uncertainty of TTD predictions

Figure 63-D view of flow pathlines of some particles in realization R5K1. Note that only a limited number of particle pathlines are displayed here.

Download

Flow paths of particle tracers in a random parameter realization (R5K1) are displayed in Fig. 6, serving as a visual reference for the regional groundwater flow pattern and the residence time distributions. In this realization, the deep low-permeability geological layers (e.g., mm2 and mu2) act as low-permeability aquitards. Therefore, the majority of streamlines do not enter these geological layers (Fig. 6).

Figure 7 displays both the simulated TTDs using 50 K_s fields for eight recharge realizations (orange solid lines) and the reference TTDs, represented by fit blue dash-dot curves using the exponential model (Eq. 6). The ensemble average (μ) and the coefficient of variation (c_v) of MTTs for each recharge realization are also calculated and shown in Fig. 7. Note that if the number of parameter realizations is sufficiently large, the ensemble average of MTTs (μ) will converge to the simulation result using the best-fit parameters achieved through model calibration (Doherty, 2015). Noticeable variability of TTDs can be observed with respect to different recharge realizations. Generally, the μ values show a decreasing trend from 166.5 years in recharge realization R1 to 110.9 years in recharge realization R8, with only two exceptions (R3 and R6), which is not surprising based on the inversely linear dependency between recharge J and μ, indicated by Eq. (8). In each recharge realization, the different realizations of K_s fields manipulate the mean travel time. c_v varies from 7.81 % (for R5) to 15.56 % (for R3), indicating a modest degree of uncertainty propagated from K_s estimation to TTD prediction.

The exponential model under the RS assumption is fit to the ensemble-averaged TTD of numerical solutions (see the black lines in Fig. 7) using Eq. (6). As shown in Fig. 7, the shapes of numerically simulated TTDs significantly deviate from the exponential distribution under the RS assumption, indicating a nonuniform sampling behavior of different water ages. The TTDs of numerical simulations are more right-skewed than the analytical TTDs under the RS assumption. This phenomenon reveals that the catchment TTD cannot be replicated by the single random sampling store.

Based on Eq. (8), we can approximate the “effective volume” of water involved in the transport process in the aquifer. The effective volume of groundwater storage related to the transport process is calculated and shown in Table 2. The effective volumes of storage (S^eff) estimated by the numerical solutions range from 9.8 to 12.0 m, whereas the S^eff estimated by the analytical solution range from 6.8 to 7.5 m. The groundwater storage that contributes to the streamflow is significantly smaller than the total groundwater storage (48.3 m). This difference is because most of the released particles only exist in the upper permeable layers rather than being spread evenly over the whole aquifer–aquitard system. We are aware that this is only a first-order approximation because the analytical solution is only rigorously valid for the idealized homogeneous aquifer system (Haitjema, 1995).

Moreover, we assess the propagation to the MTT predictions from input and parameter uncertainty yielded by the eight recharge realizations and the Monte Carlo realizations of hydraulic conductivities. Figure 8a depicts the distribution of MTTs of the ensemble simulations. The MTTs of the 400 realizations range from 87 to 212 years. At the same time, the ensemble average of MTTs over all realizations of recharges and K_s fields is 135.1 years, and the coefficient of variation is 18.93 %. Figure 8b depicts the relationship between the ensemble average of MTTs and the spatially averaged recharge rates. We observe a roughly inversely proportional relationship between the ensemble average of MTTs and the spatially averaged recharge rates (Fig. 8b). The standard deviations (σ) of simulated MTTs range from 12.9 years (for R6) to 24.7 years (for R3).

Figure 7Travel time distributions of ensemble simulations and analytical solutions categorized by recharge realization. The orange lines show the simulated TTDs of all realizations of K_s fields for each recharge realization. The black lines denote the ensemble-averaged TTDs of each recharge realization. The blue dash-dot line is the fit analytical TTD under the RS assumption. The analytical MTT, the mean (μ), and the coefficient of variance (c_v) of the simulated MTTs are also shown in this figure.

Download

Table 2Effective groundwater storages related to the transport process for each recharge realization.

Download Print Version | Download XLSX

Figure 8Uncertainty quantification: Monte Carlo simulations of MTT predictions categorized by recharge realization. (a) shows a histogram of MTT predictions. (b) shows the relationship between the ensemble-averaged MTT (μ) and the spatially averaged recharge. The error bars represent the standard deviation of MTTs (σ) for each recharge realization.

Download

4.2 Uncertainty of young/old water preference

Figure 9Cumulative rank SAS functions as a function of normalized age-ranked storage. (a) Schematic of cumulative rank SAS functions parameterized by gamma distribution with the shape parameter a=0.5, 1, and 2. (b) Cumulative rank SAS functions of the ensemble simulations (light grey lines) and the ensemble average for each recharge realization.

Download

Figure 9a provides an intuitive illustration of the relationship between the cumulative rank SAS functions and the systematic preference for discharging water with different ages. Figure 9b shows the cumulative rank SAS functions of all ensemble simulations (obtained from 400 realizations of K_s fields in eight recharge scenarios). The figure is categorized into eight groups by different colors and line styles, each representing a recharge realization. Generally, the system has a weak preference to select younger water as discharge, despite different recharge realizations and K_s realizations. Nevertheless, a moderate variation in SAS functions for the ensemble simulations can be observed. The variation of interest is introduced by the spatial distribution and velocity of flow pathlines that are controlled by different K_s fields. For example, a more permeable shallow aquifer layer will activate more flow pathways in this layer and thus will introduce a stronger preference for young water. In particular, the significant variation in hydraulic conductivities in the deepest geological layer (e.g., Lower Muschelkalk) has a pronounced impact on the discharge of old water. With a thickness of up to 100 m, the hydraulic conductivity of this layer controls how many water parcels can enter into this layer and how deep the flow paths can develop. This effect can be evidenced by a large difference in the SAS functions related to old ages and a relatively smaller difference in those related to young ages in Fig. 9b.

4.3 Sensitivity to the spatial pattern of recharge

Figure 10Sensitivity of (a) TTDs and MTTs and (b) SAS functions to the spatial pattern of recharge.

Download

Figure 10a depicts the sensitivity of simulated TTDs and MTTs to the spatial distribution of recharge, while Fig. 10b shows the sensitivity of the cumulative SAS function to the spatial pattern of recharge. The reference simulation is set up using the spatially uniform recharge that is equal to the spatial average of the spatially distributed recharge, while all of the other parameters in these two simulations are held identical. The different spatial distributions of recharge have a clear effect on the shape of TTDs. It appears that the most evident difference between the TTDs of the two recharge scenarios occurs at the early period. Additionally, the simulated MTT using the uniform recharge appears to be smaller than that using the spatially distributed recharge. Figure 10b indicates that the simulation using uniform recharge has a consistently stronger preference for sampling young water than the simulation using spatially distributed recharge. Nevertheless, both scenarios show a general preference for young water.

The difference in TTDs and SAS functions is not induced by the variability in internal hydraulic properties, since the two simulations share the same K_s field. Rather, it is mainly induced by the different flow paths of particle tracers in the two recharge scenarios. The spatially distributed recharge simulated by mHM reveals that the upstream mountainous area has higher recharge rates than those in the lowland plain. By contrast, the uniform recharge scenario neglects this spatial nonuniformity. This difference results in the following: (a) under a uniform recharge scenario, more particle tracers enter the system from locations near the streams at lowland plains (Fig. 3b), indicating that more particle tracers are transported in the local flow system rather than in the regional flow system (Toth, 1963); and (b) higher recharge rates at lowland plains accelerate the particles' movement in this area and shorten their travel times. As such, local particle flow paths within the shallow aquifer layer at lowland plains are activated, leading to a stronger preference for sampling local flow paths in shallow aquifer layers and therefore a stronger preference for young ages. Our findings are in line with the observations by Kaandorp et al. (2018), wherein the authors found a relatively higher preference for the selection of older water in the upstream area than that in the downstream area of the Springendalse Beek catchment.

5.1 Uncertainty of external forcing, internal property, and TTD predictions

In the idealized aquifers where groundwater flow is of the Dupuit–Forchheimer type, the recharge is uniform and the aquifer is locally homogeneous. TTD is controlled by recharge, saturated aquifer thickness, and porosity, and it is independent of hydraulic conductivity (Haitjema, 1995). In a real-world catchment with complex geometry and topography, stratigraphic aquifers, and nonuniform recharge, our numerical exploration demonstrates that the groundwater TTD is dependent on both recharge and the hydrostratigraphic K_s field.

The mechanisms behind the effects of the recharge rate and the K_s field are different. Given that the spatial pattern of recharge remains the same, a higher recharge rate intensifies the fluxes for the complete range of flow pathways. This process forces water downward in the recharge zones and upward to the discharge area. As such, flow rates through all flow pathways are increased equally, and the spatial distribution of flow pathways is not changed. The corresponding SAS function is also not changed (Botter et al., 2010; Cartwright and Morgenstern, 2015; Kaandorp et al., 2018). In contrast, a different K_s field activates flow pathways in more permeable layers, deactivates flow pathways in less permeable layers, changes the spatial distribution of flow pathways, and therefore changes the shape of the SAS function (Harman et al., 2016; Kim et al., 2016).

We also underline the value of observational data in reducing predictive uncertainty in simulated TTDs. In this study, the majority of model parameters can be adequately conditioned by spatially distributed groundwater head observations (Fig. 4). Provided that most of the hydraulic conductivities are constrained to the model-to-measurement misfit and reality, the TTD predictions can also be effectively bounded. This is evidenced by Fig. 7, from which moderate values of c_v ranging from 7.81 % to 15.56 % in different recharge realizations can be observed. The ensemble-averaged MTTs for different recharge realizations also have a high c_v (15.70 %), implying that the TTD prediction appears highly sensitive to recharge. Our findings are in line with Danesh-Yazdi et al. (2018), in which the interplay between recharge and subsurface heterogeneity was investigated and a strong dependency of TTDs on the recharge was observed.

5.2 Analytical model and SAS functions

The analytical solution of TTD, assuming a random sampling of water, cannot properly replicate the TTD of numerical simulation in the study domain. In the stratigraphic aquifer with complex topography and diffuse recharge, the analytical solution using Eq. (6) underestimates the MTT. Note that this conclusion holds when the simulated TTD has a relatively larger long-tail behavior than the exponential distribution. Such observations have also been reported for other real-world aquifers (Basu et al., 2012; Eberts et al., 2012; Kaandorp et al., 2018). This finding can be seen as an extension of Basu et al. (2012), in which the authors found a moderate discrepancy between the analytical solution using Eq. (6) and the numerical solution in a small catchment.

It is obvious that the analytical solution under the RS assumption cannot explicitly include the impact of the distributed hydraulic properties of stratigraphic aquifers and the spatially nonuniform recharge. The above limitations of analytical models may introduce a significant predictive error for the TTD predictions, as shown in Fig. 7. Moreover, the effective volume of storage related to the transport process S^eff is calculated utilizing the numerical solutions of distributed models. As a first-order approximation of effective storage volume, it suggests that only a small fraction of total groundwater volume is involved in the active hydrologic cycle.

The SAS function provides a good interpretability for simulation results using the fully distributed model in terms of characterizing the preference for releasing water of different ages. We find that the SAS functions are weakly dependent on the K_s fields in the stratigraphic aquifer system, but the overall preference for young water does not change. This weak dependency can be explained by the fact that different realizations of K_s fields modify the spatial distribution of particle flow paths. The overall tendency for young groundwater in the saturated aquifer has been reported by many past studies (Danesh-Yazdi et al., 2018; Kaandorp et al., 2018; Yang et al., 2018). Our study links the explicit simulations of travel times and the analytical SAS functions and offers original insights into the uncertainty propagated from recharge and the K_s fields to the SAS functions.

5.3 Dependency of TTDs and SAS functions on the spatial pattern of forcings

The sensitivity of the TTDs and SAS functions to the spatial pattern of recharge forcings can be explained mainly by the different flow paths of particle tracers, resulting primarily from the spatially heterogeneous fields of recharge across the study catchment. For the regional groundwater system, the spatial variation of recharge determines the distribution of starting points of the flow pathlines of tracer particles. For example, more particles will be injected from recharge zones that are typically located in high-elevation regions, resulting in a higher weight of flowlines starting from high-elevation regions. The pronounced spatial variability of recharge also controls the systematic (water age) preference for particles existing from the system (to river discharge) that originated from different regions and therefore exerts a strong control on the shape of the SAS function.

In the study catchment, an oversimplified spatially uniform recharge results in a smaller MTT and a stronger preference for discharging young water compared to those taking the spatial variability of recharge. Such observations are conditioned to site-specific features of the study catchment, which is noticed only when the following apply: (a) a site is located in a headwater catchment under a humid climate condition; (b) the recharge rate in areas close to the drainage network is generally lower than that in areas far away from the drainage network; and (c) the system is under (near) natural conditions, meaning that artificial drainage and pumping do not dominate the groundwater budget.

The assumption of spatially uniform input forcing has been widely applied in regional-scale subsurface hydrologic models (Yang et al., 2018; Zghibi et al., 2015). Our study indicates that a reasonable characterization of spatial patterns of recharge is crucial for reliable TTD prediction. Unfortunately, it is quite challenging to confidently quantify the groundwater recharge at the regional scale under today's technique due to the lack of reliable measurements (Cheng et al., 2017; Healy and Scanlon, 2010; Zink et al., 2017). Appropriate techniques should be chosen to estimate groundwater recharge according to the study goals and the spatial and temporal scales.

5.4 Implications for the applied groundwater modeling

Uncertainty limits the applicability of groundwater models. Most of the applied groundwater models are deterministic models that use direct values of inputs and parameters instead of probabilistic distributions of them. Specifically, both the model inputs and the inversion process are deterministic, leading to a deterministic best-fit parameter set achieved during model calibration. Our study reveals limitations of the above modeling procedure and suggests that the probabilistic distribution of inputs and parameters should be considered for the applied modeling. The main limitation is that the single exclusive assignment of recharge is inadequate for the simulation of transport processes because the error can be propagated from inputs to the conditioned parameters (e.g., hydraulic conductivities) through the calibration process (Fig. 4). This accumulated error will further lead to a seriously biased prediction of travel times. Additionally, the modeling workflow used in this study is computationally more efficient than the Bayesian approach and is suitable for complex real-world applications.

The degree of predictive uncertainty is highly dependent on the parameterization scheme. Some highly parameterized models are potentially ill posed due to the paucity of data and therefore cannot be constrained by the available calibration dataset. In this case, the predictive uncertainty of TTDs is potentially high (Danesh-Yazdi et al., 2018; Weissmann et al., 2002). Stratigraphic aquifer models with zoned parameters are still widely used for applied groundwater modeling because the field representation of local-scale heterogeneity is difficult. Given that the aquifer model is stratigraphic and the number of parameters is less than the number of observations, most of the adjustable parameters can be effectively bounded. In this case, the uncertainty of input data (e.g., recharge) appears to have a primary influence on the TTD predictions. Note that here we do not account for the error caused by model structural deficiency.