# Solute transport in dual conduit structure: Effects of aperture and flow rate

Chaoqi Wang, Samer Majdalani, Vincent Guinot, Hervé Jourde

► To cite this version:

Chaoqi Wang, Samer Majdalani, Vincent Guinot, Hervé Jourde. Solute transport in dual conduit structure: Effects of aperture and flow rate. Journal of Hydrology, 2022, 613, Part A, pp.128315. 10.1016/j.jhydrol.2022.128315. hal-03751970

HAL Id: hal-03751970

https://hal.science/hal-03751970v1

Submitted on 25 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 **Solute transport in dual conduit structure: effects of aperture and flow rate**

2 Chaoqi Wang¹, Samer Majdalani¹, Vincent Guinot¹,², Hervé Jourde¹

3 1. HSM, Univ. Montpellier, IRD, CNRS, Montpellier, France

4 2. Inria Lemon, Univ. Montpellier, Montpellier, France

5 Abstract

6 A typical karst aquifer configuration is the multiple conduit structure. However, it remains to be investigated
7 how the aperture distributions and the flow rate should influence the transport process in the multiple conduit
8 structures.

9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
To better understand the transport process in the multiple conduit structures, 11 lab-scale dual-conduit
structures are manufactured by varying the apertures of the two conduits (h₁ and h₂ denote the aperture width of the
shorter conduit and the longer conduit respectively). Solute transport experiments of three different flow rates are
conducted on these structures. As the flow rate increases, the dual-conduit structures are more likely to present
dual-peaked BTCs. The 11 structures make one exhaustive representation of the possible aperture combination of
the dual-conduit structures and the transport experiments have been conducted by three flow rates (varying by 2
degrees of magnitude), so the experimental results constitute a detailed material that should improve the
understanding of transport processes in such structures.

Two numerical models, Weighted Sum Advection–Dispersion Equation (WSADE) and Dual Region Mobile Immobile Model (DRMIM), are applied to fit the experimental BTCs in order to obtain some insight into the actual solute-transport processes by exploring the calibrated model parameters. Considering the possible effect of solute detention, we initially applied the DRMIM model. This DRMIM better replicated the experimental BTCs than the WSADE. This study suggests the karst community shall take the DRMIM as one candidate transport model for characterizing the dual-peaked BTCs obtained in karst aquifers.

23 **Keywords:** Dual-conduit structure, Flow rate, Aperture width, Dual-peaked BTC
24

**Highlights:**

25

1. This study provides experiment material to enhance understanding of transport processes in dual-conduit structures.

26

27

2. For dual-conduit structures, the influence of aperture combination and flow rate are investigated.

28

3. Experimental BTCs were well simulated by the DRMIM transport model.
29

# 1. Introduction

30

Artificial tracing is a useful tool for investigating karst aquifers (Goldscheider et al., 2008a). Injecting a known quantity of tracer at an upstream position, and monitoring the variations in the tracer concentration at a downstream location yields a BreakThrough Curve (BTC). BTCs can bring information on underground karst geometry (Wang et al., 2020; Zhao et al., 2021) and effective transport parameters (Goldscheider and Drew, 2007).

31

32

33

34

Laboratory experiments have been widely applied to investigate the flow and transport processes through complex karst aquifers. The major advantage of lab-scale experiments over field experiments is their flexibility, in that they allow for a flexible selection of sampling points, easy implementation of different hydraulic conditions, as well as detailed knowledge of the geometry of the studied media (Mohammadi et al., 2021).

35

36

37

38

According to Mohammadi et al. (2021), the common laboratory karst models can be categorized into three groups: rock block, pipe/fracture network, and pipe-matrix coupling model. For the rock block model, artificial fractures or conduits can be created on a rock block to make a physical model. The rock block model has supported various research works, including hydraulic tomography (Brauchler et al., 2003, 2013; Sharmeen et al., 2012), investigation of the flow process (Anaya et al., 2014; Cherubini et al., 2012; Develi & Babadagli, 2015; DiFrenna et al., 2008; Leven et al., 2004), investigation of solute (Cherubini et al., 2012) and heat transport process (Cherubini et al., 2017; Pastore et al., 2015). Concerning the pipe/fracture network model, researchers have manufactured single-fracture models (Qian et al., 2007, 2011; Tzelepis et al., 2015), fracture networks (Hull et al., 1987; Karay & Hajnal, 2015) and conduit networks (Anger & Alexander, 2013; Field & Leij, 2012; Wu et al., 2019; Zhao et al., 2019, 2021). The experiments on these lab-scale models have revealed some fundamental mechanisms controlling the flow and transport process in karst media. For example, Qian et al. (2011) found that the solute transport process in a single fracture exhibit BTCs with heavy tailing; Field and Leij (2012) found that three karst geometries may exhibit multi-peaked BTCs: an auxiliary conduit is connected to the main conduit, waterfall, and pools. While the effect of matrix is not considered in the pipe/fracture network models. Several

39

40

41

42

43

44

45

46

47

48

49

50

51
studies have focused on the effect of the transport exchange between conduits and matrix. To this end, they have used pipe-matrix coupling models (Castro, 2017; Faulkner et al., 2009; Florea & Wicks, 2001; Gallegos et al., 2013; Li et al., 2008; Mohammadi et al., 2019; Xiao et al., 2018). For example, Li et al. (2008) have studied contaminant retention processes in the matrix of pipe-matrix coupling models; Mohammadi et al. (2019) built bench-scale karst models and investigated the influence of hydraulic gradient and cave pattern on tracing test BTCs.

The classical Advection-Dispersion-Equation (ADE) (Bear, 1988; Majdalani et al., 2015) model often fails to reproduce the BTCs of the tracing tests (Goldscheider et al., 2008b; Hauns et al., 2001; Massei et al., 2006; Moreno & Tsang, 1991; Perrin & Luetscher, 2008). The ADE is valid for the BTCs obtained in homogeneous porous media under some conditions (Berkowitz et al., 2006). However, all aquifers are heterogeneous to different degrees (Bakalowicz, 2005; Ghodrati & Jury, 1992; Levy & Berkowitz, 2003) and among them, fractured and karst aquifers are particularly heterogeneous (Bakalowicz, 2005). Heterogeneous media entail ‘anomalous” or “non-Fickian” transport behaviors, which may exhibit BTCs with early arrivals and late time tails or multi-peaked BTCs. These particularities cannot be well represented by the ADE model.

A variety of models accounting for anomalous transport have been developed. Mobile-Immobile Model (MIM) (van Genuchten & Wierenga, 1976) and Transient Storage Model (TSM) (Runkel, 1998) have a similar structure: one mobile region (the main channel) governed by ADE and one immobile (stagnant) region; the two regions exchange mass due to concentration difference. The TSM has been applied to reproduce BTCs from solute transport experiments in lab-scale karst models (Dewaide et al., 2016; Zhao et al., 2019, 2021). The multi-region ADE model (Majdalani et al., 2015, 2018) consists of two or more ADE regions flowing in parallel and these regions exchange mass due to concentration differences. These multi-region ADE models were successfully used to model lab-scale solute transport experiments (Field & Leij, 2012; Majdalani et al., 2015, 2018; Wang et al., 2020). The so-called Weighted-Sum ADE (WSADE) model (Field & Leij, 2012) also considers several parallel ADE regions while the mass exchange between the regions is not considered. Dual-peaked BTCs can be reproduced by
75

the two-region ADE model (Wang et al., 2020) or the WSADE model (Field & Leij, 2012).

76

The tracing BTCs can be influenced by different properties of karst aquifers. Via numerical simulation of solute transport process in a conceptual karst model of dual flow system (i.e. conduit and matrix), Peely et al. (2021) have studied the effect of six properties on the BTCs, including (i) tracing distance and karst tortuosity, (ii) hydraulic gradient, (iii) matrix hydraulic conductivity, (iv) exchange flow, (v) conduit diameter and (vi) conduit roughness height. For example, increasing tracing distance causes the peak concentration ($C_p$) to decrease and the peak time ($T_p$) to become longer; increasing tortuosity also causes $C_p$ to decrease and $T_p$ to lengthen. For another example, increasing the exchange coefficient (the degree of hydraulic connection between the conduit and the host matrix) causes $C_p$ to increase and $T_p$ to shorten. The diameter variation of a single conduit may cause BTCs with strong tailing (Hauns et al., 2001; Zhao et al., 2019, 2021). The different karst network patterns should influence the BTCs' shape (Mohammadi et al., 2019), and multi-peaked BTCs are frequently observed in karst aquifers. In a previous study (Wang et al., 2020), we carried out solute transport experiments on different dual-conduit karst configurations and investigated how the dual-peaked BTCs can be influenced by some properties of the dual-conduit structure. Specifically, increasing the length ratio causes the two peaks to be more separated; increasing the total length causes the two peaks to be more separated and the concentration value of both peaks to be lower; the connection angles between the divergence conduits and the convergence conduit also influence the BTCs.

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

However, the study presented by Wang et al. (2020) has two main limitations. First, we assumed an ideal situation for all the structures: the two conduits have the same diameter. Since the objective of (Wang et al., 2020) is to study the effect of length and connection angle of the conduits, this ideal situation does not harm the previous experimental conclusions. According to the field observations, the different conduits usually bear different diameters. It is unclear how the different conduit diameter combinations may influence the transport process. Thus, it is important to study the structures that have different conduit diameter combinations in order to better understand the transport process in such dual-conduit structures. Second, the previous transport experiments have

92

93

94

95

96

97
98

99

been conducted under a constant flow rate, while in natural karst aquifers flow rate is not constant. It is also unclear
how the flow rate may influence the transport process.

100

101

102

103

104

105

106

107

108

The objective of this paper is to answer two research questions that arose from the one-dimensional
experiments reported in Wang et al. (2020), namely the influence of conduits aperture and flow rate on the
transport process in dual-conduit structures. 11 bench-scale physical configurations with different aperture
combinations are built. Solute transport experiments are conducted under three different flow rates. Then, two
transport models are applied in this study: the Weighted Sum Advection–Dispersion Equation (WSADE) model
and the Dual-Region Mobile-Immobile model (DRMIM). We use both models to fit the observed BTCs and check
how the different transport processes may influence the calibrated model parameters. Along with Wang et al.
(2020), the purpose of the present paper is to enrich and complement the available experimental and modeling
material on the transport process in dual-conduit structures.
109

## 2. Experimental setup

110

### 2.1 Dual-conduit structures

111

112

113

114

The dual-conduit experimental models are made of transparent PVC boards assembled in three layers. The top and bottom layers (Figure 1a-c) are the same: a single PVC board (150 × 25 × 2 cm). The middle layer (Figure 1a-c) includes five boards of various sizes with a thickness of 0.5 cm. The three layers are assembled (Figure 1d) with transparent silicone. The void in the middle layer serves as the space for fluid flow in this physical model.

<figure><img src="image_2.png" /><figcaption>Figure 1. Geometry of dual-conduit structure (All lengths in cm).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. Void in the physical model (All lengths in cm).</figcaption></figure>

115

116

117

118

119

120

121

