Quantifying time-variant travel time distribution by multi-fidelity model in hillslope under nonstationary hydrologic conditions

The travel time distribution (TTD) is a lumped representation of groundwater discharge and solute export responding to rainfall. It reflects the mixing process of water parcels and solute particles of different ages and characterizes reactive transport progress in hillslope aquifers. As a result of the mixing process, groundwater leaving the system at a certain time is an integration of multiple water parcels of different ages from different historical rainfall events. Under nonstationary rainfall input condition, the TTD varies with transit groundwater flow, leading to the time-variant TTD. Most methods for estimating time5 variant TTD are constrained by requiring either the long-term continuous hydrogeochemical data or the intensive computations. This study introduces a multi-fidelity model to overcome these limitations and evaluate time-variant TTD numerically. In this multi-fidelity model, groundwater age distribution model is taken as the high-fidelity model, and particle tracking model without random walk is taken as the low-fidelity model. Non-parametric regression by non-linear Gaussian process is applied to correlate the two models and then build up the multi-fidelity model. The advantage of the multi-fidelity model is that it 10 combines the accuracy of high-fidelity model and the computational efficiency of low-fidelity model. Moreover, in groundwater and solute transport model with low Péclet number, as the spatial scale of the model increases, the number of particles required for multi-fidelity model is reduced significantly compared to random walk particle tracking model. The correlation between high and low-fidelity models is demonstrated in a one dimensional pulse injection case. In a two dimensional hypothetical model, convergence analysis indicates that the multi-fidelity model converges well when increasing the number of high-fidelity 15 models. Error analysis also confirms the good performance of the multi-fidelity model.


