Model description paper 01 Jul 2020
Model description paper  01 Jul 2020
MFIT 1.0.0: MultiFlow Inversion of Tracer breakthrough curves in fractured and karst aquifers
 Université de Poitiers, CNRS, UMR 7285 IC2MP, 40 Avenue du Recteur Pineau, 86022 Poitiers, CEDEX, France
 Université de Poitiers, CNRS, UMR 7285 IC2MP, 40 Avenue du Recteur Pineau, 86022 Poitiers, CEDEX, France
Correspondence: Jacques Bodin (jacques.bodin@univpoitiers.fr)
Hide author detailsCorrespondence: Jacques Bodin (jacques.bodin@univpoitiers.fr)
More than half of the Earth's population depends largely or entirely on fractured or karst aquifers for their drinking water supply. Both the characterization and modeling of these groundwater reservoirs are therefore of worldwide concern. Artificial tracer testing is a widely used method for the characterization of solute (including contaminant) transport in groundwater. Tracer experiments consist of a twostep procedure: (1) introducing a conservative tracerlabeled solution into an aquifer, usually through a sinkhole or a well, and (2) measuring the concentration breakthrough curve (BTC) response(s) at one or several downstream monitoring locations, usually spring(s) or pumping well(s). However, the modeling and interpretation of tracer test responses can be a challenging task in some cases, notably when the BTCs exhibit multiple local peaks and/or extensive backward tailing. MFIT (MultiFlow Inversion of Tracer breakthrough curves) is a new opensource Windowsbased computer package for the analytical modeling of tracer BTCs. This software integrates four transport models that are all capable of simulating single or multiplepeak and/or heavytailed BTCs. The four transport models are encapsulated in a general multiflow modeling framework, which assumes that the spatial heterogeneity of an aquifer can be approximated by a combination of independent onedimensional channels. Two of the MFIT transport models are believed to be new, as they combine the multiflow approach and the doubleporosity concept, which is applied at the scale of the individual channels. Another salient feature of MFIT is its compatibility and interface with the advanced optimization tools of the PEST suite of programs. Hence, MFIT is the first BTC fitting tool that allows for regularized inversion and nonlinear analysis of the postcalibration uncertainty of model parameters.
 Article
(12467 KB) 
Supplement
(100 KB)  BibTeX
 EndNote