We make 11 different physical configurations by varying the position of the central board (marked by the red dot in Figure 1b) in the middle layer. The apertures of the main and parallel conduits are denoted respectively by $h_1$ and $h_2$ (Figure 2). $h_1$ and $h_2$ are varied by ± 1 cm increments such that $h_1 + h_2 = 10$ cm for all configurations (Table
1) leaving the total void volume $V_m = 675 \text{ cm}^3$ unchanged. The various geometric configurations are labeled in the form “$h_1-h_2$” (e.g. Configuration 4-6 stands for $h_1 = 4 \text{ cm}$ and $h_2 = 6 \text{ cm}$). For the central board in the middle layer, the spacing from the longitudinal axis is denoted as $M$ ($M = h_2 - 5 \text{ cm} = 5 \text{ cm} - h_1$). The 11 configurations are divided into two groups: Group 1 includes Configurations 1-6, the shorter conduit has a smaller aperture than the longer conduit ($h_1 \le h_2$); Group 2 includes Configurations 7-11, the shorter conduit has a larger aperture ($h_1 > h_2$).

With these configurations, we focus on the effect of conduit apertures, while the other properties of the dual-conduit structures remain constant. First, the lengths of the dual-conduits are constant for these configurations: 105 cm for the shorter conduit and 125 cm for the longer conduit. Second, the connection angle combinations, $\theta_1 - \theta_2$, of the 11 configurations are considered to be the same: $\theta_1 = 0^\circ$ and $\theta_2 = 90^\circ$ (Figure 2).

<table><caption>Table 1. Aperture setting of physical models.</caption><thead><tr><th>Group</th><th>Configurations</th><th>Label</th><th>$h_1$ (cm)</th><th>$h_2$ (cm)</th><th>$M$ (cm)</th><th>General remark</th></tr></thead><tbody><tr><td rowspan="6">1</td><td>1</td><td>5-5</td><td>5</td><td>5</td><td>0</td><td rowspan="6">$h_1 \le h_2$</td></tr><tr><td>2</td><td>4-6</td><td>4</td><td>6</td><td>1</td></tr><tr><td>3</td><td>3-7</td><td>3</td><td>7</td><td>2</td></tr><tr><td>4</td><td>2-8</td><td>2</td><td>8</td><td>3</td></tr><tr><td>5</td><td>1-9</td><td>1</td><td>9</td><td>4</td></tr><tr><td>6</td><td>0-10</td><td>0</td><td>10</td><td>5</td></tr><tr><td rowspan="5">2</td><td>7</td><td>6-4</td><td>6</td><td>4</td><td>-1</td><td rowspan="5">$h_1 > h_2$</td></tr><tr><td>8</td><td>7-3</td><td>7</td><td>3</td><td>-2</td></tr><tr><td>9</td><td>8-2</td><td>8</td><td>2</td><td>-3</td></tr><tr><td>10</td><td>9-1</td><td>9</td><td>1</td><td>-4</td></tr><tr><td>11</td><td>10-0</td><td>10</td><td>0</td><td>-5</td></tr></tbody></table>

We have applied a rather large scale for this experimental model ($150 \times 15 \times 0.5 \text{ cm}$), aiming at reducing the error in the model geometry, which can be inevitable to a certain degree in the manufacturing process. The error in the aperture of the conduits may cause especially detrimental damage to the experimental results because the flow resistance of a single conduit is highly sensitive to the aperture (Barnard et al., 1966). To reduce the effect caused by the possible aperture error, we apply as large a scale as possible for the experimental model. On a large scale, the influence caused by the aperture error would be relatively minor.

## 2.2 Step tracing experiments

The experimental setup of solute tracing experiments is shown in Figure 3. The physical model is placed on a
horizontal platform to eliminate density effects. The physical model is first filled with deionized water. The inlet is connected to a solution tank (deionized water + NaCl, concentration $C_0 = 0.06$ mol/L). A peristaltic pump (Lead Fluid™) supplies the flow rate, $Q$ (L/s). At the outlet, a sensor (WTW TetraCon 325™ conductimeter, resolution $10^{-6}$ mol/L) is used to measure the tracer concentration and infer the BTC. The outflowing water falls into a glass placed on a weighing machine (Mettler Toledo™) that allows the flow rate to be monitored. The sensor and weighing machine are connected to a data logger (Campbell CR1000™) for automatic data sampling.

<figure><img src="image_4.png" /><figcaption>Figure 3. Tracing experiments setup.</figcaption></figure>

A step tracing test is conducted under a constant flow rate until the concentration at the outlet reaches $C_0$. To reduce measurement error, each experiment is replicated three times. The final data is obtained as the average of the three experimental BTCs, denoted as $C(t)$. Three different flow rates, spanning two orders of magnitudes, are used. For the three flow rates, the different data sampling intervals and total experimental time are shown in Table 2.

<table><caption>Table 2. Experimental flow rates.</caption><thead><tr><th>Pump speed (rpm)</th><th>Q (L/s)</th><th>Data sampling interval (s)</th><th>Experiment time (min)</th></tr></thead><tbody><tr><td>200</td><td>6.62 × 10<sup>-3</sup></td><td>1</td><td>10</td></tr><tr><td>20</td><td>6.84 × 10<sup>-4</sup></td><td>1</td><td>80</td></tr><tr><td>2</td><td>6.69 × 10<sup>-5</sup></td><td>5</td><td>600</td></tr></tbody></table>
# 3. Numerical modeling

## 3.1 Transport models

### a) Weighted Sum Advection-Dispersion Equation (WSADE) model with two regions

We consider applying the two-region WSADE model because the physical structures have two flow conduits.

The definition sketch of the two-region WSADE is shown in Figure 4, two regions are assumed to flow in parallel without any exchange; in both regions, the ADE model is assumed valid:

$$ \frac{\partial c_i}{\partial t} + u_i \frac{\partial c_i}{\partial x} - D_i \frac{\partial^2 c_i}{\partial x^2} = 0, (i = 1, 2) \qquad (1a) $$

$$ w_1 + w_2 = 1 \qquad (1b) $$

The total concentration at the outlet is obtained from a flux balance:

$$ C = \frac{w_1 u_1 c_1 + w_2 u_2 c_2}{w_1 u_1 + w_2 u_2} \qquad (2) $$

where $C_i$ is the outlet concentration (ML⁻³) of Region i, $D_i$ is the dispersion coefficient (L²T⁻¹), L = 1.5 m is the model length, t is time (T), $u_i$ is the flow velocity (LT⁻¹), x is the space coordinate (L), and $w_i$ is the volume fraction of Region i (-). Owing to Equation (1b), only one of the two fractions ($w_1, w_2$) is a free parameter. The parameter set ($u_1, u_2, D_1, D_2, w_1$) is calibrated to fit the experimental BTCs. To reduce the size of the parameter space, Regions 1 and 2 are distinguished by imposing $u_1 > u_2$ by definition.

<figure><img src="image_5.png" /><figcaption>Figure 4. Definition sketch of the WSADE model.</figcaption></figure>

### b) Dual Region Mobile ImMobile (DRMIM) model

The DRMIM model is based on the Mobile ImMobile (MIM) model (van Genuchten & Wierenga, 1976). The
MIM consists of a mobile region governed by the conventional ADE and an immobile region; the two regions are coupled by an exchange term. For the DRMIM model, the total concentration at the outlet is the superposition of two regions of MIM (Figure 5). The governing equations of the DRMIM model are:

$$ w_{m,i} \frac{\partial C_{m,i}}{\partial t} + w_{im,i} \frac{\partial C_{im,i}}{\partial t} = w_{m,i} D_{m,i} \frac{\partial^2 C_{m,i}}{\partial x^2} - w_{m,i} u_{m,i} \frac{\partial C_{m,i}}{\partial x}, \quad (i = 1, 2) \quad (3a) $$

$$ w_{im,i} \frac{\partial C_{im,i}}{\partial t} = k_i (C_{m,i} - C_{im,i}), \quad (i = 1, 2) \quad (3b) $$

$$ w_{im,i} + w_{m,i} = w_{MIM,i}, \quad (i = 1, 2) \quad (3c) $$

$$ w_{MIM,1} + w_{MIM,2} = 1 \quad (3d) $$

where $C_{m,1}$ and $C_{m,2}$ are the concentration (ML⁻³) of the two mobile regions, $C_{im,1}$ and $C_{im,2}$ are the concentration (ML⁻³) of the two immobile regions, $D_{m,1}$ and $D_{m,2}$ (L²T⁻¹) are the dispersion coefficients in the two MIM regions, $k_1$ and $k_2$ (T⁻¹) are the mass transfer coefficients; $u_{m,1}$ and $u_{m,2}$ (LT⁻¹) are the flow velocities; $w_{m,1}$ (-) and $w_{m,2}$ (-) are the volume fractions of the mobile regions; $w_{im,1}$ (-) and $w_{im,2}$ (-) are the volume fractions of the immobile regions; $w_{MIM,1}$ and $w_{MIM,2}$ are the volumetric fractions of the two MIM regions. The total concentration at the outlet is obtained from a flux balance over the two mobile regions:

$$ C = \frac{w_{m,1} u_{m,1} C_{m,1} + w_{m,2} u_{m,2} C_{m,2}}{w_{m,1} u_{m,1} + w_{m,2} u_{m,2}} \quad (4) $$

Owing to Equations (3c-d), the four fractions ($w_{m,1}$, $w_{m,2}$, $w_{im,1}$, $w_{im,2}$) have three degrees of freedom. In practice, nine independent parameters are calibrated to fit the experimental BTCs: $u_{m,1}$, $D_{m,1}$, $k_1$, $w_{m,1}$, $u_{m,2}$, $D_{m,2}$, $k_2$, $w_{m,2}$, $w_{MIM,1}$. As for the WSADE model, by definition $u_{m,1} > u_{m,2}$. The DRMIM model has more parameters and thus more flexibility than the WSADE model, then the DRMIM should be expected to better reproduce some BTCs than the WSADE.
<figure><img src="image_6.png" /><figcaption>Figure 5. Definition sketch of DRMIM model.</figcaption></figure>

## 3.2 Parameter calibration technique

The model parameters are calibrated with the Markov Chain Monte Carlo (MCMC) method (Haario et al., 2006). The MCMC method is a Bayesian approach that evaluates the posterior distributions of parameters (Vrugt et al., 2006). To achieve a fast convergence of the chain, we first perform a manual fitting of the model response to the experimental BTCs to find the initial parameter set. Then the Metropolis-Hastings algorithm (Robert, 2015) successively draws samples from the posterior distribution by forming a Markov chain of the model parameter set. In this work, we assume the posterior probability density is of a Gaussian type. For each chain, 10000 iterations are executed. The mean value of the last 50 sampled parameter sets is determined as the calibrated parameter value. The Root-Mean-Square-Error (RMSE) objective function is applied to evaluate the deviation and the convergence of the MCMC chains:

$$ \text{RMSE} = \sqrt{\frac{\sum_{i=1}^{n} (X_i - Y_i)^2}{n}} \qquad (5) $$

where X is the experimental data, n is the data quantity, Y is the simulated data.
205

# 4. Results

206

## 4.1 Experiment results

207

Some treatments are conducted to the concentration variations of the step tracing experiments, $C(t)$. First, the derivative of $C(t)$ to $t$, denoted as $C_{der}(t)$, is calculated and used for the exhibition of the experiment results and model benchmarking. If the injection signal is a concentration step, using $C_{der}(t)$ is equivalent to performing a tracing experiment using a Dirac (pulse) input signal (Majdalani et al., 2015), which is more approximate to the common practice in natural aquifers. Second, the $C_{der}(t)$ curves are normalized with pore volume:

208

209

210

211

212

$$PV(t) = \frac{tQ}{V_m} \qquad (6a)$$

213

$$C_N(t) = \frac{C_{der}(t)}{\int_{0}^{+\infty} C_{der}(t) \, dPV(t)} \qquad (6b)$$

214

215

216

217

where, $C_N$ is the normalized concentration (-), PV is the injected pore volume fraction (-), $V_m = 675$ cm³ is the volume of the physical model (L³). All of the experimental BTCs are exhibited in the form of $C_N$-PV. Because this transformation makes it convenient for us to compare the experiment BTCs of different flow rates. After the normalization, the area below the curve of $C_N$-PV is 1 (unitless).

218

### 4.1.1 Effect of flow rate

219

Figure 6 shows the experimental $C_N$-PV curves for the three flow rates. The 11 configurations are divided into two groups for plotting: for Group 1, $h_1 \le h_2$ (Figure 6a); for Group 2, $h_1 > h_2$ (Figure 6b).

220
<figure><img src="image_7.png" /><figcaption>Figure 6. Experimental C~N~-PV curves (Three flow rates).</figcaption></figure>

For the lowest flow rate of 6.69 × 10⁻⁵ L/s (2 rpm), the curves are right-skewed. And the curves of 2 rpm have the lowest peak value when compared to the other two flow rates. The curves of the middle flow rate (6.84 × 10⁻⁴ L/s, 20 rpm) have the largest peak value and sharpest peak. All of these curves of the middle flow rate are still single-peaked, while there is a bump on the falling limb of the 8-2 curve, which indicates the existence of two separate tracer plumes.

The curves of the highest flow rate (6.62×10⁻³ L/s, 200 rpm) have a medium peak value. Under this flow rate, there are three types of curves: 1) some curves show obvious dual peaks (5-5, 4-6, 9-1); 2) other curves have one peak while there is an obvious bump on the falling limb (6-4, 7-3, 8-2), these BTCs are similar to the dual-peaked
231

BTCs, but the two peaks are so close that they superpose; 3) some configurations exhibit curves with only one peak (3-7, 2-8, 1-9, 0-10, 10-0). When compared to lowest and middle flow rates, more curves of the highest flow rate exhibit dual-peaked curves.

232

233

234

### 4.1.2 Effect of aperture

235

236

237

The effect of aperture variation is analyzed using the curves of the highest flow rate because this is the flow rate for which the contrast between the BTCs is the most pronounced. The 11 BTCs are divided into two groups for plotting: for Group 1, $h_1 \le h_2$ (Figure 7a); for Group 2, $h_1 \ge h_2$ (Figure 7b).

<figure><img src="image_8.png" /><figcaption>Figure 7. Experiment CN-PV curves of highest flow rate.</figcaption></figure>

238

239

240

241

242

243

244

245

In Figure 7a, Configurations 5-5 and 4-6 yield BTCs with two peaks, which indicates a significant difference between the travel times in the main and diverging conduit. For 5-5, the early peak is larger than the later peak. The BTC for Configuration 4-6 is rather unusual in that the early peak is slightly smaller than the late peak, while dual-peaked BTCs are known to mostly exhibit a larger early peak than the late peak (Field & Leij, 2012; Florea & Wicks, 2001; Mohammadi et al., 2019; Wang et al., 2020). The remaining configurations (3-7, 2-8, 1-9, 0-10) exhibit single-peaked BTCs. Those BTCs are similar.

246

247

248

In Figure 7b, the BTCs of Configurations 6-4, 7-3 and 8-2 exhibit a single peak, but a bump is visible within the falling limb. These BTCs appear as intermediate between the dual-peaked and single-peaked BTCs. As the configuration changes from 5-5 to 8-2, the BTCs peak (the early plume) occurs increasingly earlier; as the
configuration changes from 8-2 to 10-0, the BTCs peak occurs increasingly later. As the configuration changes from 6-4 to 10-0, the late plume (bump or the second peak) occurs increasingly delayed.

Among all of the configurations, the BTCs obtained by 5-5 and 4-6 exhibit lower peak values than the other configurations. The falling limbs of these two BTCs are also more delayed than the others. Hauns et al., (2001) showed that different karst geometry may cause different streamline distributions, thus different degrees of solute detention effects. The solute transport processes in the experiments should also be different due to the different structure of the configurations. There may be a stronger detention effect for the transport processes in Configurations 5-5 and 4-6, which has caused these two BTCs to be more delayed.

Comparing Group 1 to Group 2, we find that Configurations 4-6 and 6-4 exhibit very different BTCs: for 4-6, the early peak is smaller than the late peak; for 6-4, there is a bump on the falling limb. Similarly, Configuration 3-7 exhibits different results from 7-3, etc. The difference between the two groups suggests that both conduit lengths and apertures influence the transport process.

In this study, the double peaks overlap to a large extent (5-5, 4-6). In contrast with dual-peaked BTCs reported by (Field & Leij, 2012; Wang et al., 2020), the present experiments exhibit only a slight separation between the two peaks. This can be explained by a small contrast between the lengths of the conduits (length ratio = 125/105 ≈ 1.19) for our dual-conduit structures. As discussed in (Wang et al., 2020), the closer the length ratio to unity, the closer the two peaks are.

## 4.2 Modeling results

In this section, we present the modeling results to the BTCs of the highest flow rate (6.62×10⁻³ L/s, 200 rpm), and investigate how the aperture setup may influence the model parameters. The modeling results of the 2-rpm curves and the 20-rpm curves can be found in Appendix A.
<figure><img src="image_9.png" /><figcaption>Figure 8. WSADE model fitting 200-rpm experimental BTCs.</figcaption></figure>

271

272

273

274

275

276

277

278

279

280

281

The WSADE model can reproduce most of the BTCs (Figure 8). The RMSE values range between $2.59 \times 10^{-4}$ and $6.36 \times 10^{-4}$ m²/s (Table 4). For Configuration 5-5, the dual-peaked BTC is well reproduced (RMSE = $4.87 \times 10^{-4}$): Region 1 accounts for the first peak of the BTC, and Region 2 accounts for the second one. For Configuration 4-6, WSADE exhibits a degraded fitting performance (RMSE = $5.69 \times 10^{-4}$), and the WSADE model fails to reproduce the double peak behavior. For Configurations 3-7, 2-8, 1-9, 0-10, and 10-0, the WSADE well reproduces the single-peaked BTCs (RMSE ranges between $2.59 \times 10^{-4}$ and $3.96 \times 10^{-4}$). For Configurations 6-4, 7-3, and 8-2, the fitting performances are also good (RMSE ranges between $2.81 \times 10^{-4}$ and $6.36 \times 10^{-4}$): Region 1 accounts for the first peak, and Region 2 accounts for the bump on the falling limb. For Configuration 9-1, the fitting performance is degraded, as the second peak is not well reproduced (RMSE = $4.70 \times 10^{-4}$).
<table><caption>Table 3. Calibrated parameters of WSADE to 200-rpm experiments</caption><thead><tr><th>Pump Speed (rpm)</th><th>h<sub>1</sub>-h<sub>2</sub> (cm)</th><th>M (cm)</th><th>u<sub>1</sub> (m/s)</th><th>u<sub>2</sub> (m/s)</th><th>D<sub>1</sub> (m<sup>2</sup>/s)</th><th>D<sub>2</sub> (m<sup>2</sup>/s)</th><th>w<sub>1</sub> (-)</th><th>RMSE</th></tr></thead><tbody><tr><td rowspan="10">200</td><td>5-5</td><td>0</td><td>1.77×10<sup>-2</sup></td><td>1.38×10<sup>-2</sup></td><td>6.34×10<sup>-5</sup></td><td>5.47×10<sup>-4</sup></td><td>0.099</td><td>4.87×10<sup>-4</sup></td></tr><tr><td>4-6</td><td>1</td><td>1.47×10<sup>-2</sup></td><td>9.71×10<sup>-3</sup></td><td>4.13×10<sup>-4</sup></td><td>6.16×10<sup>-4</sup></td><td>0.793</td><td>5.69×10<sup>-4</sup></td></tr><tr><td>3-7</td><td>2</td><td>1.61×10<sup>-2</sup></td><td>1.28×10<sup>-2</sup></td><td>2.00×10<sup>-4</sup></td><td>7.04×10<sup>-4</sup></td><td>0.609</td><td>2.59×10<sup>-4</sup></td></tr><tr><td>2-8</td><td>3</td><td>1.64×10<sup>-2</sup></td><td>1.34×10<sup>-2</sup></td><td>1.45×10<sup>-4</sup></td><td>6.34×10<sup>-4</sup></td><td>0.515</td><td>3.96×10<sup>-4</sup></td></tr><tr><td>1-9</td><td>4</td><td>1.69×10<sup>-2</sup></td><td>1.36×10<sup>-2</sup></td><td>2.00×10<sup>-4</sup></td><td>8.94×10<sup>-4</sup></td><td>0.556</td><td>2.97×10<sup>-4</sup></td></tr><tr><td>0-10</td><td>5</td><td>1.64×10<sup>-2</sup></td><td>1.32×10<sup>-2</sup></td><td>1.71×10<sup>-4</sup></td><td>7.47×10<sup>-4</sup></td><td>0.527</td><td>3.37×10<sup>-4</sup></td></tr><tr><td>6-4</td><td>-1</td><td>1.83×10<sup>-2</sup></td><td>1.54×10<sup>-2</sup></td><td>9.48×10<sup>-5</sup></td><td>5.48×10<sup>-4</sup></td><td>0.144</td><td>6.36×10<sup>-4</sup></td></tr><tr><td>7-3</td><td>-2</td><td>1.85×10<sup>-2</sup></td><td>1.49×10<sup>-2</sup></td><td>1.38×10<sup>-4</sup></td><td>7.92×10<sup>-4</sup></td><td>0.294</td><td>4.18×10<sup>-4</sup></td></tr><tr><td>8-2</td><td>-3</td><td>1.85×10<sup>-2</sup></td><td>1.39×10<sup>-2</sup></td><td>1.46×10<sup>-4</sup></td><td>8.32×10<sup>-4</sup></td><td>0.385</td><td>2.81×10<sup>-4</sup></td></tr><tr><td>9-1</td><td>-4</td><td>1.74×10<sup>-2</sup></td><td>1.17×10<sup>-2</sup></td><td>2.04×10<sup>-4</sup></td><td>1.22×10<sup>-3</sup></td><td>0.597</td><td>4.70×10<sup>-4</sup></td></tr><tr><td>10-0</td><td>-5</td><td>1.70×10<sup>-2</sup></td><td>1.34×10<sup>-2</sup></td><td>1.69×10<sup>-4</sup></td><td>6.56×10<sup>-4</sup></td><td>0.590</td><td>3.72×10<sup>-4</sup></td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>Figure 9. Calibrated parameters of WSADE for 200-rpm experiments.</figcaption></figure>

We examine the calibrated WSADE parameters for the configurations with $M \ge 1$ cm, the velocity $u_1$ is larger than $u_2$, the dispersion $D_1$ is smaller than $D_2$, and the volume fraction $w_1$ is larger than 0.5. Configuration 4-6 ($M = 1$ cm) presents the smallest $u_1$ and the smallest $u_2$ among all configurations. The small values of $u_1$ and $u_2$ indicate that $t_{peak1}$ and $t_{peak2}$ (the peak time of WSADE Region 1 and Region 2, respectively) are larger than the other configurations. Configuration 4-6 presents the largest $D_1$ (4.13×10⁻⁴ m²/s) and the largest $w_1$ (0.793). This is

283

284

285

286

287

288

289
because the WSADE Region 1 covers most of the BTC while Region 2 covers only a small part of the tail for 4-6 (Figure 8). From Configurations 3-7 to 0-10 ($M = 2$ to $5$ cm), $u_1$ ranges between $1.61 \times 10^{-2}$ and $1.69 \times 10^{-2}$ m/s, $u_2$ ranges between $1.28 \times 10^{-2}$ and $1.36 \times 10^{-2}$ m/s; $D_1$ ranges between $1.45 \times 10^{-4}$ and $2.00 \times 10^{-4}$, $D_2$ ranges between $6.34 \times 10^{-4}$ to $8.94 \times 10^{-4}$ m$^2$/s. $w_1$ shows decreasing trend on the whole: it decreases from 0.609 to 0.527.

Then, we introduce the calibrated WSADE parameters for the configurations with negative $M$. From Configuration 5-5 to 8-2 ($M=0, -1, -2, -3$), $u_1$ increases from $1.77 \times 10^{-2}$ to $1.85 \times 10^{-2}$; from 8-2 to 10-0 ($M=-3, -4, -5$) $u_1$ decreases to $1.70 \times 10^{-2}$. From 5-5 to 6-4, $u_2$ increases from $1.38 \times 10^{-2}$ to $1.54 \times 10^{-2}$; then from 6-4 to 9-1, $u_2$ decreases from $1.54 \times 10^{-2}$ to $1.17 \times 10^{-2}$. The values $u_1$ and $u_2$ reflect the variations in the arrival times of the peaks in Figure 7b. From 5-5 to 9-1 ($M=0 \sim -4$), $D_1$ increases from $6.34 \times 10^{-5}$ to $2.04 \times 10^{-4}$ and $D_2$ increases from $5.47 \times 10^{-4}$ to $1.22 \times 10^{-3}$; $w_1$ (the volume fraction of the faster region) increases from 0.099 to 0.597.

$D_2$ is larger than $D_1$ for all flow configurations. Larger $D_2$ values are needed so that WSADE Region 2 can generate a wider peak to reproduce the later tail of the experimental BTCs. Further explanations are provided in Section 5.3.
<figure><img src="image_11.png" /><figcaption>Figure 10. DRMIM model fitting 200-rpm experimental BTCs.</figcaption></figure>

The DRMIM can well reproduce all BTCs (Figure 10), capturing all of the features of the BTCs: the peaks, the bumps, and the tails. Visually, the DRMIM exhibits better characterization for the dual-peaked curves (5-5, 4-6, 9-1) and the bump-tail curves (6-4, 7-3, 8-2) than the WSADE (Figure 8). We calculate the RMSE values of the two models fitting results in Figure 8 and Figure 10. Quantitatively, the RMSE values of the DRMIM are much lower than the WSADE (Figure 11).
<figure><img src="image_12.png" /><figcaption>Figure 11. Comparison of the fitting performance of the two models to the 200-rpm curves.</figcaption></figure>

The calibrated DRMIM parameters are presented in Table 4 and shown in Figure 12. As shown in Figure 12a, $u_{m,1}$ value has exhibited a relatively smaller variance than $u_{m,2}$. For the single-peaked BTCs ($M=2, 3, 4, 5, -5$ cm), the contrasts between $u_{m,1}$ and $u_{m,2}$ are relatively smaller; for the dual-peaked BTCs ($M=0, 1, -4$ cm) or the bump-tailed BTCs ($M=-1, -2, -3$ cm), the contrasts between $u_{m,1}$ and $u_{m,2}$ are relatively larger. As the Configuration changes from 5-5 to 8-2 ($M$ decreases from 0 to -3), the early peak arrives increasingly earlier, thus $u_{m,1}$ increases from $1.74 \times 10^{-2}$ to $1.83 \times 10^{-2}$ m/s; from 8-2 to 10-0 ($M$ decreases from -3 to -5), the early peak arrives increasingly delayed, thus $u_{m,1}$ decreases from $1.83 \times 10^{-2}$ to $1.70 \times 10^{-2}$ m/s. As the Configuration changes from 6-4 to 9-1 ($M$ decreases from -1 to -4), the late peak arrives increasingly delayed, thus $u_{m,2}$ decreases from $1.46 \times 10^{-2}$ to $8.27 \times 10^{-3}$ m/s.

While parameters $D_{m,1}$, $D_{m,2}$, $k_1$, $k_2$, $w_{im,1}$, $w_{im,2}$ do not show continuous variations (Figure 12b-e). The variations of these parameters may be caused by parameter equifinality, which means that more than one parameter set will reproduce the observations equally well concerning a particular performance measure (Cobelli & Distefano, 1980). Equifinality can be caused by parameters that are insensitive (Johnston & Pilgrim, 1976) and/or interactive (Ibbitt & O'Donnell, 1974). In Appendix B, we show that the problem of parameter equifinality does exist in the parameter calibration process. As $M$ increases from -4 to 5, $w_{m,1}$ decreases from 0.675 to 0.370. $w_{m,2}$ increases from 0.206 to 0.423.

<table><caption>Table 4. Calibrated parameters of DRMIM to 200-rpm experiments</caption><thead><tr><th></th><th>h<sub>1</sub>-h<sub>2</sub></th><th>M (cm)</th><th>u<sub>m,1</sub> (m/s)</th><th>D<sub>m,1</sub> (m<sup>2</sup>/s)</th><th>k<sub>1</sub> (1/s)</th><th>w<sub>m,1</sub> (-)</th><th>u<sub>m,2</sub> (m/s)</th><th>D<sub>m,2</sub> (m<sup>2</sup>/s)</th><th>k<sub>2</sub> (1/s)</th><th>w<sub>m,2</sub> (-)</th><th>w<sub>MIM,1</sub> (-)</th><th>RMSE</th></tr></thead><tbody><tr><td>200</td><td>5-5</td><td>0</td><td>1.74×10<sup>-2</sup></td><td>1.24×10<sup>-4</sup></td><td>9.88×10<sup>-3</sup></td><td>0.491</td><td>1.27×10<sup>-2</sup></td><td>1.11×10<sup>-4</sup></td><td>2.26×10<sup>-3</sup></td><td>0.285</td><td>0.619</td><td>1.49×10<sup>-4</sup></td></tr></tbody></table>
<table><thead><tr><th>rpm</th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>4-6</td><td>1</td><td>1.69×10<sup>-2</sup></td><td>1.37×10<sup>-4</sup></td><td>1.09×10<sup>-2</sup></td><td>0.477</td><td>1.32×10<sup>-2</sup></td><td>7.06×10<sup>-5</sup></td><td>4.91×10<sup>-3</sup></td><td>0.357</td><td>0.561</td><td>1.31×10<sup>-4</sup></td></tr><tr><td>3-7</td><td>2</td><td>1.67×10<sup>-2</sup></td><td>1.74×10<sup>-4</sup></td><td>9.35×10<sup>-3</sup></td><td>0.530</td><td>1.57×10<sup>-2</sup></td><td>1.39×10<sup>-4</sup></td><td>4.25×10<sup>-3</sup></td><td>0.284</td><td>0.597</td><td>6.21×10<sup>-5</sup></td></tr><tr><td>2-8</td><td>3</td><td>1.65×10<sup>-2</sup></td><td>1.02×10<sup>-4</sup></td><td>4.46×10<sup>-3</sup></td><td>0.449</td><td>1.62×10<sup>-2</sup></td><td>2.67×10<sup>-4</sup></td><td>1.36×10<sup>-2</sup></td><td>0.377</td><td>0.589</td><td>7.70×10<sup>-5</sup></td></tr><tr><td>1-9</td><td>4</td><td>1.69×10<sup>-2</sup></td><td>1.60×10<sup>-4</sup></td><td>2.28×10<sup>-3</sup></td><td>0.434</td><td>1.50×10<sup>-2</sup></td><td>4.72×10<sup>-4</sup></td><td>1.54×10<sup>-3</sup></td><td>0.424</td><td>0.542</td><td>8.17×10<sup>-5</sup></td></tr><tr><td>0-10</td><td>5</td><td>1.72×10<sup>-2</sup></td><td>1.80×10<sup>-4</sup></td><td>4.61×10<sup>-3</sup></td><td>0.370</td><td>1.64×10<sup>-2</sup></td><td>8.22×10<sup>-5</sup></td><td>1.22×10<sup>-2</sup></td><td>0.423</td><td>0.519</td><td>3.95×10<sup>-5</sup></td></tr><tr><td>6-4</td><td>-1</td><td>1.78×10<sup>-2</sup></td><td>1.59×10<sup>-4</sup></td><td>4.39×10<sup>-3</sup></td><td>0.557</td><td>1.46×10<sup>-2</sup></td><td>6.79×10<sup>-5</sup></td><td>7.92×10<sup>-3</sup></td><td>0.268</td><td>0.712</td><td>1.24×10<sup>-4</sup></td></tr><tr><td>7-3</td><td>-2</td><td>1.78×10<sup>-2</sup></td><td>2.18×10<sup>-4</sup></td><td>2.78×10<sup>-3</sup></td><td>0.610</td><td>1.40×10<sup>-2</sup></td><td>3.22×10<sup>-5</sup></td><td>1.55×10<sup>-2</sup></td><td>0.199</td><td>0.775</td><td>1.30×10<sup>-4</sup></td></tr><tr><td>8-2</td><td>-3</td><td>1.83×10<sup>-2</sup></td><td>1.59×10<sup>-4</sup></td><td>5.66×10<sup>-3</sup></td><td>0.624</td><td>1.22×10<sup>-2</sup></td><td>1.25×10<sup>-4</sup></td><td>4.64×10<sup>-3</sup></td><td>0.212</td><td>0.714</td><td>9.55×10<sup>-5</sup></td></tr><tr><td>9-1</td><td>-4</td><td>1.77×10<sup>-2</sup></td><td>1.47×10<sup>-4</sup></td><td>7.36×10<sup>-3</sup></td><td>0.675</td><td>8.27×10<sup>-3</sup></td><td>1.19×10<sup>-4</sup></td><td>4.84×10<sup>-3</sup></td><td>0.206</td><td>0.750</td><td>1.52×10<sup>-4</sup></td></tr><tr><td>10-0</td><td>-5</td><td>1.70×10<sup>-2</sup></td><td>1.30×10<sup>-4</sup></td><td>3.46×10<sup>-3</sup></td><td>0.480</td><td>1.51×10<sup>-2</sup></td><td>3.63×10<sup>-4</sup></td><td>2.70×10<sup>-3</sup></td><td>0.337</td><td>0.522</td><td>3.31×10<sup>-5</sup></td></tr></tbody></table>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>Figure 12. Calibrated parameters of DRMIM to 200-rpm experiments.</figcaption></figure>
# 5. Discussion

## 5.1 Representativeness of the experiments for karst tracer tests

The purpose of the present subsection is to investigate how representative of field scale configurations the present experiments are. The ratio of conduit length to diameter ratio, the Péclet number, and the Reynolds number is computed from the model calibrated parameters and given in Table 5.

<table><caption>**Table 5. Calibrated and field-scale dimensionless parameters.**</caption><thead><tr><th>Parameter</th><th>WSADE calibration range</th><th>Field-scale range</th></tr></thead><tbody><tr><td>L / d<sub>c</sub></td><td>[158, 225]</td><td>[25, 25000]</td></tr><tr><td>Pé</td><td>Pe<sub>1</sub> ∈ [53, 418]<br>Pe<sub>2</sub> ∈ [14, 83]</td><td>[12, 331]</td></tr><tr><td>Re</td><td>Re<sub>1</sub> ∈ [1.5, 168]<br>Re<sub>2</sub> ∈ [1.1, 140]</td><td>[6.5×10<sup>3</sup>, 8×10<sup>4</sup>]</td></tr></tbody></table>

The characteristic diameter of the conduit in rectangle shape is calculated as:

$$d_c = \frac{2hw}{(h+w)} \qquad (7)$$

where *h* is the conduit aperture, *w* = 0.005 m is the thickness of the flow channel. Since we know the model length is *L* = 1.5 m. The calculated values of *L* / *d*<sub>c</sub> range from 158 to 225. The values of the length to diameter ratio in natural karst aquifers are found to range from 25 to 25 000 (according to the summary in Wang et al., 2020). Thus the experimental values are reasonable.

The Péclet number is a measure of the relative importance of advection versus diffusion. For a solute transport process, the longitudinal dispersion (*D*<sub>L</sub>, L<sup>2</sup>T<sup>-1</sup>) can be used instead of diffusion in the calculation of the Péclet number:

$$Pé = \frac{Lu}{D_L} \qquad (8)$$

We apply the calibrated parameters of the WSADE for the calculation of this number. The calibrated parameters of *u*<sub>1</sub> and *u*<sub>2</sub> are used for this *u*; values of *D*<sub>1</sub> and *D*<sub>2</sub> are used for this *D*<sub>L</sub>. With *u*<sub>1</sub> and *D*<sub>1</sub>, the calculated Pé values for the experiments range from 53.5 to 418.4; with *u*<sub>2</sub> and *D*<sub>2</sub>, the calculated Pé values range from 14.4 to 82.8. *D*<sub>1</sub> is
always much smaller than $D_2$ (Figure 9b), so the Pé calculated by ($u_1$ and $D_1$) are larger than that by ($u_2$ and $D_2$).
According to the tracing tests carried out in real karst systems, the calculated Pé values range from 12–331 (summarized in Wang et al., 2020). Most of the experimental Pé values are within the range in the field tests.

We check the flow regime of the experiments by calculating the Reynolds number:

$$Re = \frac{\rho u D_H}{\mu} = 4 \frac{\rho u A}{\mu P} = 4 \frac{\rho Q}{\mu P} \qquad (9)$$

where $A$ (L$^2$) is the cross-section area, $D_H$ (L) is the hydraulic diameter of the pipe, $P$ (L) is the wetted perimeter, $Q$ (L$^3$T$^{-1}$) is the flow rate, $u$ (LT$^{-1}$) is the average flow velocity, $\mu$ (ML$^{-1}$T$^{-1}$) is the fluid dynamic viscosity and $\rho$ (ML$^{-3}$) is the fluid density. Again the calibrated WSADE parameters of $u_1$ and $u_2$ are used for this calculation. According to Table 5, the Reynolds number of the experiments ranges from 1 to 168. While researchers have noted that Re varies between 6500 and 80000 in natural karst aquifers (Duran et al., 2016; Morales et al., 2010).

In summary, our experimental models can be representative of the natural karst systems in terms of geometric and kinematic similarities ($L$ / $d_c$ and Pé), while not in the dynamic similarity (Re). The flow regime of the experiments is laminar (Re < 2000). We are clear that the turbulent flow regime is more generally observed in natural karst conduits. In practice, it would be difficult to achieve a turbulent flow regime owing to the experimental dimensions and the availability of deionized water. Note however that lab-scale tracing experiments under a laminar flow regime have been shown to provide insights into solute transport processes in natural karst systems (Wang et al., 2020).

## 5.2 Effect of flow rate.

According to the experimental results, the flow rate influences the BTCs significantly (Figure 6). The dual-peak effect is more obvious under a higher flow rate. This result is in agreement with the experimental study by Mohammadi et al. (2019): under high hydraulic gradients (thus large flow rate), the tracing tests in the conduit network of branchwork structure generated multi-peaked BTCs; under low hydraulic gradients (thus small flow rates), the network exhibited single-peaked BTCs.
374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

For a given conduit geometry multi-peaked BTCs may bring us more information about the system investigated by the tracing test. Firstly, dual-peaked BTC suggests the existence of flow bifurcations. Secondly, they allow an average transport velocity to be identified for each of the two regions. Thirdly, we can estimate the underground conduit length of different channels as proposed by Wang et al. (2020).

## 5.3 Application of experiment results

Comparing the results of the present study to those of Wang et al. (2020) allows qualitative deductions to be derived about the geometry of the underground karst structure from dual-peaked BTCs from tracing tests in karst aquifers.

First, if the two peaks are fully separated from each other, the lengths of the two conduits are likely to be significantly different. According to Figure 7, a length ratio close to unity leads to peak superposition. Thus we can deduce a large length contrast for the two conduits when the two peaks are fully separated. For this kind of dual-peaked BTC, a method has been proposed to estimate the lengths of the conduits (Wang et al., 2020).

Second, if the two peaks overlap each other and the early peak is smaller than the late peak, the length ratio may be very small (close to unity) and the shorter conduit has a smaller aperture than the longer conduit ($h_1 < h_2$). This BTC is exhibited by Configuration 4-6 (Figure 7a). To our best knowledge, such behaviors have been very rarely reported (if at all) in the literature.

Third, if the two peaks overlap each other and the early peak is larger than the late peak. There are two possible cases. The first case is that the length ratio is small and the shorter conduit has a larger aperture than the longer conduit ($h_1 \ge h_2$), the BTCs in Figure 7b belong to this type. The second case is that the conduit is connected to a waterfall or large pool (Field & Leij, 2012).

The above deductions are summarized in Table 6. The deductions are qualitative rather than quantitative because karst structures can be highly complex. If researchers want to obtain more information about the underground karst geometry, further tests may be required.
<table><caption>Table 6. Interpretation of dual-peaked BTCs to deduce geometry of underground dual-conduit structures.</caption><thead><tr><th>BTC shape</th><th>Example</th><th>Deductions</th></tr></thead><tbody><tr><td>Two peaked are separated.</td><td>Length ratio, 6 : 1<br>Experiment data<br>-Modeling</td><td>Large contrast between the conduit lengths.</td></tr><tr><td>Two peaked are not separated.<br>Early peak is smaller.</td><td>Configuration 4-6, highest flow rate.<br>C<sub>N</sub></td><td>Length ratio is close to unity, and h<sub>1</sub> &lt; h<sub>2</sub></td></tr><tr><td>Two peaked are not separated.<br>Early peak is larger.</td><td>Configuration 9-1, highest flow rate.<br>C<sub>N</sub></td><td>Length ratio is close to unity and h<sub>1</sub> &gt; h<sub>2</sub>.<br>Or a large pool is connected to the conduit.</td></tr></tbody></table>

## 5.4 Fitting performance of the numerical models

In our application of the WSADE model to fitting the experimental BTCs. We find some WSADE model parameters do not bear the physical meaning we have assumed. For example, when we apply the WSADE to fit the BTC of 5-5, 200 rpm (Figure 8, Table 4), we find that the calibrated values of $w_1$ and $w_2$ cannot correctly represent the true tracer quantity repartition in the two conduits. Because the conceptual structure of the WSADE is similar to the experimental model, in that the WSADE has two regions and the experimental model has two conduits. We have expected that the WSADE Region 1 parameters ($u_1$, $D_1$, and $w_1$) should represent the properties of the transport process of the shorter conduit, which makes the early peak of the BTC; and the WSADE Region 2 parameters ($u_2$, $D_2$, and $w_2$) should represent the properties of the transport process of the longer conduit, which makes the late peak. The calibrated value $w_1 = 0.099$ and $w_2 = 1 - w_1 = 0.901$. $w_1$ is nine times smaller than $w_2$. And according to the experimental BTC, the early peak is higher than the late peak. Despite this peak magnitude is not sufficient to indicate the tracer plume quantity, the comparison of the two peaks still indicates that the quantity of
410

411

the early tracer plume cannot be nine times smaller than the late tracer plume. Thus the values of $w_1$ and $w_2$ are
unrealistic and cannot correctly represent the true quantity repartition between the two tracer plumes.

412

413

414

415

416

417

418

In this fitting, the values of $D_1$ and $D_2$ are also not realistic: the calibrated value of $D_2$ is around 8.6 times
larger than $D_1$ ($D_1$ = 6.34×10⁻⁵ m²/s, $D_2$=5.47×10⁻⁴ m²/s). The contrast between the two values is too large: the
length ratio of the two conduits is around 1.13 and the cross-section areas of the two conduits are the same for
Configuration 5-5. The properties of the two conduits cannot differ by 8.6 times. Thus, parameters $D_1$ and $D_2$
cannot correctly represent the dispersion effect of the two conduits. For the WSADE model fitting the experiment
of 5-5, 200 rpm, some parameters ($w_1, w_2, D_1$, and $D_2$) of WSADE do not represent the physical meanings we have
expected.

419

420

421

422

423

424

425

426

According to the characteristic of the experimental BTCs and the WSADE model, we explain the cause of the
above problem. First, according to Figure 6, the BTCs under 20-rpm and 200-rpm exhibit a significant tailing effect
for the falling limb. The tailing effect has caused the falling limb to last for a relatively long time scale. The tailing
effect has been caused by the solute detention effect introduced in Section 4.1.3. Second, the two-region WSADE
consists of two parallel-flowing regions governed by the ADE model (Figure 4). It is well acknowledged that the
simple ADE model cannot well reproduce the BTCs with the tailing effect. To properly reproduce the falling limb,
the WSADE Region 2 is more widespread and has a larger size. So, the calibrated values of $D_2$ and $w_2$ are
unrealistically large.

427

428

429

430

431

432

Compared to the two-region WSADE, the DRMIM catches better the transport process in the dual-conduit
structures (Figure 11). The DRMIM has a better fitting performance because the conceptual structure of the
DRMIM is more similar to the experimental solute transport process in the dual-conduit structures. Experimental
studies have proved that solute transport in a single fracture (Qian et al., 2011), or a single conduit (Anger &
Alexander, 2013; Field & Leij, 2012; Zhao et al., 2019, 2021) can exhibit some anomalous features like right-
skewness and tailing. The dual-conduit structures have two channels, thus one experiment can be deemed as the
combination of two anomalous transport processes. As the ADE cannot well characterize anomalous BTCs. Similarly, the WSADE (which includes two ADE regions) cannot well characterize the experimental BTCs in this study. However, it is known is the MIM can well characterize anomalous BTCs with right-skewness and tailing, thus the DRMIM (which includes two MIM regions) can well characterize the experimental BTCs from the dual-conduit structures.

## 6. Conclusions

We studied the solute transport process in dual conduit structures through lab-scale experiments. By definition, $h_1$ is the aperture of the shorter conduit and $h_2$ is the aperture of the longer conduit. If $h_1 < h_2$ and the contrast between the two apertures is low, the dual-conduit structures generate dual-peaked BTCs that exhibit a lower early peak than the late peak; if $h_1 < h_2$ and the contrast between $h_1$ and $h_2$ is large, we obtain a single-peaked BTCs; if $h_1 \ge h_2$, dual-peaked BTCs or bump-tailed BTCs are observed. As the flow rate increases, it is more likely to observe the separation of the two tracer plumes.

In a previous experimental study (Wang et al., 2020), we investigated how three factors may influence the transport process in dual-conduit structures (length ratio, total length, and connection angle). The current study is a complementary work to this previous study, revealing the effects of the aperture and flow rate. The experimental results should help the hydrologists to better interpret the dual-peaked BTCs from karst aquifers. We provide a guide (Table 6) on estimating some properties related to the geometry of the dual-conduit structure according to the shape of the dual-peaked BTCs.

The experimental BTCs are fitted by two numerical models. The WSADE model properly reproduces the bump-tailed BTCs and some of the dual-peaked BTCs, but it failed to correctly reproduce the dual-peaked BTC obtained from configuration 4-6 under 200 rpm. For the first time, the DRMIM model is applied to characterize the anomalous transport BTCs. The DRMIM exhibited better performances than the WSADE model and captured all of the symptoms of the BTCs: the second peak, the bump, and the tailing.
456

Acknowledgments

457

Chaoqi Wang is supported by the China Scholarship Council (CSC) from the Ministry of Education of P.R. China.

458

H. J. is grateful for financial support from the European Commission through the Partnership for Research and

459

Innovation in the Mediterranean Area (PRIMA) program under Horizon 2020 (KARMA project, Grant No.

460

01DH19022A).
461 References

462 Anaya, A. A., Padilla, I., Macchiavelli, R., Vesper, D. J., Meeker, J. D., & Alshawabkeh, A. N. (2014).

463 Estimating preferential flow in karstic aquifers using statistical mixed models. Groundwater, 52(4), 584–596.

464 Anger, C. T., & Alexander, E. C. (2013). Bench-scale models of dye breakthrough curves. Carbonates and
465 Evaporites, 28(1), 221–227.

466 Bakalowicz, M. (2005). Karst groundwater: a challenge for new resources. Hydrogeology Journal, 13(1), 148–160.

468 Barnard, A. C. L., Hunt, W. A., Timlake, W. P., & Varley, E. (1966). A theory of fluid flow in compliant tubes.
469 Biophysical Journal, 6(6), 717–724.

470 Bear, J. (1988). Dynamics of fluids in porous media. Courier Corporation.

471 Berkowitz, B., Cortis, A., Dentz, M., & Scher, H. (2006). Modeling non-Fickian transport in geological formations as a continuous time random walk. Reviews of Geophysics, 44(2).

473 Brauchler, R., Böhm, G., Leven, C., Dietrich, P., & Sauter, M. (2013). A laboratory study of tracer tomography. Hydrogeology Journal, 21(6), 1265–1274.

475 Brauchler, R., Liedl, R., & Dietrich, P. (2003). A travel time based hydraulic tomographic approach. Water Resources Research, 39(12).

477 Castro, R. B. P. (2017). Statistical Analysis of Karst Aquifer Pollution, Karst Flow Model Validation at Laboratory Scale, and Development of Seepage Meter. The Florida State University.

479 Cherubini, C., Giasi, C. I., & Pastore, N. (2012). Bench scale laboratory tests to analyze non-linear flow in fractured media. Hydrology and Earth System Sciences, 16(8), 2511–2522.

481 Cherubini, C., Pastore, N., Giasi, C. I., & Allegretti, N. M. (2017). Laboratory experimental investigation of heat transport in fractured media. Nonlinear Processes in Geophysics, 24(1), 23–42.

483 Cobelli, C., & Distefano, J. J. (1980). Parameter and Structural Identifiability Concepts and Ambiguities – a
484

Critical-Review and Analysis. The American Journal of Physiology, 239(1), R7-24.

485

Develi, K., & Babadagli, T. (2015). Experimental and visual analysis of single-phase flow through rough fracture replicas. *International Journal of Rock Mechanics and Mining Sciences*, 73, 139–155.

486

487

Dewaide, L., Bonniver, I., Rochez, G., & Hallet, V. (2016). Solute transport in heterogeneous karst systems: Dimensioning and estimation of the transport parameters via multi-sampling tracer-tests modelling using the OTIS (One-dimensional Transport with Inflow and Storage) program. *Journal of Hydrology*, 534, 567–578.

488

489

490

DiFrenna, V. J., Price, R. M., & Savabi, M. R. (2008). Identification of a hydrodynamic threshold in karst rocks from the Biscayne Aquifer, south Florida, USA. *Hydrogeology Journal*, 16(1), 31–42.

491

492

Duran, L., Fournier, M., Massei, N., & Dupont, J. (2016). Assessing the Nonlinearity of Karst Response Function under Variable Boundary Conditions. *Ground Water*, 54.

493

494

Faulkner, J., Hu, B. X., Kish, S., & Hua, F. (2009). Laboratory analog and numerical study of groundwater flow and solute transport in a karst aquifer with conduit and matrix domains. *Journal of Contaminant Hydrology*, 110(1–2), 34–44.

495

496

497

Field, M. S., & Leij, F. J. (2012). Solute transport in solution conduits exhibiting multi-peaked breakthrough curves. *Journal of Hydrology*, s 440–441(none), 26–35.

498

499

Florea, L. J., & Wicks, C. M. (2001). Solute transport through laboratory-scale karstic aquifers. *Journal of Cave and Karst Studies*, 63(2), 59–66.

500

501

Gallegos, J. J., Hu, B. X., & Davis, H. (2013). Simulating flow in karst aquifers at laboratory and sub-regional scales using MODFLOW-CFP. *Hydrogeology Journal*, 21(8), 1749–1760.

502

503

Ghodrati, M., & Jury, W. A. (1992). A field study of the effects of soil structure and irrigation method on preferential flow of pesticides in unsaturated soil. *Journal of Contaminant Hydrology*, 11(1–2), 101–125.

504

505

Goldscheider, N., Meiman, J., Pronk, M., & Smart, C. (2008a). Tracer tests in karst hydrogeology and speleology. *International Journal of Speleology*, 37(1), 27–40.

506
507

Goldscheider, N., Meiman, J., Pronk, M., & Smart, C. (2008b). Tracer tests in karst hydrogeology and speleology. International Journal of Speleology, 37(1), 27–40.

508

509

Goldscheider, N., & Drew, D. (2007). Methods in Karst hydrogeology. International contribution to hydrogeology, IAH. Taylor & Francis/Balkema, London.

510

511

Haario, H., Laine, M., Mira, A., & Saksman, E. (2006). DRAM: efficient adaptive MCMC. Statistics and Computing, 16(4), 339–354.

512

513

Hauns, M., Jeannin, P.-Y., & Atteia, O. (2001). Dispersion, retardation and scale effect in tracer breakthrough curves in karst conduits. Journal of Hydrology, 241(3–4), 177–193.

514

515

Hull, L. C., Miller, J. D., & Clemo, T. M. (1987). Laboratory and simulation studies of solute transport in fracture networks. Water Resources Research, 23(8), 1505–1513.

516

517

Ibbitt, R. P., & O'Donnell, T. (1974). Designing conceptual catchment models for automatic fitting methods. Jahs Publication.

518

519

Johnston, P. R., & Pilgrim, D. H. (1976). Parameter optimization for watershed models. Water Resources Research, 12(3), 477–486.

520

521

Karay, G., & Hajnal, G. (2015). Modelling of groundwater flow in fractured rocks. Procedia Environmental Sciences, 25, 142–149.

522

523

Leven, C., Sauter, M., Teutsch, G., & Dietrich, P. (2004). Investigation of the effects of fractured porous media on hydraulic tests—an experimental study at laboratory scale using single well methods. Journal of Hydrology, 297(1–4), 95–108.

524

525

526

Levy, M., & Berkowitz, B. (2003). Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. Journal of Contaminant Hydrology, 64(3/4), 203–226.

527

528

Li, G., Loper, D. E., & Kung, R. (2008). Contaminant sequestration in karstic aquifers: Experiments and quantification. Water Resources Research, 44(2).

529
530

Majdalani, S., Chazarin, J.-P., Delenne, C., & Guinot, V. (2015). Solute transport in periodical heterogeneous porous media: Importance of observation scale and experimental sampling. *Journal of Hydrology*, 520, 52–60.

531

532

Majdalani, S., Guinot, V., Delenne, C., & Gebran, H. (2018). Modelling solute dispersion in periodic heterogeneous porous media: model benchmarking against intermediate scale experiments. *Journal of Hydrology*, 561, 427–443.

534

535

Massei, N., Wang, H. Q., Field, M. S., Dupont, J. P., Bakalowicz, M., & Rodet, J. (2006). Interpreting tracer breakthrough tailing in a conduit-dominated karstic aquifer. *Hydrogeology Journal*, 14(6), 849–858.

536

537

Mohammadi, Z., Gharaat, M. J., & Field, M. (2019a). The Effect of Hydraulic Gradient and Pattern of Conduit Systems on Tracing Tests: Bench-Scale Modeling. *Groundwater*, 57(1), 110–125.

538

539

Mohammadi, Z., Gharaat, M. J., & Field, M. (2019b). The Effect of Hydraulic Gradient and Pattern of Conduit Systems on Tracing Tests: Bench-Scale Modeling. *Groundwater*, 57(1), 110–125.

540

541

Mohammadi, Z., Illman, W. A., & Field, M. (2021). Review of laboratory scale models of karst aquifers: approaches, similitude, and requirements. *Groundwater*, 59(2), 163–174.

542

543

Morales, T., Uriarte, J. A., Olazar, M., Antigüedad, I., & Angulo, B. (2010). Solute transport modelling in karst conduits with slow zones during different hydrologic conditions. *Journal of Hydrology*, 390(3–4), 182–189.

544

545

Moreno, L., & Tsang, C. F. (1991). Multiple-Peak Response to Tracer Injection Tests in Single Fractures: A Numerical Study. *Water Resources Research*, 27(8), 2143–2150.

546

547

Pastore, N., Giasi, C. I., Allegretti, N. M., Redondo, J. M., Tarquis, M., & others. (2015). Experimental study of heat transport in fractured network. *Energy Procedia*, 76, 273–281.

548

549

Peely, A. B., Mohammadi, Z., & Raeisi, E. (2021). Breakthrough curves of dye tracing tests in karst aquifers: Review of effective parameters based on synthetic modeling and field data. *Journal of Hydrology*, 602, 126604.

550

551

Perrin, J., & Luetscher, M. (2008). Inference of the structure of karst conduits using quantitative tracer tests and geological information: example of the Swiss Jura. *Hydrogeology Journal*, 16(5), 951–967.

552
553

Qian, J., Zhan, H., Chen, Z., & Ye, H. (2011). Experimental study of solute transport under non-Darcian flow in a single fracture. *Journal of Hydrology*, 399(3-4), 246-254.

554

555

Qian, J., Zhan, H., Luo, S., & Zhao, W. (2007). Experimental evidence of scale-dependent hydraulic conductivity for fully developed turbulent flow in a single fracture. *Journal of Hydrology*, 339(3-4), 206-215.

556

557

Robert, C. P. (2015). The Metropolis-Hastings algorithm. John Wiley & Sons, Ltd.

558

Runkel, R. L. (1998). One-dimensional transport with inflow and storage (otis): a solute transport model for streams and rivers.

559

660

Sharmeen, R., Illman, W. A., Berg, S. J., Yeh, T.-C. J., Park, Y.-J., Sudicky, E. A., & Ando, K. (2012). Transient hydraulic tomography in a fractured dolostone: Laboratory rock block experiments. *Water Resources Research*, 48(10).

661

662

663

Tzelepis, V., Moutsopoulos, K. N., Papaspyros, J. N. E., & Tsihrintzis, V. A. (2015). Experimental investigation of flow behavior in smooth and rough artificial fractures. *Journal of Hydrology*, 521, 108-118.

664

665

van Genuchten, M. T., & Wierenga, P. J. (1976). Mass transfer studies in sorbing porous media I. Analytical solutions. *Soil Science Society of America Journal*, 40(4), 473-480.

666

667

Vrugt, J. A., Clark, M. P., Diks, C. G. H., Duan, Q., & Robinson, B. A. (2006). Multi-objective calibration of forecast ensembles using Bayesian model averaging. *Geophysical Research Letters*, 33(19).

668

669

Wang, C., Wang, X., Majdalani, S., Guinot, V., & Jourde, H. (2020). Influence of dual conduit structure on solute transport in karst tracer tests: An experimental laboratory study. *Journal of Hydrology*, 590, 125255.

670

671

Wu, Y., Hunkeler, D., & Goldscheider, N. (2019). Laboratory and numerical study of hyporheic flow-mediated DNAPL dissolution in karst conduits. *Hydrogeology Journal*, 27(1), 335-343.

672

673

Xiao, X., Xu, M., Ding, Q., Kang, X., Xia, Q., & Du, F. (2018). Experimental study investigating deformation behavior in land overlying a karst cave caused by groundwater level changes. *Environmental Earth Sciences*, 77(3), 1-11.

674

675
576

577

578

Zhao, X., Chang, Y., Wu, J., Li, Q., & Cao, Z. (2021). Investigating the relationships between parameters in the transient storage model and the pool volume in karst conduits through tracer experiments. Journal of Hydrology, 593, 125825.

579

580

Zhao, X., Chang, Y., Wu, J., & Xue, X. (2019). Effects of flow rate variation on solute transport in a karst conduit with a pool. Environmental Earth Sciences, 78(7), 1–18.

581
<figure><img src="image_18.png" /><figcaption>Figure A1. WSADE model fitting the 2 rpm experimental BTCs.</figcaption></figure>
$Q = 6.84 \times 10^{-4}$ L/s, pumping speed 20 rpm

<figure><img src="image_19.png" /></figure>

(a) Group 1, $h_1 \le h_2$

<figure><img src="image_20.png" /></figure>

(b) Group 2, $h_1 > h_2$

<figure><img src="image_21.png" /><figcaption>Figure A2. WSADE model fitting the 20 rpm experimental BTCs.</figcaption></figure>

<table><caption>Table A1. Calibrated parameters of WSADE to all of the experiments.</caption><thead><tr><th>Pump Speed (rpm)</th><th><i>M</i> (cm)</th><th><i>u</i><sub>1</sub> (m/s)</th><th><i>u</i><sub>2</sub> (m/s)</th><th><i>D</i><sub>1</sub> (m<sup>2</sup>/s)</th><th><i>D</i><sub>2</sub> (m<sup>2</sup>/s)</th><th><i>w</i><sub>1</sub> (-)</th><th>RMSE</th></tr></thead><tbody><tr><td rowspan="10">2</td><td>0</td><td>1.60×10<sup>-4</sup></td><td>1.23×10<sup>-4</sup></td><td>1.49×10<sup>-6</sup></td><td>3.99×10<sup>-6</sup></td><td>0.214</td><td>3.62×10<sup>-6</sup></td></tr><tr><td>1</td><td>1.61×10<sup>-4</sup></td><td>1.28×10<sup>-4</sup></td><td>8.85×10<sup>-7</sup></td><td>3.86×10<sup>-6</sup></td><td>0.149</td><td>4.20×10<sup>-6</sup></td></tr><tr><td>2</td><td>1.68×10<sup>-4</sup></td><td>1.36×10<sup>-4</sup></td><td>1.10×10<sup>-6</sup></td><td>2.46×10<sup>-6</sup></td><td>0.239</td><td>4.17×10<sup>-6</sup></td></tr><tr><td>3</td><td>1.68×10<sup>-4</sup></td><td>1.35×10<sup>-4</sup></td><td>9.75×10<sup>-7</sup></td><td>2.61×10<sup>-6</sup></td><td>0.227</td><td>3.93×10<sup>-6</sup></td></tr><tr><td>4</td><td>1.74×10<sup>-4</sup></td><td>1.42×10<sup>-4</sup></td><td>1.05×10<sup>-6</sup></td><td>3.05×10<sup>-6</sup></td><td>0.198</td><td>3.80×10<sup>-6</sup></td></tr><tr><td>5</td><td>1.67×10<sup>-4</sup></td><td>1.34×10<sup>-4</sup></td><td>1.06×10<sup>-6</sup></td><td>3.00×10<sup>-6</sup></td><td>0.160</td><td>3.51×10<sup>-6</sup></td></tr><tr><td>-1</td><td>1.77×10<sup>-4</sup></td><td>1.45×10<sup>-4</sup></td><td>8.49×10<sup>-7</sup></td><td>2.91×10<sup>-6</sup></td><td>0.165</td><td>4.09×10<sup>-6</sup></td></tr><tr><td>-2</td><td>1.80×10<sup>-4</sup></td><td>1.47×10<sup>-4</sup></td><td>9.35×10<sup>-7</sup></td><td>2.88×10<sup>-6</sup></td><td>0.162</td><td>3.83×10<sup>-6</sup></td></tr><tr><td>-3</td><td>1.80×10<sup>-4</sup></td><td>1.47×10<sup>-4</sup></td><td>1.20×10<sup>-6</sup></td><td>2.74×10<sup>-6</sup></td><td>0.175</td><td>4.14×10<sup>-6</sup></td></tr><tr><td>-4</td><td>1.79×10<sup>-4</sup></td><td>1.43×10<sup>-4</sup></td><td>8.37×10<sup>-7</sup></td><td>3.04×10<sup>-6</sup></td><td>0.123</td><td>4.32×10<sup>-6</sup></td></tr><tr><td>-5</td><td>1.76×10<sup>-4</sup></td><td>1.42×10<sup>-4</sup></td><td>9.57×10<sup>-7</sup></td><td>3.25×10<sup>-6</sup></td><td>0.135</td><td>3.58×10<sup>-6</sup></td></tr><tr><td rowspan="5">20</td><td>0</td><td>1.49×10<sup>-3</sup></td><td>1.09×10<sup>-3</sup></td><td>2.35×10<sup>-5</sup></td><td>5.95×10<sup>-5</sup></td><td>0.640</td><td>7.25×10<sup>-5</sup></td></tr><tr><td>1</td><td>1.48×10<sup>-3</sup></td><td>1.07×10<sup>-3</sup></td><td>1.67×10<sup>-5</sup></td><td>5.35×10<sup>-5</sup></td><td>0.654</td><td>3.46×10<sup>-5</sup></td></tr><tr><td>2</td><td>1.61×10<sup>-3</sup></td><td>1.29×10<sup>-3</sup></td><td>8.55×10<sup>-6</sup></td><td>6.00×10<sup>-5</sup></td><td>0.570</td><td>4.74×10<sup>-5</sup></td></tr><tr><td>3</td><td>1.60×10<sup>-3</sup></td><td>1.31×10<sup>-3</sup></td><td>8.46×10<sup>-6</sup></td><td>5.36×10<sup>-5</sup></td><td>0.558</td><td>5.15×10<sup>-5</sup></td></tr><tr><td>4</td><td>1.68×10<sup>-3</sup></td><td>1.37×10<sup>-3</sup></td><td>1.02×10<sup>-5</sup></td><td>5.73×10<sup>-5</sup></td><td>0.525</td><td>5.56×10<sup>-5</sup></td></tr></tbody></table>
<table><tr><td>5</td><td>1.57×10<sup>-3</sup></td><td>1.24×10<sup>-3</sup></td><td>9.52×10<sup>-6</sup></td><td>4.73×10<sup>-5</sup></td><td>0.581</td><td>4.53×10<sup>-5</sup></td></tr><tr><td>-1</td><td>1.71×10<sup>-3</sup></td><td>1.38×10<sup>-3</sup></td><td>1.46×10<sup>-5</sup></td><td>4.49×10<sup>-5</sup></td><td>0.507</td><td>8.51×10<sup>-5</sup></td></tr><tr><td>-2</td><td>1.83×10<sup>-3</sup></td><td>1.51×10<sup>-3</sup></td><td>6.67×10<sup>-6</sup></td><td>5.28×10<sup>-5</sup></td><td>0.277</td><td>5.57×10<sup>-5</sup></td></tr><tr><td>-3</td><td>1.83×10<sup>-3</sup></td><td>1.44×10<sup>-3</sup></td><td>7.84×10<sup>-6</sup></td><td>5.97×10<sup>-5</sup></td><td>0.362</td><td>3.66×10<sup>-5</sup></td></tr><tr><td>-4</td><td>1.73×10<sup>-3</sup></td><td>1.31×10<sup>-3</sup></td><td>9.76×10<sup>-6</sup></td><td>6.68×10<sup>-5</sup></td><td>0.473</td><td>5.30×10<sup>-5</sup></td></tr><tr><td>-5</td><td>1.69×10<sup>-3</sup></td><td>1.33×10<sup>-3</sup></td><td>9.59×10<sup>-6</sup></td><td>5.90×10<sup>-5</sup></td><td>0.508</td><td>4.69×10<sup>-5</sup></td></tr></table>

<figure><img src="image_22.png" /><figcaption>Figure A3. DRMIM model fitting the 2 rpm experimental BTCs.</figcaption></figure>
$Q = 6.84 \times 10^{-4}$ L/s, pumping speed 20 rpm

<figure><img src="image_23.png" /><figcaption>(a) Group 1, $h_1 \le h_2$</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b) Group 2, $h_1 > h_2$</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure A4. DRMIM model fitting the 20 rpm experimental BTCs.</figcaption></figure>