Introduction
The travel time of water in a hillslope is defined as the time duration for a water parcel to travel from one location where rainfall enters the aquifer to another location where groundwater exits the aquifer. Distinct water parcels can have very different travel time due to the mixing processes in the hillslope. The travel time distribution (TTD) is a lumped probability density function 20 that represents the probability distribution of travel time. If water parcels are sampled from hillslope outflow, the TTD of the water parcels is defined as the backward TTD (bTTD); if the water parcels are tracked from rainfall input, the TTD of the water parcels is defined as the forward TTD (fTTD). TTD encapsulates the hydrologic dynamic and solute transport processes into a mathematic function and helps to study the time-varying transport and reaction processes in the hillslope, such as the genesis and transport of weathering products (Ameli et al., 2017;Maher, 2010Maher, , 2011 and the attenuation of contaminants (Queloz et al., 2015).
Based on the time-invariant assumption of TTD, many previous studies approximated the TTD by the Gamma (Γ) distribution or some other probability density functions (Hrachowitz et al., 2009(Hrachowitz et al., , 2010Kirchner et al., 2000;McGuire and McDonnell, 2006). The time-invariant TTD quantifies the yearly averaged behaviors of transport processes. However, when attempting to explore the time variability of water and solute transport in hillslopes under nonstationary rainfall input condition, it is nec-30 essary to adopt the time-variant assumption of TTD. Hence, several theories and approaches have been developed in the last decade to study the time-variant TTD.
The age master equation in accompany with the StorAge Selection (SAS) function (Botter et al., 2011) is one of the fundamental theories for these approaches. The age master equation is a storage mass balance equation characterizing the evolution of storage age, namely the residence time. In the age master equation, the backward TTD quantifies the loss of age-tagged stor- 35 age, and the SAS function relates the RTD with backward TTD. The SAS function illustrates how water parcels of different ages are released from storage. It characterizes the selection preference during the outflow generation. There are three principal types of selection preference: selection preference for younger water, selection preference for older water and no selection preference. In transient flow systems, the time-variant SAS function indicates the internal variability of the hydrologic system, such as the temporal variability of flow paths. If the SAS function and storage volume are known based on prior information, 40 the time-variant TTDs can be calculated by solving the age master equation analytically (Botter et al., 2011). However, in most real cases, the SAS function is unknown and is inferred from observed tracers using the age master equation.
A common approach to calculate time-variant TTDs is by coupling the age master equation with observed tracer data (Benettin et al., 2015a(Benettin et al., , b, c, 2017(Benettin et al., , 2013van der Velde et al., 2015). This approach requires long-term continuous records of rainfall rate, discharge rate and concentrations of the conservative tracer in both rainfall and discharge water. The approach 45 generally employs an inverse modelling scheme to infer parameters in the adopted SAS function and age master equation.
The long-term time series data are used to calibrate the model. When the posterior distributions of model parameters are estimated, time-variant TTD can be derived immediately from the age master equation. However, in many real catchments and hillslopes, the long-term monitored solute concentrations are not available. This approach of solving the age master equation mainly requires the information at the recharge and discharge boundaries and does not need to have extensive hydrogeological 50 information inside the hillslope, so the approach cannot be used to study the internal factors controlling the time variation of TTD and storage selection process, such as the influence of permeability architecture and preferential flow.
In addition to the approach to infer TTD based on age master equation, (Harman and Kim, 2014) developed an approach to derive the time-variant TTD directly from the breakthrough curve (BTC) of tracers under periodic hydrologic condition. No pre-assumption about the SAS function is needed in this approach. The time-variant TTD can be constructed by decomposing 55 the overprinted BTC into the contributions of multiple input events. This approach has been applied to indoor experimental hillslope successfully to explore the transport variability and storage-dependent TTDs (Kim et al., 2016;Pangle et al., 2017).
The periodic condition is easy to control in experimental hillslopes, but such a condition is rare in real-world hillslopes.
Random walk particle tracking is another approach to simulate the time-variant travel time distribution. Given the hydrogeological properties and boundary conditions in hillslope, the flow field can be determined by solving the governing equation 60 of groundwater flow. A set of particles are released at the input boundary and are tracked thereafter. The movement of each particle depends on the transient flow field and molecular diffusion and mechanical dispersion. The molecular diffusion and mechanical dispersion are approximated by random walk scheme. Several studies utilized this approach to analyze the timevariant travel time distributions and compared the results of the SAS function from numerical simulations with analytical SAS functions (Ameli et al., 2016a(Ameli et al., , b, 2017Danesh-Yazdi et al., 2018;Kaandorp et al., 2018;Yang et al., 2018). This approach 65 makes it convenient to study the relation between the internal variability of transport processes and time-variant TTD, such as the variability of flow path and preferential flow. However, if the transport is strongly influenced by molecular diffusion and mechanical dispersion indicating a low Péclet number, this approach requires releasing a large number of particles during each input event to reduce the random noise when implementing the random walk scheme.
Time-variant travel time distribution can also be evaluated by groundwater age distribution model. This approach adds an 70 age dimension to the governing equation of groundwater and solute transport (Benettin et al., 2013;Ginn et al., 2009;Gomez and Wilson, 2013;Massoudieh and Ginn, 2011;Massoudieh, 2013). In a three dimensional (3D) transient flow system, it is challenging to solve the five dimensional (5D) partial differential equation to obtain the groundwater age distribution. The five dimensions consist of three spatial dimensions, one age dimension and one time dimension. Even so, there are some methods to solve this 5D equation. Ginn et al. (2009) solved this equation numerically by directly discretizing this equation in a 5D space, 75 but most current softwares or codes do not have the module to visualize and discretize the 5D space, let alone solve the 5D equation. Gomez and Wilson (2013) solved this equation by tracking each rainfall input which requires solving breakthrough curves for each rainfall event. Cornaton (2012) solved this equation numerically in a Laplace space. The advantage of this groundwater age distribution model is that it is valid to depict advection-dispersion transport processes with any Péclet number.
In this study, a multi-fidelity model is proposed to quantify the time-variant TTD. It consists of two models with different 80 levels of fidelities. The high-fidelity model is the groundwater age distribution model, and the low-fidelity model is the particle tracking model without considering random walk. Nonlinear Gaussian process regression method is applied to relate the highfidelity with low-fidelity models. Then, the 5D governing equation of groundwater and solute transport can be solved and timevariant TTD can be calculated by the multi-fidelity model. The multi-fidelity model takes advantages of both the efficiency of low-fidelity model and the accuracy of high-fidelity model, so it is easier to be implemented numerically in real-world aquifers 85 compared with previous methods especially for those aquifers with low Pécelt number and large spatial scale.

Theoretic framework of multi-fidelity model
The multi-fidelity model is designed to calculate TTD, RTD and SAS function by solving the equation of groundwater age distribution. It combines the groundwater age distribution model and particle tracking model to evaluate the time-variant TTD in a system. The two models which have different levels of fidelities will be introduced first.

Groundwater age distribution model
Due to the mixing process undergoing along groundwater flow lines, the groundwater or solutes sampled from a certain point is a mixture of constituents (e.g., groundwater, solutes, and suspended colloids) of various ages. The groundwater age distribution describes the portion of the constituent of a certain age τ at a point. The governing equation for groundwater age distribution was primarily introduced by Ginn (1999) and is modified with variable saturation in this study as follows: where φ is porosity; S w is water saturation; v is Darcy velocity; D is the diffusion coefficient including molecular diffusion and mechanical dispersion. ρ(x, t, τ ) is the age density function of a constituent at time t and point x. If integrating ρ(x, t, τ ) along age-axis (τ -axis), the mass density c(x, t) of this constituent is obtained, that is, The definition of mass density is the ratio of mass to volume. If the concerned constituent is water, mass density is the water density and equal to 1 g/cm 3 ; if the concerned constituent is solute, mass density is the solute concentration. Equation (1) describes the evolution of the age of groundwater in the system in Eularian method. The groundwater age distribution can be obtained by normalizing the age density function ρ(x, t, τ ) over the mass density c(x, t), that is, ρ(x,t,τ ) c(x,t) . The transport of mass density obeys the rules described in the following equation (3), i.e. the advection-dispersion equation (ADE) for conservative 105 solute transport: If groundwater flow is under steady state condition, the groundwater age distribution is a time-invariant distribution, i.e., ρ(x, t, τ ) = ρ(x, 0, τ ) and c(x, t) = c(x, 0), and the governing equation for groundwater age distribution (equation (1)) becomes: Equation (4) is similar to the ADE in equation (3) and can be solved using the methods to solve the ADE.
However, if groundwater flow is under unsteady state condition, the groundwater age distribution is a time-variant distribution, which means that the age density function ρ(x, t, τ ) and groundwater age distribution ρ(x,t,τ ) c(x,t) varies in both clock-time t and field point x. As a result, the groundwater age distribution can only be obtained by solving equation (1) which is a 5D 115 equation with 3 space dimensions, 1 clock-time dimension and 1 age dimension. Equation (1) resolves the evolution of age distribution of a constituent in the system, and it is challenging to solve equation (1) either analytically or numerically.
Referring to the study of Gomez and Wilson (2013), one approach to obtain the numerical solution (i.e. ρ(x, t, τ )) of equation (1) is to decompose equation (1) into a series of independent ADEs. Each ADE has a different initial condition. In this approach, the input time η is introduced as a new variable to replace the age variable τ , such that: Then equation (1) becomes: where ρ is age density function which represents the mass density of water entering the system at the time η as is shown in equation (6). Equation (8) is similar to the ADE, i.e., equation (3), but they have different physical meanings, because the variable ρ depends on the input time η. Equation (8) must be solved for each input event occurred during the study period.
The BTC at the discharge boundary for a certain input time η can be derived from ρ : where τ is the travel time defined as the age at the discharge boundary Σ discharge . Therefore, equation (9) quantifies the mass flux along the discharge boundary for a certain travel time. Since TTD can be decomposed into a series of overprinted BTCs induced by different pulse injections (Harman and Kim, 2014), the time-variant TTDs at the outflow boundary can be constructed by assembling these sets of BTCs.
The advantage of this approach is that it reduces the dimension of the problem by changing the age variable τ into the input 135 time variable η, because the water from the same input event must have the same age. However, in real problems, the model should be simulated for each input event to acquire a good resolution of age distribution. Therefore, in this study, the ADE for a certain input event is defined as high-fidelity model. This high-fidelity model is required to run only a few times for several input events in multi-fidelity model to capture the trend of the variation of age distributions and ensure the multi-fidelity model unbiased.

Particle tracking model
Particle tracking is a popular approach to model the time-variant TTDs based on Lagrangian advective transport. As all particles are tagged with their own ages, the time-variant TTDs can be evaluated by counting the particles at the discharge boundary. To simulate the molecular diffusion and mechanical dispersion of solutes, it is necessary to add an random walk velocity to the advection velocity for each particle. Yang et al. (2018) used the random walk particle tracking model to evaluate time-variant 145 TTD in an agricultural catchment. However, the random walk particle tracking model may suffer from random noise when the number of particles released is insufficient to simulate the diffusion and mechanical dispersion behaviors. This random noise is even inevitable if using particle tracking method in the problems with large space scales. Therefore, in this multi-fidelity model, random walk velocity is excluded from particle velocity, and the particle tracking without considering random walk is defined as the low-fidelity model. Particles are released at the recharge boundary during each rainfall event. Each particle 150 represents a certain volume of water parcels based on the rainfall infiltration rate. As the random walk velocity is not included in the low-fidelity model, the random noise disappears, and the number of particles decreases drastically. Generally, as the scale of the model becomes larger, more particles are required to be released in random walk particle tracking model, so the number of particles in the low fidelity model decreases even more.

Analytical theorem for one dimensional scenario 155
Prior to establishing the multi-fidelity model, the correlation between particle tracking model and ADE is investigated analytically in 1D condition. In the particle tracking model, a particle represents a pulse signal of water or solute entering the domain at a certain injection time. The transport of the pulse signal in 1D domain is studied to illustrate the relation between the two models. Assume the advection velocity is a constant in the entire domain, the transport of the pulse injection can be expressed as: with the initial condition: where ρ is defined in equation (6); δ(x − x 0 ) is the Dirac delta function which represents a pulse signal ocurred at the point x 0 ; M is the intensity of the pulse signal, namely the input mass density. The solution of equation (10) is explicative: when considering molecular diffusion and mechanical dispersion in ADE, equation (10) is manifested as: The analytical solution of equation (13) with respect to the initial condition of equation (11) at infinite domain is: 170 Comparing solution (14) with solution (12), it can be found that as the parameter D goes to zero, the bell-shape curve of solution (14) converges to the Dirac delta function of solution (12), that is, This is consistent with the relation between the two governing equations: equation (13) has an additional diffusion term compared with equation (10). Therefore, no matter under what boundary conditions, the solution of equation (13) will always 175 converge to equation (12) when D approaches zero.
Theoretically, a mapping is capable of being setting up between the solutions (12) and (14) for a given hydrogeological model if the diffusion parameter D is known. This is the theoretical foundation of the multi-fidelity model.

Multi-fidelity model for arbitrary 2D/3D models
In the approach proposed by Gomez and Wilson (2013), the equation of groundwater age distribution is decomposed into a 180 series of ADE to track the transport of water from each rainfall event separately, because the water from the same rainfall event always has the same age. The ADE can be solved either by numerical simulation or by random walk particle tracking model.
For random walk particle tracking model, the problem of random noise has already been discussed in the previous section.
For the numerical simulation method, it should be run separately for each rainfall input, or a different chemical species need to be assigned for each rainfall event to differentiate the transport among different rainfall events. Specifically, if there are Based on the nonlinear correlation between equations (13) and (10) described previously, the solution of equation (8) can be obtained by mapping from the low-fidelity model to the high-fidelity model. The mapping is expressed as a nonlinear autoregressive model: where f H (t, η) is the solution of ADE model which has high-fidelity, f L (t, η) is the solution of particle tracking model without conidering random walk which has low-fidelity, and g represents the nonlinear mapping. Given a certain hydrogeological model, the nonlinear mapping g only depends on the diffusion coefficient D, but the mathematical form of the nonlinear mapping can vary with different hydrogeological models. The mapping of equation (17) can also be applied to link the cumulative 205 breakthrough curves (cBTC) between two models of different levels of fidelity, because the evaluation of BTC through equation (9) and the transformation of BTC to cBTC are both linear operations with respect to age mass density ρ . Therefore, in this study, f H (t, η) is the cBTC from ADE model, and f L (t, η) is the cBTC from particle tracking model without random walk. To ensure the mapping accommodate any hydrogeological model, non-parametric regression with Gaussian processes (GPs) Zhang et al., 2018) is adopted to build up the multi-fidelity model. Let m = [f L (t, η), t, η], then the 210 nonlinear mapping g (f L (t, η), t, η) in equation (17) can be written as g(m) and is assigned to a GP prior: where GP represents the Gaussian process; 0 is the prior mean value; k (m, m ; θ) is the kernel function with hyper-parameters θ. The kernel function in GP prior quantifies the pairwise correlation between f Hi and f Hj and leads to a symmetric covariance matrix K ij = K ji = k (m i , m j ; θ) between the points (m i , f Hi ) and (m j , f Hj ). If m i and m j are close to each other, the 215 correlation between f Hi and f Hj will be high. The term on the right hand side of equation (18) is the GP prior which means the probability distribution of f H is a Gaussian distribution with mean value 0 and variance k (m, m ; θ).
With reference to the work of Perdikaris et al. (2017), the kernel function is expressed as: In equation (19), the kernel function k (m, m ; θ) consists of three components: k 1 denotes the correlation dependence of the 220 modeling time and input time, k 2 denotes the correlation dependence of the results of low-fidelity model, and k 3 denotes the correlation of the regression errors. Each correlation is quantified by a squared exponential covariance function which means that the covariance between f Hi and f Hj decreases exponentially as the distance between m i and m j decreases: where θ p = σ 2 p , (w p,q ) d q=1 , and θ is the vector of hyper-parameters for the covariance function k p (p = 1, 2, 3) in equation

225
(19), σ 2 p is the variance parameter, (w p,q ) d q=1 are the automatic relevance determination weights, and x is a d-dimensional input vector which can be (t, η) or f L in equation (19).
The vector of hyper-parameters θ p (p = 1, 2, 3) in equation (19) determines the performance of the mapping. With a given regression data set (M in , F H ) = (m i , f Hi ) i=1...n obtained from particle tracking model and ADE model, the hyper-parameters Θ = [θ 1 , θ 2 , θ 3 ] can be estimated by minimizing the negative marginal log-likelihood: where where K x−y represents the covariance between x and y. Then, the posterior distribution is exactly the marginal distribution of the joint Gaussian distribution: where the predicted mean is µ = K T Hnew−H K −1 H−H F H , and the predicted variance is σ 2 out = K Hnew−Hnew −K T Hnew−H K −1 H−H K Hnew−H . In this way, the BTCs for all rainfall events can be inferred from the mapping based on the solution of the particle tracking model without random walk. For a complex hydrogeological system, it is difficult to derive the mapping analytically, so the non-parametric regression model is employed to estimate the mapping. Overall, the number of rainfall events considered in ADE model is reduced compared with the groundwater age distribution model (Gomez and Wilson, 2013); and the number of 245 particles released is reduced significantly compared with the random particle tracking model (Yang et al., 2018). The multifidelity model also makes the evaluation of time-variant TTD simpler to be implemented compared with the approach that evaluates time-variant TTD numerically in Laplace space (Cornaton, 2012). There is a tradeoff between the computational cost and accuracy. Larger the size of regression data set is, the more accurate the regression will be, but the computational cost will also be higher. The problem of how to determine the size of regression data set will be discussed in the illustrative example.

Evaluation of the time variant TTD and SAS function
Referring to the definition of forward TTD of water (fTTD), each BTC can be transformed to a fTTD by normalizing it with the total rainfall recharge volume. Let − → p T (η, τ ) represents the fTTD conditional to rainfall event occurred at time η, then where τ is the age known as travel time at the outflow site, P (η) is the rainfall recharge volume at time η, BT C(τ, η) is the 255 BTC of water which is input at time η. Equation (24) displays a one-to-one mapping between the BTC and fTTD. The fTTD quantifies the portion of water leaving the system at time t within the same rainfall event.
To study the time-variant behavior of the storage selection process, backward TTD of water (bTTD) and residence time distribution of water (RTD) are evaluated to explore the dynamic storage selection process. Similar to the evaluation of fTTD, bTTD can be calculated by normalizing the BTCs with the discharge rate as follows, where ← − p T (t, τ ) is the bTTD. Given a certain time t in equation (25), the evaluation of bTTD requires the BTCs before the time t. This is because the water in the hillslope discharging at the time t comes from the rainfall occurred before the time t.
Up to now, most studies employ BTC to calculate fTTD and bTTD directly (Harman and Kim, 2014;Kim et al., 2016), but few studies calculate RTD by the BTC. As the results obtained from multi-fidelity model are a set of BTCs conditional on each 265 rainfall event, it is necessary to put forward a method to calculate RTD using BTCs. The cBTC is defined as: The cBTC cBT C(t, η) characterizes the total volume of water or solutes from the same rainfall event at the discharge boundary before the time t. Therefore, the water or solutes that still remains in the system at the time t is equal to P (η) − cBT C (t, η).
where P (η) is the total volume of rainfall which is infiltrated at the time η. By assembling these cBTCs, the denominator in equation (27) represents the total storage volume at time t.
Having calculated the bTTD and RTD, the storage selection (SAS) function can be quantified according to its definition.
To characterize the time-variant behaviors of storage selection process conveniently, the fractional SAS function (fSAS) is 275 calculated based on the results of bTTD and RTD (van der Velde et al., 2012). Because there is a one-to-one mapping between the groundwater residence time (groundwater age) and cumulative RTD, fSAS function is defined as a function of cumulative RTD and time: where ω Q (t, P s ) is the fSAS and P s (t, τ ) is the cumulative RTD. Many studies employ beta distribution to quantify this fSAS a lower solute concentration and a higher water-rock reaction rate compared with the groundwater at the down-slope (Braun et al., 2016;Rempe and Dietrich, 2014). This leads to the wedge-shaped profile of weathered rock zone.
The numerical simulation of ADE model and particle tracking model without random walk are implemented in the software COMSOL Multiphysics (COMSOL AB, 2018). Three modules in COMSOL are employed, namely Richards' equation module, transport of dilute species in porous media module, and mathematical particle tracing module.

295
In Richards' equation module, hydraulic head and Darcy velocity are solved to couple with the other two modules. The left boundary and bottom boundary are no-flow boundaries, the right boundary is a fix-head boundary with hydraulic head equal to 4.5 meters. The top boundary is a time-variant rainfall infiltration boundary. Figure 2 shows the variation of 2-year nonstationary rainfall volume per day with 2728.5 mm rainfall volume in total during the first year and 2678.5 mm during the rate is set to be 0.2. The modeling period is 10046 days. The 2-year rainfall will be repeated every two years. As the goal of this study is to explore the time-variant TTD and storage selection process in hillslope, the spatial heterogeneity of hydrogeological properties is ignored and only one layer with constant hydraulic conductivity is considered in this hypothetic model. The detailed hydraulic parameters are shown in Table 1. Theoretically, the initial condition for Richards' equation can be arbitrary, but in order to reach a periodic quasi-steady state rapidly, the initial condition is set to be: where the initial hydraulic gradient is equal to the slope gradient. After running a spin-up period of 4 years, the solution of Richards' equation reaches a periodic quasi-steady state. Then the solution of this module is coupled with the other two modules through Darcy velocity to further solve for the BTCs.
The transport of dilute species in the porous media module in COMSOL is used to solve for the groundwater age distribution 310 model. This is because equation (8) has the same mathematical form as the solute transport equation. Therefore, the mass density in equation (8) can be regarded as a solute concentration while modeling. As is discussed in the section 2.1, this module should be run for each single rainfall event independently to obtain BTCs. The advection velocity required in this  The mathematical particle tracing module in COMSOL is employed to solve for the BTCs without considering molecular diffusion and mechanical dispersion. The particles are released along the hillslope surface uniformly during rainfall. The number of particles released at the upper boundary is the same for each rainfall event. Each particle represents a certain volume of water parcels depending on the rainfall intensity. As the diffusion effect is not considered in this module, the movements 325 of these particles are deterministic based on the groundwater velocity. A particle counter is set at the discharge boundary to keep counting the arriving particles to evaluate the BTCs. The BTCs calculated in this module are taken as the solution of low-fidelity model, because the solution is biased against the high-fidelity model due to the absence of molecular diffusion and mechanical dispersion effect. shows the spatial distribution of age density function ρ(x, t, τ ) (unit: kg/m 3 ) in high-fidelity model. These two results displayed in figure 3 are for the water entering the system during the rainfall event occurred at day 165. In figure 3, the age density function spreads around the red dots due to molecular diffusion and mechanical dispersion, and this is consistent with the analytical 1D correlation of equation (15).

Selecting regression data set for multi-fidelity model
Having set up these modules in COMSOL, there are two more factors affecting the result of the multi-fidelity model and should be specified before running the multi-fidelity model, that is, (1) the selection of specific high-fidelity models for regression data set and (2) the number of selected high-fidelity models in regression data set. To study how these two factors affect the results, correlation, error and convergence analysis are implemented. The nonlinear regression with Gaussian process used in multi-340 fidelity model is realized in python using the Gpy package (GPy, since 2012).
There are many methods available to assist in selecting the high-fidelity models for regression data set. In this study, the distance method is employed and is fully described by Scheidt and Caers (2009);Josset et al. (2015). Specifically, as the distance of selected parameters between high-fidelity models becomes short, the results of multi-fidelity model are closer to the results of high-fidelity models. In high-fidelity models, the selected parameter is the time when rainfall occurred. If the 345 two rainfall events are very close, the two BTCs at the discharge boundary will be similar due to the similarity of hydraulic conditions. Therefore, the selection of the data set for regression is based on the 'distance' between the rainfall time, i.e. the interval between two rainfall events. In this hypothetical model, there are 290 rainfall events every 2 years, and the rainfall events are selected successively from 290 events with a constant distance. The cBTCs calculated from high-fidelity models within the selected rainfall events make up the regression data set for regression.

350
The distance among the rainfall events is used to control the number of high-fidelity models in regression data set. Therefore, it is vital to determine a reasonable distance, and this will be discussed in the following error analysis and convergence analysis.
To evaluate the error of regression and study the convergence of the multi-fidelity model, the high-fidelity model is solved for all the 290 rainfall events intensively to generate the test data set. Practically, only error analysis is required during the implementation of multi-fidelity model, so the size of test data set can be much smaller .

355
The 2-norm error of a certain BTC reflects the overall deviations of the predicted BTC from the BTC in the test data set, and is calculated below: where e 2 is the 2-norm error of a certain BTC, N is the number of points in the BTC, y n and y(t n ) are the points in BTC calculated by multi-fidelity model and high-fidelity model respectively. The 2-norm error is used to quantify the goodness-of-fit 360 of the calculated BTC by multi-fidelity model.
In the convergence analysis, the size of the regression data set increases by shortening the interval between rainfall events.   Figure 4, the goodness-of-fit improves as the interval among rainfall events shortens. This confirms that multi-365 fidelity model converges as the size of regression data set increases.

As shown in
Then an interval of 10 rainfall events is employed to calculate the predicted errors of cBTCs in details. The predicted error is defined as the difference value between the predicted BTC and the BTC evaluated by high-fidelity model. Figure 5 compares the predicted errors for all 290 realizations between the multi-fidelity model and low-fidelity model. Each grey curve in Figure   5 represents the predicted errors of a BTC. The mean predicted errors in Figure 5a are very close to zero, which indicates that 370 the predicted curves are unbiased and the selected high-fidelity models with the interval of 10 days for regression data set is representative of all high-fidelity BTCs. However, the mean predicted errors in Figure 5b are biased against the negative value.
The biased errors are identified as the system errors. It implies that the solute concentration in discharge predicted by particle tracking without random walk is generally smaller than that predicted by high-fidelity model, because except for the advection process, the molecular diffusion and mechanical dispersion can also transport solute from recharge boundary to discharge 375 boundary.
The goodness-of-fit for model with 'distance' of the interval 10 is shown in Figure 6 directly. It can be seen that most of the predicted cBTCs match the high-fidelity model well. Also, these inconsistencies between predicted cBTCs and cBTCs 10 rainfall events is recommended for the application of the multi-fidelity model in other problems, because the interval of 10 rainfall events is equivalent to about one-month interval in wet season and about two-or-three month interval in dry season.

Errors of time-variant bTTD, RTD and fSAS function
With reference to equations (25) and (27) In figures 7 and 8, the bTTD and cumulative TTD are significantly distinctive among the multi-fidelity model, high-fidelity model and low-fidelity model, but the difference of the RTD and cumulative RTD among the three models is not identifiable.
This phenomenon is ascribed to the difference of kinetics (i.e., molecular diffusion and mechanical dispersion) of water parcels   Figure 6. The comparison of predicted cBTCs (multi-fidelity model) with the solute cBTC (high-fidelity model) and the particle cBTCs (low-fidelity model) on different rainy days.