Artificial tracer testing is one of the most valuable methods for the characterization of flow and solute transport in groundwater. Tracer experiments consist of a twostep procedure: (1) introducing a known mass of a tracer species into an aquifer, usually through a sinkhole or well, and (2) measuring the concentration breakthrough curve (BTC) response(s) at one or several downstream monitoring locations, usually spring(s) or pumping well(s). The analysis of a tracer BTC is best done by fitting a modelcomputed time–concentration curve to the measured values. Although spatially distributed numerical models (e.g., MODFLOW/MT3DMS or FEFLOW) can be used for this purpose, simpler (i.e., spatially lumped) models are generally used, at least in the early stages of tracer studies, either because of time constraints or because of a lack of model input data. A number of computer codes for BTC fitting have been developed in recent decades: CATTI (Sauty et al., 1992), TRACI (Käss, 1998, 2004), OTIS (Runkel, 1998), STANMOD (van Genuchten et al., 2012), TRAC (Gutierrez et al., 2013), OMMADE (Tinet et al., 2019), and OptSFDM (Gharasoo et al., 2019). Note that STANMOD integrates a number of former codes, including the widely used CXTFIT code developed by Parker and van Genuchten (1984) and Toride et al. (1999). Despite the range of possibilities offered by these programs, the fitting and interpretation of tracer BTCs remains a challenging task in some cases, notably for BTCs exhibiting multiple local peaks and extensive backward tailing. Such BTC shapes, which fall in the general category of nonFickian (or anomalous) transport (Berkowitz et al., 2006; Neuman and Tartakovsky, 2009), are frequently observed in fractured and karst aquifers (Tsang and Neretnieks, 1998; Streetly et al., 2002; Massei et al., 2006; Loefgren et al., 2007; Goldscheider et al., 2008; Field and Leij, 2012; Bertrand et al., 2015; Yang et al., 2019).
To the best of the author's knowledge, only the TRACI and OMMADE programs are able to simulate multimodal BTCs. Unfortunately, these two programs suffer from some limitations both in terms of ease of use and with respect to their modeling and calibration capabilities. For instance, the TRACI software has not been maintained since 2004 and can only be used on physical or virtual computers running Windows operating system versions from Windows 98 to Windows 7. Another drawback of TRACI is the inability of the inversion (automated calibration) algorithm included in the software to handle multimodal BTCs. Each local concentration peak must be sequentially fitted through a manual (trialanderror) calibration procedure. The OMMADE program was written as a Python script and has neither a graphical user interface (GUI) nor inverse modeling functionality. The purpose of this paper is to present a new opensource GUIbased software, named MFIT, that aims to help in the interpretation of single or multiplepeak and/or heavytailed BTCs. MFIT stands for “MultiFlow Inversion of Tracer breakthrough curves”. The MFIT software integrates four transport models that can be tested against field and laboratory tracer BTCs with the assistance of the PEST automated calibration and uncertainty analysis routines (Doherty, 2019a). In its current version, the scope of the software is limited to tracer tests involving nonreactive tracer species and performed in steady flow conditions. These assumptions are maintained throughout the paper.
The remainder of this paper is organized as follows. Section 2 discusses the possible origins of multiple peaks and long tails in tracer BTCs and presents the conceptual and mathematical framework of the transport models integrated in MFIT. The code implementation and coupling with PEST for automated BTC fitting are discussed in Sect. 3. In Sect. 4, the accuracy of MFITcomputed BTCs is verified against CXTFIT and TRACI simulations for five test cases. An additional test is presented in Sect. 5 to assess the reliability of a new multistart method that was specifically developed to improve the automatic optimization of the model parameters. In Sect. 6, we illustrate the use of the software by analyzing tracer BTCs obtained in the karst aquifer of the Hydrogeological Experimental Site (HES) in Poitiers, France. The summary and conclusions are presented in Sect. 7. A number of acronyms, model abbreviations, and model parameters are employed throughout this paper. Two glossaries, Tables A1 and A2, are provided in Appendix A for easy reference.
Under the already mentioned assumptions (nonreactive tracer, steadystate flow), multimodal BTCs unequivocally indicate that a number of tracerplume splittings occurred somewhere between the injection site and the monitoring point. Although injection artifacts may be involved in some cases (see Guvanasen and Guvanasen, 1987), tracer splitting most commonly originates from the spreading (transverse dispersion) of the solute into areas of contrasting flow velocities; see Moreno and Tsang (1991), SiirilaWoodburn et al. (2015), and Boon et al. (2017). More precisely, assuming a singlepulse tracer injection signal, multimodal BTCs reflect a threestep process: (1) tracer spreading into different flowing or nonflowing aquifer subdomains characterized by different transit or residence times, (2) tracer motion within each subdomain with little or no exchange between the different subdomains, and (3) convergence (mixing) of the subtracer fluxes somewhere upstream from, or at, the monitoring point. The different models that have been proposed in the literature for simulating multimodal tracer BTCs share a common “multiflow” approach initially proposed by Zuber (1974) for the modeling of layered aquifers. In this approach, which is depicted in Fig. 1, the flow system is described as a juxtaposition of a number of onedimensional (1D) channels that are connected by a single common diverging (splitting) node at the entrance to the system and a single common converging (mixing) node at the outlet.
In the MultiDispersion Model (MDM) proposed by Maloszewski et al. (1992) and implemented in the TRACI software, the transport along each channel is assumed to obey the onedimensional (1D) advection–dispersion equation (ADE), and no mass exchange is allowed between different channels. In the dualadvection dispersion equation (DADE) proposed by Field and Leij (2012), only two channels are considered. The tracer is transported by advection and dispersion along each channel, and mass exchanges between the two domains are possible. These exchanges are assumed to be governed by a firstorder process. The transport model implemented in the OMMADE code can be viewed as a generalization of the DADE model, where (i) a larger number of channels can be used, (ii) each channel can be discretized to a number of subelements with different hydraulic and transport properties, and (iii) some channels can be specified as nonflowing (stagnant) water volumes. Mass exchanges between the different channels (either flowing or nonflowing) are likewise modeled as a firstorder process. As pointed out above, the production of a multimodal BTC requires little or no exchange between the subtransport domains; otherwise, the mixing of the mass fluxes would rehomogenize the subtracer plumes. In accordance with this principle, small exchange coefficient values must be used in the DADE and OMMADE models for simulating multimodal BTCs, and this approach makes these models converge toward the MDM.
The interpretation of the longtail behavior of a BTC may be more difficult than that of multiple peaks, as different mechanisms can be involved. The possible sources of extensive BTC tailing can be listed as follows: (i) tracer retention that produces a decaying boundary condition at the injection site; (ii) tracer splitting into wellseparated flow paths and then downstream convergence, mixing, and overlapping of the individual pathway responses; and (iii) mass exchanges between flow domains characterized by different transit or residence times. The abovelisted processes are referred to below as “injection decay”, “multiflow overlapping”, and “multiflow exchanges”, respectively. The MDM can simulate longtailed BTCs as a result of multiflow overlapping. Multiflow exchanges are the core of the DADE model, and both multiflow overlapping and multiflow exchanges can be combined in the OMMADE model. A number of other models have been proposed in the literature for simulating unimodal longtailed BTCs; see reviews in Bodin et al. (2003b), Neuman and Tartakovsky (2009), Zhang et al. (2009), and Dentz et al. (2011) and examples of recent works in Field and Leij (2014) and Labat and Mangin (2015). The two most commonly used models for the analysis of artificial tracer tests are the tworegion nonequilibrium (2RNE) model by Toride et al. (1993), implemented in the CXTFIT code, and the SingleFracture Dispersion Model (SFDM) by Maloszewski and Zuber (1990), implemented in TRACI and OptSFDM software. Both the 2RNE model and SFDM assume mass exchange between a single mobile (flowing) domain and a single immobile domain. A key distinction between the 2RNE model and SFDM is the formulation of mass exchange, which is described as a firstorder process in the 2RNE model (as in the DADE and OMMADE models) and as a secondorder (diffusion) process in the SFDM.
As already noted, multimodal and longtailed BTCs are typical of tracer tests performed in fractured and karst aquifers. A common feature of both aquifer types is the existence of lowhydraulicresistance pathways provided by the fractures and karst conduits (Tsang and Neretnieks, 1998; Worthington and Ford, 2009). A generic multiflow modeling approach is therefore intuitively appealing for the interpretation of tracer tests in fractured and karst aquifers. Of course, the actual (and generally unknown) geometry of the discrete flow network experienced by the tracer is likely more complex than that depicted in Fig. 1. The channels are therefore not assumed to represent individual fractures or karst conduits but are lumped submodels of the main flow routes used by the tracer through the fractures or karst conduit network. The four transport models integrated in the MFIT software are based on the multiflow approach. The first model is a reimplementation of the MDM. The second model is a variant of the MDM that assumes an exponentially decaying injection of the tracer concentration at the inlet of the flow system. In the third and fourth models, the doubleporosity concept (2RNE model and SFDM) is applied at the scale of the individual channels. It is unclear whether this idea of combining multiflow and doubleporosity systems is new. In the TRACI software, it is technically possible to fit a series of SFDM curves to a multimodal tracer BTC and then calculate the mean combined model curve, but to the best of the author's knowledge, this method has never been discussed or applied in the literature. A possible reason is the increasing number of fitting parameters, which makes the inverse problem more complicated. Among the challenges related to the inversion of a multiflow model is the inherent problem of nonuniqueness (or equifinality). A variety of parameter sets can yield nearly identical simulated BTCs, because the change in the value of a parameter of a given channel can be compensated by modifying at least one other parameter that pertains to this same channel or the parameters of the other channels. This nonuniqueness causes the inverse problem to be illposed in the sense of Hadamard (1902) and requires the use of advanced optimization methods, such as regularization, to make the inverse problem tractable (Tikhonov and Arsenin, 1977; Moore and Doherty, 2006; Zhou et al., 2014).
In this article, the combination of multiflow and doubleporosity systems is referred to as the multidouble porosity (MDP). The immobile domain that is assigned to each flow channel is assumed to describe the porous rock matrix in contact with the fractures or karst conduits and/or any other stagnant water zones (e.g., pool volumes) adjacent to the main tracer pathways. For each of the four MFIT models, the channels are assumed to be independent of each other, i.e., no mass exchange is allowed between the channels. Actually, this assumption is mathematically convenient rather than physically motivated. As already indicated, the channels are abstractions of the real main tracer pathways, which may cross (and therefore exchange between) each other between the injection site and the monitoring point. Assuming fully separated channels allows for analytical modeling of mass fluxes in the multiflow system, and this approach makes the inversion of model parameters computationally more efficient (see discussion in Sect. 3).
The governing equations of the transport models are given as follows. The concentration at the outlet of a multiflow system as depicted in Fig. 1 can be calculated from the mass flux balance as follows:
where Q [L^{3} T^{−1}] is the total system flow rate; C [M L^{−3}] is the outflow concentration; N is the number of flow channels; the subscript j denotes the flow channel index; and Q_{j} [L^{3} T^{−1}] and C_{j} [M L^{−3}] are the flow rate and concentration in the jth channel, respectively.
The mathematical equations that have been used by Maloszewski et al. (1992) in the MDM to describe the solute transport in each flow channel are the 1D ADE as follows:
and its analytical solution for the case of an instantaneous solute injection in a semiinfinite medium with both injection and detection in flux, which is expressed as (Kreft and Zuber, 1978)
where t [T] is the time variable; x_{j} [L] is the spatial coordinate along the jth flow channel; u_{j} [L T^{−1}] and D_{j} [L^{2} T^{−1}] are the advection velocity and the dispersion coefficient, respectively; m_{j} [M] is the part of the solute mass that flows through the jth channel; and T_{0j} [T] and Pe_{j} [–] are the mean transit time and Péclet number, respectively, which are expressed as
where L_{j} [L] is the length of the jth pathway. Substituting Eq. (3) into Eq. (1) yields
The calibration of Eq. (6) against a tracer test BTC requires determination of the total system flow rate Q; the number N of flow channels; and for each flow channel, the values of m_{j}, T_{0j}, and Pe_{j}. In this work, we generalize the abovedescribed method by considering alternative models for the transport in individual channels and substituting the related analytical expressions of C_{j} into Eq. (1). The analytical transport models that are considered are (i) the solution of Eq. (2) for the case of a decaying injection boundary condition, (ii) the SFDM, and (iii) the 2RNE model.
The analytical solution of Eq. (2) for the case of a decaying injection boundary condition ${C}_{j}({x}_{j}=\mathrm{0},t)={C}_{\mathrm{0}}\mathrm{exp}({\mathit{\lambda}}_{j}t)$ was derived by Marino (1974) and can be written in the following form:
where C_{0} [M L^{−3}] is the initial (maximum) injection concentration at the inflow boundary, λ_{j} [T^{−1}] is the time decay constant, and
The SFDM developed by Maloszewski and Zuber (1990) describes solute transport in a doubleporosity fracturematrix system. The considered transport mechanisms are advection–dispersion in the fracture and diffusion in the surrounding rock matrix. The fracture is idealized as a parallelplate channel, and the matrix diffusion is assumed to be unlimited, i.e., not influenced by the fluxes from other fractures. The transport equations can be written as follows:
where C_{j} [M L^{−3}] and C_{pj} [M L^{−3}] are the solute concentrations in the flow channel and in the rock matrix, respectively; θ_{pj} [–] is the matrix porosity; D_{pj} [L^{2} T^{−1}] is the molecular diffusion coefficient in the matrix; b_{j} [L] is the halfaperture of the flow channel; and y_{j} [L] is the spatial coordinate perpendicular to the channel extension. The solution to Eqs. (9) and (10) for the case of an instantaneous injection is
where ξ [T] is the integration variable and β_{j} [T${}^{\mathrm{1}/\mathrm{2}}$] is the socalled diffusion parameter defined as
Coats and Smith (1964) proposed a different mathematical formulation of solute mass exchange between flowing and stagnant water regions in doubleporosity media, which is well known in the literature either as the mobile–immobile (MIM) model or the 2RNE model as
where θ_{j} [–] and θ_{imj} [–] are the mobile and immobile volumetric water contents, respectively; C_{imj} [M L^{−3}] is the concentration in the immobile domain; and α_{j} [T^{−1}] is a firstorder mass transfer coefficient. The two main differences with respect to the SFDM are (i) the dualdomain formulation of the problem (mobile and immobile regions are assumed to coexist at each point in space, and this assumption differs from the parallelplate channel geometry in the SFDM) and that (ii) the solute mass exchange between mobile and immobile domains is assumed to be governed by a firstorder process, whereas the SFDM refers to the secondorder diffusion Eq. (10). Building on a general set of analytical solutions developed by Toride et al. (1993), the solution of the 2RNE model for the case of an instantaneous injection can be written as follows:
where I_{1} is the modified Bessel function of the first kind, τ [L] is the integration variable,
It is notable that when ψ_{j}=1 and ω_{j}=0, Eq. (15) simplifies to Eq. (3). Table 1 summarizes the parameters of the four MFIT transport models.
The four analytical models described in the previous section have been implemented in C$++$ and compiled as executable Windows programs named MDMi.exe (for MDM, instantaneous injection), MDMed.exe (for MDM, exponentially decaying injection), MDP_SFDM.exe, and MDP_2RNE.exe. The code and executable files are freely available on the public Zenodo repository: https://doi.org/10.5281/zenodo.3470751. In the MDP_SFDM and MDP_2RNE programs, the numerical evaluation of the integrals in Eqs. (11) and (15) is performed using the QAG adaptive integration routine from the GNU Scientific Library with a 61point Gauss–Kronrod rule and a relative error convergence criterion of 10^{−2}. These four programs can be run as console applications to solve a direct (forward) problem, i.e., computing a series of time–concentration values for a given set of model parameters. Both the input and output files are in ASCII format and can be edited with any text editor program for pre and/or postprocessing. A convenient alternative is to use the MFIT software as a GUI for these applications. The MFIT software has been developed using the C$++$ Builder environment (Embarcadero RAD Studio 10.1 Berlin) and provides a GUI for (i) importation and graphic visualization of userprovided BTC data; (ii) parameterization, direct running, and graphical output of the analytical transport models; (iii) inversion (automatic calibration) of model parameters for optimal curve fitting; and (iv) assessment of the uncertainty of calibrated parameter values.
The optimization and uncertainty analysis of the model parameters for a given number of flow channels are carried out using PEST routines (Doherty, 2019a, b). The influence of the number of channels on the model fitting performance can be analyzed once a series of calibrations has been performed for a variety of channel numbers, as illustrated below. PEST is a publicdomain modelindependent program suite that has been widely used over the past 2 decades, notably in the field of surface and subsurface hydrology (e.g., Long, 2015; Woodward et al., 2016; Gaudard et al., 2017; Wang et al., 2019). The theoretical framework and full range of capabilities of the PEST software are well documented (Doherty et al., 2010; Doherty, 2015, 2019a, b) and are not repeated here. Only the concepts and methods that were deemed to be the most relevant to the multiflow modeling approach and that have been made accessible through the MFIT GUI software are briefly reviewed below.
PEST is based on a gradient optimization method and, as such, requires the derivatives of model outputs with respect to the adjustable model parameters to be calculated in each iteration for implementing the Jacobian (sensitivity) matrix. As pointed out by Doherty (2015), the accuracy of these derivative calculations is critical to the performance of the PEST optimization algorithm. In the MFIT program suite, most of the model partial derivatives are calculated analytically and externally provided to PEST. This approach ensures both the accuracy and speed of this part of the optimization process. Less straightforward partial derivative expressions were derived using MAPLE and exported as C code using the MAPLE code generation routine. The partial derivative functions were implemented in the MDMi, MDMed, MDP_SFDM, and MDP_2RNE programs and are processed during the PEST system calls to these programs by providing an optional “/d” commandline argument to the program name. In a few cases, however, the partial derivatives cannot be calculated analytically, as they involve undefined limits. Such is the case for the derivatives of Eq. (15) with respect to the parameters ψ_{j}, L_{j}, and T_{0j}. In these cases, the partial derivatives are computed by PEST using finite differences.
The calibration of a multiflow transport model against a tracer BTC is hampered by two wellknown issues in inverse modeling: (i) model nonlinearity and (ii) solution nonuniqueness. Both issues may cause numerical instabilities that can prevent the inversion algorithm from converging to the optimal solution. PEST includes two regularization methods that can be used either individually or together to guide the optimization process. The singular value decomposition (SVD) method subtracts parameter combinations for which the tracer BTC is uninformative. The inversion is conducted on the basis of a reduced set of orthogonal linear combinations of the model parameters rather than attempting to estimate the parameters individually. The Tikhonov regularization method provides a different but complementary strategy, where the information content of the tracer BTC is supplemented with expert knowledge pertaining to the model parameters. When using Tikhonov regularization, the objective function that is minimized by PEST is defined as the sum of two terms. The first term is the “measurement objective function” and is defined as the sum of the squared weighted differences between the real tracer BTC and the modelsimulated curve. The second term is referred to as the “regularization objective function” and acts as a penalty function for deviations from some preferred parameter conditions. Two Tikhonov regularization options have been implemented in MFIT. The first option, referred to as “preferred homogeneity”, promotes a solution of minimum variance for the model parameters pertaining to the different channels. In the second option, referred to as “preferred value”, the optimization process seeks the solution that is the closest to some prior estimates of the model parameters.
Unfortunately, neither SVD nor Tikhonov regularization can guarantee that the PEST optimization algorithm will converge to the global optimal solution in the parameter space. Where local minima exist in the objective function, which is the rule rather than the exception with nonlinear models, the optimization process may become trapped and fail to identify existing better solutions (Singh et al., 2012; Espinet and Shoemaker, 2013; Abdelaziz et al., 2019). A central issue in this case is the sensitivity to initial parameter values, i.e., different initial parameter sets may lead to different optimized solutions. Global optimization methods have been proposed in the literature to overcome this issue; see Arsenault et al. (2014) for a review and comparison of various algorithms. The PEST program suite includes two such global optimizers based on the SCEUA method (Duan et al., 1992) and the CMAES method (Hansen and Ostermeier, 2001). The corresponding programs are named SCEUA_P and CMAES_P, respectively. It must be noted, however, that global optimization methods suffer from their own drawbacks, including sensitivity to tuning parameters and low computational efficiency. An alternative strategy to improve the chances of convergence toward the global optimum with gradientbased methods is the “multistart” approach, which consists of repeating the optimization process starting from different initial parameter value sets (Skahill and Doherty, 2006; Piotrowski and Napiorkowski, 2011). Such a strategy has been implemented in the MFIT software. The key principle of the proposed algorithm is that rather than conducting the optimization for a fixed number N of channels only, a series of automatic tracer BTC fittings is performed for a decreasing number of channels ranging from N_{max} to 1. The main steps of the MFIT multistart algorithm are detailed as follows:

The first optimization is performed by considering the maximum number of flow channels N_{max}. The initial transport parameters are automatically tuned by MFIT to obtain N_{max} wellseparated concentration peaks. For this goal, the tracer BTC is first analyzed to determine the times T_{5} and T_{95}, which are defined as
$$\begin{array}{}\text{(18)}& {T}_{\mathrm{5}}=max\left({T}_{\text{5th}},\mathrm{1.1}\times {T}_{min}\right),\end{array}$$$$\begin{array}{}\text{(19)}& {T}_{\mathrm{95}}=min\left({T}_{\text{95th}},\mathrm{0.9}\times {T}_{max}\right),\end{array}$$where T_{5th} and T_{95th} are the earliest and latest times at which the concentration values are above and below 5 % of the maximum concentration value, respectively, and T_{min} and T_{max} are the minimum and maximum time values of the userprovided BTC, respectively. The mean travel times, T_{0j}, are then uniformly distributed between the times T_{5} and T_{95}. Next, the initial Péclet number, Pe_{j}, of each channel is calculated as
$$\begin{array}{}\text{(20)}& {\mathit{Pe}}_{j}={\left(\mathrm{15}{\displaystyle \frac{{N}_{max}{T}_{\mathrm{0}j}}{{T}_{\mathrm{95}}{T}_{\mathrm{5}}}}\right)}^{\mathrm{2}}.\end{array}$$Equation (20) is based on a semiempirical relationship between the standard deviation of travel times for transport by advection and dispersion, ${\mathit{\sigma}}_{j}={T}_{\mathrm{0}j}\sqrt{\mathrm{2}/{\mathit{Pe}}_{j}}$ (see Bodin et al., 2003a; Eq. 10), and the time span of the jth concentration peak, which is on the order of 6σ_{j}. The constraint of wellseparated concentration peaks may be formulated as $\mathrm{6}{\mathit{\sigma}}_{j}{N}_{max}\ll ({T}_{\mathrm{95}}{T}_{\mathrm{5}})$, which is verified by Eq. (20). The initial values of the other transport parameters in Eqs. (7), (11), and (15) are chosen to minimize the tailing effect due to noninstantaneous injection or solute mass exchange between flowing and stagnant water regions as follows: γ_{j}=0.1, β_{j}=0.001, ψ_{j}=0.9, and ω_{j}=0.05.