<table><caption>Table A2. Calibrated parameters of DRMIM to all of the experiments.</caption><thead><tr><th>Pump Speed (rpm)</th><th>M (cm)</th><th>u<sub>m,1</sub> (m/s)</th><th>D<sub>m,1</sub> (m²/s)</th><th>k<sub>1</sub> (1/s)</th><th>w<sub>m,1</sub> (-)</th><th>u<sub>m,2</sub> (m/s)</th><th>D<sub>m,2</sub> (m²/s)</th><th>k<sub>2</sub> (1/s)</th><th>w<sub>m,2</sub> (-)</th><th>w<sub>MIM,1</sub> (-)</th><th>RMSE</th></tr></thead><tbody><tr><td rowspan="12">2 rpm</td><td>0</td><td>1.77×10<sup>-4</sup></td><td>5.30×10<sup>-7</sup></td><td>2.36×10<sup>-4</sup></td><td>0.624</td><td>1.39×10<sup>-4</sup></td><td>2.78×10<sup>-6</sup></td><td>4.34×10<sup>-4</sup></td><td>0.110</td><td>0.876</td><td>3.72×10<sup>-6</sup></td></tr><tr><td>1</td><td>1.76×10<sup>-4</sup></td><td>2.56×10<sup>-7</sup></td><td>3.06×10<sup>-4</sup></td><td>0.434</td><td>1.34×10<sup>-4</sup></td><td>2.16×10<sup>-6</sup></td><td>4.83×10<sup>-4</sup></td><td>0.370</td><td>0.551</td><td>4.80×10<sup>-6</sup></td></tr><tr><td>2</td><td>1.84×10<sup>-4</sup></td><td>2.21×10<sup>-7</sup></td><td>4.00×10<sup>-4</sup></td><td>0.225</td><td>1.63×10<sup>-4</sup></td><td>1.45×10<sup>-6</sup></td><td>5.18×10<sup>-4</sup></td><td>0.598</td><td>0.255</td><td>6.13×10<sup>-6</sup></td></tr><tr><td>3</td><td>2.13×10<sup>-4</sup></td><td>5.11×10<sup>-11</sup></td><td>1.03×10<sup>-3</sup></td><td>0.134</td><td>1.68×10<sup>-4</sup></td><td>3.10×10<sup>-6</sup></td><td>8.72×10<sup>-4</sup></td><td>0.665</td><td>0.176</td><td>5.93×10<sup>-6</sup></td></tr><tr><td>4</td><td>1.95×10<sup>-4</sup></td><td>3.23×10<sup>-7</sup></td><td>3.87×10<sup>-4</sup></td><td>0.597</td><td>1.66×10<sup>-4</sup></td><td>7.62×10<sup>-8</sup></td><td>4.30×10<sup>-4</sup></td><td>0.169</td><td>0.767</td><td>5.00×10<sup>-6</sup></td></tr><tr><td>5</td><td>2.00×10<sup>-4</sup></td><td>4.06×10<sup>-8</sup></td><td>6.39×10<sup>-4</sup></td><td>0.358</td><td>1.58×10<sup>-4</sup></td><td>1.85×10<sup>-7</sup></td><td>5.17×10<sup>-4</sup></td><td>0.398</td><td>0.468</td><td>2.76×10<sup>-6</sup></td></tr><tr><td>-1</td><td>1.95×10<sup>-4</sup></td><td>7.05×10<sup>-7</sup></td><td>4.34×10<sup>-4</sup></td><td>0.511</td><td>1.65×10<sup>-4</sup></td><td>2.25×10<sup>-9</sup></td><td>5.01×10<sup>-4</sup></td><td>0.282</td><td>0.628</td><td>5.54×10<sup>-6</sup></td></tr><tr><td>-2</td><td>2.00×10<sup>-4</sup></td><td>3.19×10<sup>-10</sup></td><td>4.25×10<sup>-4</sup></td><td>0.283</td><td>1.94×10<sup>-4</sup></td><td>5.97×10<sup>-7</sup></td><td>4.40×10<sup>-4</sup></td><td>0.457</td><td>0.368</td><td>4.60×10<sup>-6</sup></td></tr><tr><td>-3</td><td>2.26×10<sup>-4</sup></td><td>3.32×10<sup>-9</sup></td><td>1.08×10<sup>-3</sup></td><td>0.192</td><td>1.88×10<sup>-4</sup></td><td>2.59×10<sup>-6</sup></td><td>6.95×10<sup>-4</sup></td><td>0.555</td><td>0.265</td><td>5.98×10<sup>-6</sup></td></tr><tr><td>-4</td><td>2.18×10<sup>-4</sup></td><td>1.64×10<sup>-8</sup></td><td>7.09×10<sup>-4</sup></td><td>0.326</td><td>1.45×10<sup>-4</sup></td><td>1.52×10<sup>-6</sup></td><td>1.98×10<sup>-4</sup></td><td>0.481</td><td>0.446</td><td>5.47×10<sup>-6</sup></td></tr><tr><td>-5</td><td>1.97×10<sup>-4</sup></td><td>2.72×10<sup>-7</sup></td><td>3.79×10<sup>-4</sup></td><td>0.621</td><td>1.29×10<sup>-4</sup></td><td>1.51×10<sup>-6</sup></td><td>1.52×10<sup>-3</sup></td><td>0.157</td><td>0.840</td><td>3.79×10<sup>-6</sup></td></tr><tr><td rowspan="12">20 rpm</td><td>0</td><td>1.65×10<sup>-3</sup></td><td>6.33×10<sup>-6</sup></td><td>7.01×10<sup>-4</sup></td><td>0.420</td><td>1.48×10<sup>-3</sup></td><td>6.13×10<sup>-6</sup></td><td>1.91×10<sup>-3</sup></td><td>0.375</td><td>0.578</td><td>2.82×10<sup>-5</sup></td></tr><tr><td>1</td><td>1.60×10<sup>-3</sup></td><td>7.67×10<sup>-6</sup></td><td>6.08×10<sup>-4</sup></td><td>0.389</td><td>1.38×10<sup>-3</sup></td><td>7.96×10<sup>-6</sup></td><td>4.04×10<sup>-4</sup></td><td>0.457</td><td>0.438</td><td>2.18×10<sup>-5</sup></td></tr><tr><td>2</td><td>1.63×10<sup>-3</sup></td><td>5.83×10<sup>-6</sup></td><td>3.16×10<sup>-4</sup></td><td>0.372</td><td>1.55×10<sup>-3</sup></td><td>1.38×10<sup>-5</sup></td><td>4.92×10<sup>-4</sup></td><td>0.488</td><td>0.458</td><td>1.87×10<sup>-5</sup></td></tr><tr><td>3</td><td>1.62×10<sup>-3</sup></td><td>4.55×10<sup>-6</sup></td><td>6.77×10<sup>-4</sup></td><td>0.420</td><td>1.57×10<sup>-3</sup></td><td>1.43×10<sup>-5</sup></td><td>2.26×10<sup>-4</sup></td><td>0.452</td><td>0.500</td><td>2.19×10<sup>-5</sup></td></tr><tr><td>4</td><td>1.77×10<sup>-3</sup></td><td>3.77×10<sup>-6</sup></td><td>1.95×10<sup>-3</sup></td><td>0.393</td><td>1.70×10<sup>-3</sup></td><td>1.90×10<sup>-5</sup></td><td>1.05×10<sup>-3</sup></td><td>0.447</td><td>0.428</td><td>1.54×10<sup>-5</sup></td></tr><tr><td>5</td><td>1.60×10<sup>-3</sup></td><td>5.18×10<sup>-6</sup></td><td>7.72×10<sup>-4</sup></td><td>0.441</td><td>1.54×10<sup>-3</sup></td><td>1.38×10<sup>-5</sup></td><td>1.79×10<sup>-4</sup></td><td>0.413</td><td>0.529</td><td>1.54×10<sup>-5</sup></td></tr><tr><td>-1</td><td>1.80×10<sup>-3</sup></td><td>3.73×10<sup>-6</sup></td><td>1.28×10<sup>-3</sup></td><td>0.408</td><td>1.62×10<sup>-3</sup></td><td>2.07×10<sup>-5</sup></td><td>6.37×10<sup>-5</sup></td><td>0.442</td><td>0.543</td><td>1.98×10<sup>-5</sup></td></tr><tr><td>-2</td><td>2.10×10<sup>-3</sup></td><td>8.44×10<sup>-6</sup></td><td>4.56×10<sup>-3</sup></td><td>0.351</td><td>1.84×10<sup>-3</sup></td><td>4.81×10<sup>-6</sup></td><td>7.14×10<sup>-4</sup></td><td>0.388</td><td>0.481</td><td>2.14×10<sup>-5</sup></td></tr><tr><td>-3</td><td>1.92×10<sup>-3</sup></td><td>1.31×10<sup>-5</sup></td><td>2.85×10<sup>-3</sup></td><td>0.377</td><td>1.83×10<sup>-3</sup></td><td>6.72×10<sup>-6</sup></td><td>3.70×10<sup>-4</sup></td><td>0.402</td><td>0.513</td><td>2.14×10<sup>-5</sup></td></tr></tbody></table>
<table><tr><td>-4</td><td>1.78×10<sup>-3</sup></td><td>4.84×10<sup>-6</sup></td><td>1.06×10<sup>-3</sup></td><td>0.458</td><td>1.69×10<sup>-3</sup></td><td>1.41×10<sup>-5</sup></td><td>1.64×10<sup>-4</sup></td><td>0.332</td><td>0.620</td><td>2.40×10<sup>-5</sup></td></tr><tr><td>-5</td><td>1.71×10<sup>-3</sup></td><td>6.22×10<sup>-6</sup></td><td>8.02×10<sup>-4</sup></td><td>0.516</td><td>1.65×10<sup>-3</sup></td><td>1.58×10<sup>-5</sup></td><td>2.03×10<sup>-4</sup></td><td>0.323</td><td>0.639</td><td>1.70×10<sup>-5</sup></td></tr></table>

# Appendix B. Parameter identifiability

We check the parameter identifiability of the DRMIM model by plotting the parameter histogram and correlation plots from the parameter sets on the MCMC chain. Correlation plots reveal the interaction between the parameters (Rana et al., 2019; Wang et al., 2019). Histograms allow direct estimation of uncertainty (relative standard deviation, RSD) and identifiability (whether the distribution is Gaussian). The Gaussian distribution of the posterior parameter sets is deemed essential to parameter identifiability (Younes et al., 2016) because the Gaussian distribution allows us to uniquely identify the parameter estimates from the experimental data.

The correlation between two data is evaluated by the correlation coefficient, $r$. Given paired data $\{(x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\}$, $r$ is defined as:

$$r = \frac{\sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{n} (x_i - \bar{x})^2} \sqrt{\sum_{i=1}^{n} (y_i - \bar{y})^2}}, \quad (10)$$

where $n$ is the sample size; $\bar{x}$ and $\bar{y}$ are the mean value for $x$ and $y$, respectively.

The parameter variability (or uncertainty) is represented by the relative standard deviation (RSD) value:

$$RSD = \frac{\sqrt{\frac{\sum_{i=1}^{n} (x_i - \bar{x})^2}{n-1}}}{|\bar{x}|}, \quad (11)$$

To better address parameter equifinality, it is optimal to investigate as many as possible the parameter sets which can reproduce the experimental data by certain criteria. However, the common practice of launching MCMC calibration for once may not be enough. Because, in our experience of applying the MCMC for optimization, the different initial parameter sets may attribute to different results: the starting point influences the path and the endpoint of the automatical searching algorithm. Different from the common practice, we carry out 20 times of MCMC calibrations for this identifiability study, to fully address the possible parameter space. Each MCMC calibration starts with different randomly generated initial parameter sets and contains 18000 iterations. The last 400 iterations of the 20 MCMC chains, thus 8000 parameter sets in total, are used in this identifiability study.
<figure><img src="image_26.png" /><figcaption>Figure B1. Parameter histograms and correlation plots of the MCMC chain of DRMIM fitting experiment of 5-5, 200 rpm.</figcaption></figure>

The results shown in Figure 17 suggest that different parameter sets can make the DRMIM model well reproduce the experimental BTC. The parameter histograms exhibit scattered distributions. The multi-solution phenomenon is especially severe for $u_{m,2}$ and $k_2$: there exist two separated groups of values (marked by the yellow bracket). With these different parameter sets, the DRMIM Region 1 generates various BTCs (the purple curves in Figure 17); Region 2 also generates various BTCs (the green curves in Figure 17). When the simulated BTCs of the two Regions are summed up, the simulated BTCs of the DRMIM model show little variance (the red curves in Figure 17) and all of the simulated BTCs are approximate to the experimental BTC.

Due to this equifinality problem, the calibrated DRMIM parameters may fall into a wide possible range. This has caused some calibrated DRMIM parameters to exhibit discontinuous variations with the *M* value (Figure 12b-e).