Once the optimization has been performed for the N_{max} channel model, the next step is to optimize the transport parameters for ${N}_{max}\mathrm{1}$ channels. The multistart optimization approach begins here as not only one but N_{max} optimizations are performed in this step. The initial parameter values for the ${N}_{max}\mathrm{1}$ channels are initialized from the previously optimized N_{max} channel solution by sequentially removing one of the channels. Only the solution corresponding to the lowest sum of the squared weighted differences between the tracer BTC and modelsimulated curve is retained.

This procedure is repeated up to the singlechannel solution. The total number of PEST optimizations is ${N}_{max}({N}_{max}+\mathrm{1})/\mathrm{2}$.
Calling the multistart algorithm has been made optional in MFIT, as this algorithm significantly increases the computational cost and running time of the optimization process. However, experience has shown that the multistart approach can truly improve the model fit results and can be worth the effort in many circumstances. A comparison between optimizations conducted by the PEST multistart algorithm and the global SCEUA and CMAES methods was conducted in this study and is discussed in Sect. 6.
Because of the nonuniqueness of the inverse problem, some uncertainties may be associated with the PESToptimized model parameter values. A nonlinear analysis method has been implemented in MFIT for the assessment of postcalibration parameter uncertainty. The method is essentially similar to that described by Fang et al. (2019) and relies on the use of the PREDUNC7 and RANDPAR utilities documented in the PEST manual (Doherty, 2019b). The algorithm can be described by the following steps: (1) compute a linear approximation to the posterior parameter covariance matrix using PREDUNC7; (2) sample the posterior parameter covariance matrix and generate multiple calibrationconstrained random parameter sets with RANDPAR; (3) recalibrate each parameter set with PEST up to achieving a level of fit fairly similar to the original calibration result (a tolerance of +5 % for the measurement objective function is allowed by MFIT); and (4) compute histograms of the recalibrated parameter values. The following two assumptions underlie this method: (i) the upper and lower parameter bounds specified by the user for the PEST inversion reflect the prior (expert knowledge) parameter uncertainty, and (ii) the model parameters are statistically independent from a prior point of view. This second assumption is relaxed through the recalibration process.
The robustness of the PEST inversion program has been demonstrated in a number of studies (see Anderson et al., 2015; and Hunt et al., 2019) and is not reassessed here. The purpose of this section is to assess the accuracy of MFIT direct simulations through five synthetic test cases. Tests 1 and 2 address the case of a single flow channel described as a singleporosity medium in which the transport is governed by advection–dispersion. An instantaneous injection of the tracer is assumed in test 1, whereas test 2 addresses the case of an exponentially decaying concentration at the inlet. A doubleporosity medium and single flow channel are assumed in tests 3 and 4, which conform to the assumptions of the SFDM and 2RNE model, respectively. In test 5, the tracer is transported by advection–dispersion in a multiflow system composed of two channels. This scenario corresponds to the MDM. The input parameters for the five test cases are listed in Table 2. The BTCs simulated by MFIT for tests 1, 2, and 4 are compared to those obtained by CXTFIT. The MFIT simulations for tests 3 and 5 are compared against those obtained by TRACI. As shown in Fig. 2, very good agreement was obtained in each case.
The purpose of this section is to assess the automatic multistart method described in Sect. 3 using a new synthetic test case. A multimodal BTC that corresponds to three channels has been simulated using the MDMi program with the parameters listed in Table 3. A “blind” inversion of this BTC has been performed using the automatic multistart method with a maximum number of flow channels ${N}_{max}=\mathrm{6}$. The only model parameter that has been fixed prior to the inversion process was the total flow rate Q to simplify the postcomparison of the inverted mass values in each channel with the “true” mass values. Otherwise, a degree of freedom would persist for the pairs of the optimized Q and m_{j} values, i.e., multiplying or dividing these parameters by the same constant would yield the same BTCs; refer to Eq. (6). The parameters m_{j}, T_{0j}, and Pe_{j} of the different flow channels were optimized with virtually no upper and lower bound constraints (minimum and maximum allowed parameter values of $\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{10}}$ and $\mathrm{1.0}\times {\mathrm{10}}^{+\mathrm{10}}$, respectively). As shown in Fig. 3, the inverted BTCs that correspond to N=3, 4, 5, and 6 channels overlap perfectly with each other and with the original simulated BTC; as shown in Table 4, the optimized values for the parameters of the threechannel model are equal to the true parameter values.
The HES is a field research facility operated by the University of Poitiers, France. The facility consists of 32 wells that have been drilled within an overall area of 0.2 km^{2} (Fig. 4) and fully penetrate a 100 m thick confined limestone aquifer. The interwell flow and transport connectivity have been shown to be mainly related to karst conduits, 0.01–3 m in diameter, that develop preferentially within specific lithostratigraphic horizons interbedded with nonkarstified limestone units. The karstified layers may contribute to the connectivity from one well to another either directly (e.g., the wells intersect with the same karst network in a single layer) or indirectly (e.g., the wells intersect with different karst network layers that are interconnected by either a third well or a subvertical fracture); see Audouin et al. (2008) and Chatelier et al. (2011).
A large number of pumping test experiments have been conducted at the HES since 2002. As discussed in a number of studies (Delay et al., 2007, 2011; Riva et al., 2009; Bodin et al., 2012; SanchezVila et al., 2016; Le Coz et al., 2017), the drawdown responses exhibit complex behaviors, which are likely due to the strong aquifer heterogeneity induced predominantly by the presence of karst features. In addition to the pumping test experiments, a number of crossborehole tracer tests have been performed at the HES since 2011. The standard experimental protocol of HES tracer experiments can be summarized as follows:

Starting a pumping experiment and waiting for the establishment of a pseudosteadystate flow regime (i.e., stabilization of interwell piezometric head gradients) is the first step, which typically takes approximately 6 h at the HES.

Performing flow log measurements in the candidate injection well to identify the main inflow and outflow levels along the well bore is the second step.

Connecting a series of 2.5 m length and 1.5 cm inner diameter PVC pipes in the injection well, from the ground down to the tracer injection depth (usually chosen to be as close as possible to a main outflow level) is the third step. The pipeline is terminated by a 5 cm length screened cap that ensures a horizontal outflow of the tracer solution in the injection well.

Injecting a tracer solution (typically 2 L of Uranine solution at 1 g L^{−1}) in the pipe and flushing with 40 L of clean groundwater is the fourth step. The total duration of an injection is typically less than 3 min.

Monitoring the tracer BTC at the pumped well using a flowthrough fluorometer (Albillia GGUNFL22) connected to a branch pipe extending from the discharge line at ground level is the last step. The fluorometer is periodically calibrated in the laboratory with solutions of 10 and 100 µg L^{−1}.
To date, more than 70 crosswell tracer experiments have been performed at the HES. The purpose here is not to interpret each of these experiments but to pick a few examples for illustrating the application of the MFIT software. The selected data correspond to three tracer experiments that were performed in 2016 and 2017 using well M22 as the pumped well and M16, MP6, and P2 as injection wells. Figure 5 shows the experimental BTCs and a collection of calibrated MFIT curves for different numbers of channels. The selected experiments were chosen for their representativeness of the BTC shapes observed at the HES, which exhibit either a single peak followed by a more or less pronounced tailing, e.g., P2M22; overlapping doublepeak responses, e.g., M16M22; or wellmarked multimodal responses, e.g., MP6M22. The mass recovery ratios for these three tracer experiments were 58 %, 79 %, and 60 %, respectively. Note that these recovery data cannot be included in the model, because the flow structure assumption that underlies the multiflow approach (Fig. 1) implies that all the mass that enters the system flows out after a certain lapse of time. The same holds for any single or doubleporosity modeling approach based on a 1D flow assumption. For tracer tests that are performed in steady state conditions and involve nonreactive tracers, an incomplete recovery of the injected mass indicates a diverging flow structure between the injection site and the monitoring point. Unfortunately, no additional information can be obtained about this flow divergence from the tracer data only. Therefore, the total mass in a multiflow model must be consistent with the recovered tracer mass rather than the injected mass.
The model fit results shown in Fig. 5 were obtained using the multistart method discussed in Sect. 3 and only SVD as a regularization tool for the inversion. None of the model parameters were fixed, and all were optimized within realistic upper and lower limits. The optimized parameter values and their composite sensitivities at the end of the optimization process are provided in the Supplement (Table S1). Unsurprisingly, the model parameters that influence the spreading of transit or residence times in the individual flow channels, while accounting for different processes (Pe, γ, β, ψ, and ω) are sensitive to the number of channels. For instance, when comparing single with multiplechannel models, the former requires lower Pe values to compensate for the coarser description of the flow system heterogeneity (recall that the dispersion coefficient integrated in the Péclet number reflects the unresolved variability of the flow velocity below the modeling scale). The same observation holds when comparing single and doubleporosity models with the same number of flow channels, i.e., the Pe values of singleporosity models are lower than the Pe values of doubleporosity models, because part of the spreading of transit or residence times in the latter case is implicitly captured by solute mass exchanges between the mobile and immobile domains. A noticeable exception is the diffusion parameter β of the SFDM model, whose values are mostly around $\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{3}}$ h${}^{\mathrm{1}/\mathrm{2}}$. This value corresponds to the upper bound of the optimization range set for this parameter, which is based on a matrix porosity of 30 %, a molecular diffusion coefficient of $\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{9}}$ m^{2} s^{−1}, and a flowchannel aperture of $\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{2}}$ m. Beta values larger than $\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{3}}$ h${}^{\mathrm{1}/\mathrm{2}}$ would be physically unrealistic. The fact that the beta value is limited by its upper bound during the optimization process indicates that the SFDM model is not suitable for describing the HES tracer experiments, as further discussed below. All other parameters have converged to values far from their optimization bounds.
Beyond what can be visually inferred from Fig. 5, the assessment of the relative fitting performance of the different models can be analyzed through the evolution of the measurement objective function, hereafter named PHI, with respect to the number, N, of channels and/or the number, P, of optimized model parameters. Figure 6 displays the PHI(N) and PHI(P) curves summarizing the bestfitting results achieved with the multistart PEST optimization and the SCEUA_P and CMAES_P global optimization routines. A number of observations can be made from this figure. As a first remark, the SCEUA curves for the two MDP models are missing in Fig. 6. The reason is that the SCEUA_P program has no “forgive error” capability, i.e., if a set of trial parameters causes the numerical evaluation of Eq. (11) or Eq. (15) to crash, the optimization process is stopped instead of moving to a new set of parameter values. Such a forgiveness option is available in the PEST and CMAES_P programs. The next observations that can be made from Fig. 6 are that the CMAES and SCEUA curves are (i) more irregular, (ii) always above or equal to their PESTcomputed counterparts, and (iii) do not always follow the expected decreasing trend in the PHI value (meaning a better model fit) as the number of channels rises, as depicted by the PEST curves. However, it must be mentioned that the number of optimization runs was much greater for the CMAES_P and SCEUA_P programs, and various optimization options were tested (e.g., changing the upper and lower parameter bounds and logtransformation of parameters). The CMAES and SCEUA curves shown in Fig. 6 are actually the “best results” obtained after several days of computation time. It is clear that the multistart PEST optimization method performs better in each case.
The PHI curves obtained by PEST can be viewed as Pareto curves, illustrating the tradeoff between the model fitting quality and the number of channels or the number of calibration parameters. It must be noted that since no Tikhonov regularization was used in this illustration example, the model inversion results for higher N values are likely affected by overfitting. More reliable parameter values could be obtained by adding Tikhonov regularization constraints to the optimization process.
According to the PHI(N) curves shown in Fig. 6, the MDMi model and MDPSFDM perform similarly for the three tracer tests, and the related PHI(N) curves are hardly differentiable. This result was expected, as the short duration of the HES tracer tests, typically from a few hours to a few days, makes the matrix diffusion process unlikely to be significant. Assuming exponential decaying (MDMed) instead of instantaneous (MDMi) injection gives slightly better fitting results for a low number of channels but provides no benefit for a moderatetohigh number of channels. According to the PHI(N) curves, the fitting performance of the MDP2RNE model seems significantly better than that of the three other models. However, this observation must be counterbalanced by the larger number of calibration parameters in the MDP2RNE model (see Table 1). A twochannel MDP2RNE model involves 13 parameters, which corresponds to the number of parameters in a fourchannel MDMi model. The PHI(P) curves shown in Fig. 6 provide a fairer assessment of the fitting performance of the different models. According to these curves, the MDP2RNE model performs slightly better than the MDMi model for the P2M22 tracer test (single peak slightly tailed BTC), almost equally well for the M16M22 tracer test (overlapping double peaks), and worse for the MP6M22 tracer test (wellmarked multimodal BTC). It must be appreciated that these two models should not be opposed to each other. Both models likely provide an equally valid description of the tracer transport in the HES aquifer while relying on different conceptualizations of the medium heterogeneity.
The Pareto curves in Fig. 6 indicate that the final choice of a model, if one is to be made, relies on a tradeoff between the desired fitting accuracy and the desired degree of simplification or complexity with respect to the model structure (number of channels and/or number of model parameters). Beyond this subjective (expert) decision, which may depend on the goal of the study, and therefore will not be discussed further in the present application case, uncertainty remains in the inverted model parameters as a consequence of the nonuniqueness of the inverse problem. This uncertainty is related to both the equifinality of the model parameters, which is partly due to the multiflow framework structure, and the measurement noise in the tracer BTCs. Figures 7 and 8 illustrate the postcalibration uncertainty analysis capabilities of MFIT via an assessment of the MDMi and MDP2RNE model fittings of the M16M22 tracer BTC with 1, 2, and 3 flow channels. Owing to the balance between the Q and m_{j} terms in the model equation (Eqs. 6 and 15), at least one of these parameters must be fixed to assess the uncertainty of the other parameters. Here, the value of Q was set to 25 m^{3} h^{−1}, which ensures the consistency of the model against the recovered tracer mass that was independently calculated from the experimental data (refer to Table S1). Following the PEST optimization of the different model parameters, 500 calibrationconstrained parameter fields were stochastically generated and recalibrated by PEST. Depending on the model (MDMi or MDP2RNE) and number of flow channels, between 483 and 500 recalibration runs successfully achieved a level of fit that is fairly similar (i.e., within a tolerance of +5 % for the PHI value; refer to Sect. 3) to that associated with the original calibration parameter field. The histograms shown in Figs. 7 and 8 were constructed from these recalibrated parameter sets and illustrate the multitude of parameter combinations that are equally good for a given number of flow channels, in terms of fitting the M16M22 tracer BTC. As shown in these figures, the confidence intervals are quite narrow for most parameters but tend to widen as the number of channels increases, which reflects the equifinality of the multiflow modeling approach. Although not shown here, it has been established that the tailed behaviors of the parameters L_{j} and ω_{j} in Fig. 8 are due to a partial correlation between these two parameters (refer to Eq. 15), i.e., fixing the value of one parameter prior to the inversion drastically reduces the uncertainty of the other parameter. As previously discussed, the higher Pe values in Fig. 8 compared to Fig. 7 are due to the fact that the distribution of the transit or residence times with the 2RNE model is primarily controlled by the solute mass exchanges between the mobile and immobile domains.
Multiple flow path transport is likely the rule rather than the exception in most transport problems in fractured and karst aquifers. The main aim of this paper was to present a new curvefitting tool for the analytical modeling of BTCs from tracer tests performed in such media. The MFIT software is a free opensource Windowsbased GUI that provides access to four multiflow transport models. The multiflow approach assumes that the transport from the injection site to the monitoring point takes place in a number of independent 1D channels. The channels are not assumed to represent individual fractures or karst conduits but are lumped submodels of the main flow routes used by the tracer through the fractures or karst conduit network. The multiflow modeling framework allows for the simulation of multimodal BTCs, which are frequently observed in fractured and karst aquifers. Two of the MFIT transport models combine the multiflow framework and the doubleporosity concept, which is applied at the scale of the individual channels. This modeling approach, which has been named MDP, is believed to be new and versatile for the fitting of BTCs with multiple local peaks and/or extensive backward tailing. The accuracy of the MFITcomputed BTCs was verified against two other wellaccepted simulation tools for five synthetic test cases.
An important feature of MFIT is its compatibility and interface with the advanced calibration tools of the PEST suite of programs. Hence, MFIT is the first BTC fitting tool that allows for regularized inversion and nonlinear analysis of the postcalibration uncertainty of model parameters. Given the nonlinearity of the MFIT model equations, an original multistart algorithm was implemented to maximize the chances for PEST to converge to the global optimal solution in the parameter space during a BTC fitting procedure. The main drawback of the multistart optimization method is that the processing time can be long (up to a few hours) if a large number of channels is assumed in the model. Time reduction for this method is one of the development perspectives of the MFIT code, as the multistart process is computationally parallelizable. Other development perspectives are the management of more complex injection signals, e.g., described as multiple steps, and the implementation of additional analytical transport models for the simulation of reactive transport processes.
Three tracer test BTCs from the HES in Poitiers, France, were used for illustrating the application of the MFIT software. An analysis of the Pareto curves between the model fitting quality and the number of model calibration parameters suggests that the MDMi and MDP2RNE models are the most appropriate for the interpretation of HES tracer tests. This preliminary result needs to be refined or confirmed by the analysis of additional HES tracer BTCs.
The source codes of the MFIT program suite version 1.0.0 are available from https://doi.org/10.5281/zenodo.3470751 (Bodin, 2020) under the terms of the CeCILL Free Software License Agreement v2.1 (https://spdx.org/licenses/CECILL2.1.html, last access: 24 June 2020). An “EXE” installation package compiled with Inno Setup (http://www.jrsoftware.org/isinfo.php, last access: 24 June 2020) and a user's guide are provided along with the source codes. The following numerical libraries are required for the compilation of the MFIT suite of codes: Boost (https://www.boost.org/, last access: 24 June 2020), GSLGNU (https://www.gnu.org/software/gsl/, last access: 24 June 2020) and Spline (https://github.com/ttk592/spline, last access: 24 June 2020). The PEST program package is also required for running MFIT. PEST is distributed by default using the MFIT software installer or can be independently downloaded from http://www.pesthomepage.org/Downloads.php (last access: 24 June 2020). The data of the HES tracer experiments processed in Sect. 6 of this study are available from the H+ database (http://hplus.ore.fr/en/poitiers/datapoitiers, last access: 24 June 2020, de Dreuzy et al., 2020) with registration of a free account.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1329052020supplement.
The author declares that there is no conflict of interest.
The HES field research facility is managed by Gilles Porel. The experimental protocol used for the HES tracer experiments was developed by Benoit Nauleau and Gilles Porel. The assistance of Benoit Nauleau, Gilles Porel, and Denis Paquet in conducting the tracer experiments is gratefully acknowledged. The author would like to thank the two anonymous reviewers for their valuable comments and feedback, which helped with improving the article.
This research was supported by the French National Observatory H+, the European Union (ERDF), and Région Nouvelle Aquitaine.
This paper was edited by Jeffrey Neal and reviewed by two anonymous referees.
Abdelaziz, R., Merkel, B. J., ZambranoBigiarini, M., and Nair, S.: Particle swarm optimization for the estimation of surface complexation constants with the geochemical model PHREEQC3.1.2, Geosci. Model Dev., 12, 167–177, https://doi.org/10.5194/gmd121672019, 2019.
Anderson, M. P., Woessner, W. W., and Hunt, R. J.: Applied groundwater modeling: simulation of flow and advective transport, Academic Press, 2015.
Arsenault, R., Poulin, A., Cote, P., and Brissette, F.: Comparison of stochastic optimization algorithms in hydrological model calibration, J. Hydrol. Eng., 19, 1374–1384, https://doi.org/10.1061/(ASCE)HE.19435584.0000938, 2014.
Audouin, O., Bodin, J., Porel, G., and Bourbiaux, B.: Flowpath structure in a limestone aquifer: multiborehole logging investigations at the hydrogeological experimental site of Poitiers, France, Hydrogeol. J., 16, 939–950, https://doi.org/10.1007/s1004000802754, 2008.
Berkowitz, B., Cortis, A., Dentz, M., and Scher, H.: Modeling nonFickian transport in geological formations as a continuous time random walk, Rev. Geophys., 44, RG2003, https://doi.org/10.1029/2005RG000178, 2006.
Bertrand, G., CelleJeanton, H., Huneau, F., Baillieux, A., Mauri, G., Lavastre, V., Undereiner, G., Girolami, L., and Moquet, J. S.: Contaminant transfer and hydrodispersive parameters in basaltic lava flows: artificial tracer test and implications for longterm management, Open Geosci., 7, 513–526, https://doi.org/10.1515/geo20150037, 2015.
Bodin, J.: MFIT 1.0.0, Zenodo, https://doi.org/10.5281/zenodo.3470751, 2020.
Bodin, J., Porel, G., and Delay, F.: Simulation of solute transport in discrete fracture networks using the time domain random walk method, Earth Planet. Sc. Lett., 208, 297–304, https://doi.org/10.1016/S0012821X(03)000529, 2003a.
Bodin, J., Delay, F., and de Marsily, G.: Solute transport in a single fracture with negligible matrix permeability: 2. mathematical formalism, Hydrogeol. J., 11, 434–454, https://doi.org/10.1007/s1004000302691, 2003b.
Bodin, J., Ackerer, P., Boisson, A., Bourbiaux, B., Bruel, D., de Dreuzy, J.R., Delay, F., Porel, G., and Pourpak, H.: Predictive modelling of hydraulic head responses to dipole flow experiments in a fractured/karstified limestone aquifer: Insights from a comparison of five modelling approaches to realfield experiments, J. Hydrol., 454, 82–100, https://doi.org/10.1016/j.jhydrol.2012.05.069, 2012.
Boon, M., Bijeljic, B., and Krevor, S.: Observations of the impact of rock heterogeneity on solute spreading and mixing, Water Resour. Res., 53, 4624–4642, https://doi.org/10.1002/2016WR019912, 2017.
Chatelier, M., Ruelleu, S., Bour, O., Porel, G., and Delay, F.: Combined fluid temperature and flow logging for the characterization of hydraulic structure in a fractured karst aquifer, J. Hydrol., 400, 377–386, 2011.
Coats, K. H. and Smith, B. D.: Dead end pore volume and dispersion in porous media, Soc. Pet. Eng. J., 23, 73–84, 1964.
de Dreuzy, J.R., Bodin, J., Le Grand, H., Davy, P., Boulanger, D., Battais, A., Bour, O., Gouze, P. and Porel, G.: H+ database, available at: http://hplus.ore.fr/en/poitiers/datapoitiers, last access: 24 June 2020.
Delay, F., Kaczmaryk, A., and Ackerer, P.: Inversion of interference hydraulic pumping tests in both homogeneous and fractal dual media, Adv. Water Resour., 30, 314–334, 2007.
Delay, F., Ackerer, P., and Guadagnini, A.: Theoretical analysis and field evidence of reciprocity gaps during interference pumping tests, Adv. Water Resour., 34, 592–606, https://doi.org/10.1016/j.advwatres.2011.02.006, 2011.
Dentz, M., Le Borgne, T., Englert, A., and Bijeljic, B.: Mixing, spreading and reaction in heterogeneous media: A brief review, J. Contam. Hydrol., 120–121, 1–17, https://doi.org/10.1016/j.jconhyd.2010.05.002, 2011.
Diersch, H.J.: FEFLOW: finite element modeling of flow, mass and heat transport in porous and fractured media, SpringerVerlag, Berlin Heidelberg, https://doi.org/10.1007/9783642387395, 2014.
Doherty, J.: Calibration and uncertainty analysis for complex environmental models, Watermark Numer. Comput. Brisb. Aust., 2015.
Doherty, J.: PEST, modelindependent parameter estimation – User manual part I: PEST, SENSAN and global optimisers, Watermark Numer. Comput. Brisb. Aust., 2019a.
Doherty, J.: PEST, modelindependent parameter estimation – User manual part II: PEST utility support software, Watermark Numer. Comput. Brisb. Aust., 2019b.
Doherty, J. E., Hunt, R. J., and Tonkin, M. J.: Approaches to highly parameterized inversion: A guide to using PEST for modelparameter and predictiveuncertainty analysis, U.S. Geological Survey Scientific Investigations Report 2010–5211, 71 pp., 2010.
Duan, Q., Sorooshian, S., and Gupta, V.: Effective and efficient global optimization for conceptual rainfallrunoff models, Water Resour. Res., 28, 1015–1031, https://doi.org/10.1029/91WR02985, 1992.
Espinet, A. J. and Shoemaker, C. A.: Comparison of optimization algorithms for parameter estimation of multiphase flow models with application to geological carbon sequestration, Adv. Water Resour., 54, 133–148, https://doi.org/10.1016/j.advwatres.2013.01.003, 2013.
Fang, Q., Ma, L., Harmel, R. D., Yu, Q., Sima, M. W., Bartling, P. N. S., Malone, R. W., Nolan, B. T., and Doherty, J.: Uncertainty of CERESMaize calibration under different irrigation strategies using PEST optimization algorithm, Agronomy, 9, 241, https://doi.org/10.3390/agronomy9050241, 2019.
Field, M. S. and Leij, F. J.: Solute transport in solution conduits exhibiting multipeaked breakthrough curves, J. Hydrol., 440–441, 26–35, https://doi.org/10.1016/j.jhydrol.2012.03.018, 2012.
Field, M. S. and Leij, F. J.: Combined physical and chemical nonequilibrium transport model for solution conduits, J. Contam. Hydrol., 157, 37–46, https://doi.org/10.1016/j.jconhyd.2013.11.001, 2014.
Gaudard, A., Schwefel, R., Vinnå, L. R., Schmid, M., Wüest, A., and Bouffard, D.: Optimizing the parameterization of deep mixing and internal seiches in onedimensional hydrodynamic models: a case study with Simstrat v1.3, Geosci. Model Dev., 10, 3411–3423, https://doi.org/10.5194/gmd1034112017, 2017.
Gharasoo, M., Wietzke, L. M., Knorr, B., Bakkour, R., Elsner, M., and Stumpp, C.: A robust optimization technique for analysis of multitracer experiments, J. Contam. Hydrol., 224, 103481, https://doi.org/10.1016/j.jconhyd.2019.04.004, 2019.
Goldscheider, N., Meiman, J., Pronk, M., and Smart, C.: Tracer tests in karst hydrogeology and speleology, Int. J. Speleol., 37, 27–40, https://doi.org/10.5038/1827806X.37.1.3, 2008.
Gutierrez, A., Klinka, T., Thiery, D., Buscarlet, E., Binet, S., Jozja, N., Defarge, C., Leclerc, B., Fecamp, C., Ahumada, Y., and Elsass, J.: TRAC, a collaborative computer tool for tracertest interpretation, in: EPJ Web of Conferences, vol. 50, no. 03002, EDP Sciences, 2013.
Guvanasen, V. and Guvanasen, V. M.: An approximate semianalytical solution for tracer injection tests in a confined aquifer with a radially converging flow field and finite volume of tracer and chase fluid, Water Resour. Res., 23, 1607–1619, https://doi.org/10.1029/WR023i008p01607, 1987.
Hadamard, J.: Sur les problèmes aux dérivées partielles et leur signification physique, Princet. Univ. Bull., 49–52, 1902.
Hansen, N. and Ostermeier, A.: Completely derandomized selfadaptation in evolution strategies, Evol. Comput., 9, 159–195, https://doi.org/10.1162/106365601750190398, 2001.
Hunt, R. J., Fienen, M. N., and White, J. T.: Revisiting “An Exercise in Groundwater Model Calibration and Prediction” after 30 years: Insights and New Directions, Groundwater, 58, 168–182, https://doi.org/10.1111/gwat.12907, 2019.
Käss, W.: Evaluation and interpretation of tracing tests, in: Tracing technique in geohydrology, Balkema, Rotterdam, the Netherlands, 372–379, 1998.
Käss, W.: Geohydrologische Markierungstechnik, Borntraeger, Stuttgart, Germany, 2004.
Kreft, A. and Zuber, A.: On the physical meaning of the dispersion equation and its solutions for different initial and boundary conditions, Chem. Eng. Sci., 33, 1471–1480, https://doi.org/10.1016/00092509(78)851963, 1978.
Labat, D. and Mangin, A.: Transfer function approach for artificial tracer test interpretation in karstic systems, J. Hydrol., 529, 866–871, https://doi.org/10.1016/j.jhydrol.2015.09.011, 2015.
Langevin, C. D., Hughes, J. D., Banta, E. R., Niswonger, R. G., Panday, S., and Provost, A. M.: Documentation for the MODFLOW 6 Groundwater Flow Model, USGS Numbered Series, U.S. Geological Survey, Reston, VA, 2017.
Le Coz, M., Bodin, J., and Renard, P.: On the use of multiplepoint statistics to improve groundwater flow modeling in karst aquifers: A case study from the Hydrogeological Experimental Site of Poitiers, France, J. Hydrol., 545, 109–119, https://doi.org/10.1016/j.jhydrol.2016.12.010, 2017.
Leibundgut, C., Maloszewski, P., and Külls, C.: Mathematical modelling of experimental data: Artificial tracer experiments in multiflow systems, in: Tracers in hydrology, WileyBlackwell, 140–144, 2009.
Loefgren, M., Crawford, J., and Elert, M.: Tracer tests – possibilities and limitations. Experience from SKB fieldwork: 1977–2007, Swedish Nuclear Fuel and Waste Management Co., available at: http://inis.iaea.org/Search/search.aspx?orig_q=RN:39028075 (last access: 24 June 2020), 2007.
Long, A. J.: RRAWFLOW: RainfallResponse Aquifer and Watershed Flow Model (v1.15), Geosci. Model Dev., 8, 865–880, https://doi.org/10.5194/gmd88652015, 2015.
Maloszewski, P. and Zuber, A.: Mathematical modeling of tracer behaviour in shortterm experiments in fissured rocks, Water Resour. Res., 26, 1517–1528, 1990.
Maloszewski, P., Harum, T., and Benischke, R.: Mathematical modelling of tracer experiments in the karst of Lurbach system, Steirische Beitraege Zur Hydrogeol., 43, 116–136, 1992.
Marino, M. A.: Distribution of contaminants in porous media flow, Water Resour. Res., 10, 1013–1018, https://doi.org/10.1029/WR010i005p01013, 1974.
Massei, N., Wang, H. Q., Field, M. S., Dupont, J. P., Bakalowicz, M., and Rodet, J.: Interpreting tracer breakthrough tailing in a conduitdominated karstic aquifer, Hydrogeol. J., 14, 849–858, https://doi.org/10.1007/s1004000500103, 2006.
Moore, C. and Doherty, J.: The cost of uniqueness in groundwater model calibration, Adv. Water Resour., 29, 605–623, https://doi.org/10.1016/j.advwatres.2005.07.003, 2006.
Moreno, L. and Tsang, C. F.: Multiplepeak response to tracer injection tests in single fractures: a numerical study, Water Resour. Res., 27, 2143–2150, https://doi.org/10.1029/91WR00507, 1991.
Neuman, S. P. and Tartakovsky, D. M.: Perspective on theories of nonFickian transport in heterogeneous media, Adv. Water Resour., 32, 670–680, https://doi.org/10.1016/j.advwatres.2008.08.005, 2009.
Parker, J. and van Genuchten, M. T.: Determining transport parameters from laboratory and field tracer experiments, Va. Agric. Exp. Stn. Bull., 843, 1–96, 1984.
Piotrowski, A. P. and Napiorkowski, J. J.: Optimizing neural networks for river flow forecasting – Evolutionary Computation methods versus the LevenbergMarquardt approach, J. Hydrol., 407, 12–27, https://doi.org/10.1016/j.jhydrol.2011.06.019, 2011.
Riva, M., Guadagnini, A., Bodin, J., and Delay, F.: Characterization of the Hydrogeological Experimental Site of Poitiers (France) by stochastic well testing analysis, J. Hydrol., 369, 154–164, https://doi.org/10.1016/j.jhydrol.2009.02.040, 2009.
Runkel, R. L.: Onedimensional transport with inflow and storage (OTIS): a solute transport model for streams and rivers, U.S. Geological Survey WaterResources Investigations Report 984018, https://doi.org/10.3133/wri984018, 1998.
SanchezVila, X., Ackerer, P., Delay, F., and Guadagnini, A.: Characterization of reciprocity gaps from interference tests in fractured media through a dual porosity model, Water Resour. Res., 52, 1696–1704, https://doi.org/10.1002/2015WR018171, 2016.
Sauty, J. P., Kinzelbach, W., and Voss, A.: CATTI: Computer aided tracer test interpretation, BRGM Rep., 1992.
SiirilaWoodburn, E. R., SanchezVila, X., and FernandezGarcia, D.: On the formation of multiple local peaks in breakthrough curves, Water Resour. Res., 51, 2128–2152, https://doi.org/10.1002/2014WR015840, 2015.
Singh, S. K., Liang, J., and Bárdossy, A.: Improving the calibration strategy of the physicallybased model WaSiMETH using critical events, Hydrol. Sci. J., 57, 1487–1505, https://doi.org/10.1080/02626667.2012.727091, 2012.
Skahill, B. E. and Doherty, J.: Efficient accommodation of local minima in watershed model calibration, J. Hydrol., 329, 122–139, https://doi.org/10.1016/j.jhydrol.2006.02.005, 2006.
Streetly, H. R., Hamilton, A. C. L., Betts, C., Tellam, J. H., and Herbert, A. W.: Reconnaissance tracer tests in the Triassic sandstone aquifer north of Liverpool, UK, Q. J. Eng. Geol. Hydrogeol., 35, 167–178, https://doi.org/10.1144/14709236/200030, 2002.
Tikhonov, A. N. and Arsenin, V. Y.: Solutions of illposed problems, Winston & Sons, Washington, 1977.
Tinet, A.J., Collon, P., Philippe, C., Dewaide, L., and Hallet, V.: OMMADE: An opensource program to simulate onedimensional solute transport in multiple exchanging conduits and storage zones, Comput. Geosci., 127, 23–35, https://doi.org/10.1016/j.cageo.2019.03.001, 2019.
Toride, N., Leij, F. L., and van Genuchten, M. T.: A comprehensive set of analytical solutions for nonequilibrium solute transport with firstorder decay and zeroorder production, Water Resour. Res., 29, 2167–2182, 1993.
Toride, N., Leij, F. J., and Van Genuchten, M. Th.: The CXTFIT code for estimating transport parameters from laboratory or field tracer experiments, version 2.1, Res. Rep. No. 137 U Salin. Lab. USDA ARS Riverside CA, 1999.
Tsang, C. F. and Neretnieks, I.: Flow channeling in heterogeneous fractured rocks, Rev. Geophys., 36, 275–298, 1998.
van Genuchten, M. T., Simunek, J., Leij, F. J., Toride, N., and Sejna, M.: Stanmod: Model Use, Calibration, and Validation, Trans. ASABE, 55, 1353–1366, 2012.
Wang, J., Wang, C., Rao, V., Orr, A., Yan, E., and Kotamarthi, R.: A parallel workflow implementation for PEST version 13.6 in highperformance computing for WRFHydro version 5.0: a case study over the midwestern United States, Geosci. Model Dev., 12, 3523–3539, https://doi.org/10.5194/gmd1235232019, 2019.
Woodward, S. J. R., Wöhling, T., and Stenger, R.: Uncertainty in the modelling of spatial and temporal patterns of shallow groundwater flow paths: The role of geological and hydrological site information, J. Hydrol., 534, 680–694, https://doi.org/10.1016/j.jhydrol.2016.01.045, 2016.
Worthington, S. R. H. and Ford, D. C.: Selforganized permeability in carbonate aquifers, Ground Water, 47, 326–336, https://doi.org/10.1111/j.17456584.2009.00551.x, 2009.
Yang, P., Li, Y., Groves, C., and Hong, A.: Coupled hydrogeochemical evaluation of a vulnerable karst aquifer impacted by septic effluent in a protected natural area, Sci. Total Environ., 658, 1475–1484, https://doi.org/10.1016/j.scitotenv.2018.12.172, 2019.
Zhang, Y., Benson, D. A., and Reeves, D. M.: Time and space nonlocalities underlying fractionalderivative models: Distinction and literature review of field applications, Adv. Water Resour., 32, 561–581, https://doi.org/10.1016/j.advwatres.2009.01.008, 2009.
Zheng, C., Hill, M. C., Cao, G., and Ma, R.: MT3DMS: model use, calibration, and validation, Trans. ASABE, 55, 1549–1559, https://doi.org/10.13031/2013.42263, 2012.
Zheng, C. and Bennett, G. D.: Applied contaminant transport modeling, 2nd edn., WileyInterscience, New York, 2002.
Zhou, H., GómezHernández, J. J., and Li, L.: Inverse methods in hydrogeology: Evolution and recent trends, Adv. Water Resour., 63, 22–37, https://doi.org/10.1016/j.advwatres.2013.10.014, 2014.
Zuber, A.: Theoretical possibilities of the twowell pulse method, Isot. Tech. Groundw. Hydrol. Vol. II IAEA, 277–294, 1974.
 Abstract
 Introduction
 Multimodal and heavytailed BTCs: causes and modeling
 Code implementation and inversion
 Code verification
 Assessment of the multistart optimization method
 Application example: analysis of tracer BTCs from the Hydrogeological Experimental Site in Poitiers, France
 Summary and conclusions
 Appendix A: Glossary
 Code and data availability
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Multimodal and heavytailed BTCs: causes and modeling
 Code implementation and inversion
 Code verification
 Assessment of the multistart optimization method
 Application example: analysis of tracer BTCs from the Hydrogeological Experimental Site in Poitiers, France
 Summary and conclusions
 Appendix A: Glossary
 Code and data availability
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement