Introduction

Multipeaked breakthrough curves (BTCs) are not unusual when dealing with dye tracing. They can be explained by factors unrelated to the hydrogeological system itself, such as the injection mode (Brouyère and Rentier 1997) or varying flow conditions, precipitation events, weather changes, snow melt, etc. (Aydin et al. 2014; Goldscheider 2005; Göppert and Goldscheider 2008; Lauber et al. 2014; Liu et al. 2020; Perrin and Luetscher 2008; Poulain et al. 2018; Richter et al. 2022; Werner et al. 1997). Many research papers explain the occurrence of multipeaked BTCs in porous aquifers by heterogeneous hydraulic conductivity (Li et al. 2022; Siirila-Woodbrun et al. 2015), such as that found in stratified aquifers (Huyakorn et al. 1986), which allows for contrasted transit times. The combination of a rapid flow, e.g., along fractures (Li et al. 2023; Moreno and Tsang 1991), and a slow flow in the porous matrix allows the occurrence of a multimodal BTC. Thus, this explanation is also valid for karstic aquifers where hydraulic conductivity can be extremely variable. Authors mention the occurrence of multipeaked BTCs in karstic aquifers when conduit flow combines with a porous and/or fractured matrix flow (Goldscheider 2005; Lauber and Goldscheider 2014; Mudarra et al. 2014). Papers also report the occurrence of multipeaked BTCs in karstic rivers in case of large storage zones along the conduits (Dewaide et al. 2018; Sun et al. 2021; Tam et al. 2004; Tinet et al. 2019), which is supported by laboratory experiments (Field and Leij 2012; Zhao et al. 2017, 2021) and numerical simulations (Hauns et al. 2001). Additionally, Richter et al. (2022) found that tracer tests with particles in a karstic diffluent river can show a pronounced multipeaked BTC, while the same test with solute shows a single peak.

In this paper, the focus is on the occurrence of a multipeaked BTC downstream of a diffluence-confluence (DC) system where the tracer cloud is shared between multiple conduits of various lengths and flow velocities. This kind of system may generate two or more peaks on a BTC downstream, depending on the number of conduits, and this has been demonstrated either by field studies (Assunção et al. 2023; Cen et al. 2021; Goldscheider et al. 2008; Perrin 2004; Perrin and Luetscher 2008; Qi et al. 2018; Smart 1988; Wang et al. 2022b; Werner et al. 1997) or by laboratory experiments (Chu et al. 2021; Field and Leij 2012; Wang et al. 2020, 2022a, 2023). Similar behavior has been observed in glaciers where the drainage system can be complex (Fyffe et al. 2019). This behavior is strongly dependent on the characteristics of the DC system such as the respective lengths of the conduits, or the flow velocities and transport properties of each conduit (Field and Leij 2012; Wang et al. 2023). It is possible that a DC system cannot be highlighted by tracer tests if the multiple tracer clouds coming from the different streams arrive at the same time at the sampling location. In this case, the multiple peaks are merged, resulting in a single-peaked BTC (Chu et al. 2021; Field and Leij 2012; Wang et al. 2020, 2022a, 2023). A point of inflection can sometimes remain visible and is the last indication that a BTC is composed of multiple underlying peaks (Massei et al. 2006; Wang et al. 2022a, 2023).

As the studies about the influencing factors of the occurrence of a multipeaked BTC in divergent karstic rivers are mostly performed with laboratory designs, such studies with actual cave tracer tests are scarce. This study presents results from field tracer tests in the divergent karstic river of the Han-sur-Lesse Cave, in Wallonia, Belgium. It will focus on the influence of the DC system on a BTC downstream of the confluence and will attempt a quantification of the influence of various characteristics of the DC system on the occurrence of a multipeaked BTC. The methodology combines a dye tracing campaign and the use of a one-dimensional (1D) solute transport model. The rest of the paper is structured as follows. First, the methodology is described and is followed by the presentation of the results from the field tracer tests campaign and the solute transport model. Then, a discussion on the relative impact of various characteristics of the DC system on the resulting BTCs is provided. Finally, the main conclusions of the paper are presented.

Materials and methods

Two methods have been used: field tracer tests and solute transport modelling. This section first describes the study site, then the methodology for the field tracer tests campaign, and finally the solute transport modelling methodology.

Site: the Han-sur-Lesse Cave

The cave of Han-sur-Lesse is an internationally famous touristic cave (almost 300,000 visitors each year), located in the UNESCO Global Geopark “Famenne-Ardenne” in Wallonia, Belgium. The network has a total development of 10.7 km inside a folded Givetian limestones massif called “Massif de Boine” (Quinif and Hallet 2018). The karst in the massif is generally well-developed with the presence of large rooms and conduits (Bonniver 2011). The river Lesse flows inside the cave at the opening “Gouffre de Belvaux” and circulates inside the cave for about 1.7 km before it resurfaces at the opening “Trou de Han” (Fig. 1a). The annual average discharge of the river is 6.24 m3/s and can rise to more than 140 m3/s during storm events. The maximum capacity of the “Gouffre de Belvaux” is ~25 m3/s. When the discharge is higher, the old surface meander of the river is reactivated and quickly flooded. About 1.3 km after the “Gouffre de Belvaux”, a small portion of the river (about 0.2 m3/s in March 2022) flows west into a network called “Réseau Sud”, while most of the flow (~5.3 m3/s in March 2022) continues to reach the “Salles d’Armes” room. It is in this same room, 1.52 km after the “Gouffre de Belvaux”, that the river splits into two different streams, creating the main diffluence of the system. The eastern stream, called the Styx stream, flows into three sumps before reaching the confluence downstream at the “Salle d’Embarquement” room. The western stream, called the Dérivation stream (3.9 m3/s in March 2022), only activates when the Styx Stream reaches its maximum capacity (~1.4 m3/s in March 2022) and can reach a higher flow rate depending on the initial discharge of the river Lesse. Both streams have about the same length (~330 m) from the diffluence to the confluence. The river coming from the “Réseau Sud” (~0.2 m3/s in March 2022) flows into the Dérivation stream 170 m downstream of the diffluence.

Fig. 1
figure 1

a The diffluence-confluence system of the Han-sur-Lesse Cave in Belgium, and location of the fluorometers for the tracer tests campaign; b OTIS model reaches and their respective lengths and cell numbers

Full size image

Field tracer tests campaign

In March 2022, a tracer tests campaign was performed in the Han-sur-Lesse Cave aiming at studying the effect of a diffluence on a breakthrough curve downstream of the confluence. Four dye injections were done at various locations in the river, and seven field fluorometers were placed in the streams, as shown in Fig. 1a. One fluorometer IN (“inlet”) was placed a few meters upstream of the diffluence, and this will represent the “input” BTC entering the system for the injection 4. Two fluorometers were placed in each stream—FD1 and FD2 in the Dérivation stream, and FS1 and FS2 in the Styx stream, as shown in Fig. 1a—which will obtain two BTCs in each stream at various distances and will allow for calibration of the OTIS model reaches. One fluorometer was placed downstream of the confluence at the outlet of the system (OUT), which is where a double peak induced by the diffluence-confluence system might occur on the BTC. An additional fluorometer (RS) was placed at the outlet of the “Réseau Sud” network where a small portion of the river Lesse was flowing. This fluorometer is a control to check for any signal coming from this secondary system and is not directly related to the objectives of this study.

For the field tracer tests, the STREAM fluorometers were used (Poulain et al. 2017). They are designed with a RGB sensor and a blue LED (470 nm wavelength), which is the excitation wavelength of the fluorescein molecule. While the blue LED was on, the fluorometers were measuring the light intensity around the emission wavelength of fluorescein (520 nm) at an average frequency of one measurement every 8 s and saving the data on a SD card. All fluorometers were calibrated in the laboratory using the high-resolution methodology described by Deleu et al. (2021). The fluorometers are placed in a bath of water while measuring every minute. Uranine is slowly added to the bath using a peristaltic micropump, resulting in a high-resolution calibration curve, which is then fitted with a third-order polynomial with a R2 of more than 0.9999. A manual calibration is also done to account for any error linked to a bad calibration of the pump itself. This correction is done by a slight shift of the curve if necessary. The error of the pump is usually around 1%.

Four injections were done during three consecutive days, while the flow rate of the river Lesse was fairly constant with an average of 5.5 m3/s. For injection 1, 50 g of uranine were injected in the Dérivation stream. For injection 2, 70 g were injected in the Styx stream, whereas for injection 3, 20 g were injected in the Dérivation stream and 30 g were injected in the Styx stream simultaneously. The injections were done at the exact same locations as injections 1 and 2, just downstream of the diffluence. For injection 4, 200 g of uranine were injected in the “Gouffre de Belvaux” ~1.52 km upstream of the diffluence.

The BTC from the field fluorometers are then analyzed in terms of time of first detection, peak time, modal concentration, and calculated velocities. The velocities are calculated by dividing the estimated distance from the injection point by the time of first detection and the peak time, for the maximum and modal velocities respectively. The recovery rates cannot be calculated as no stream discharge can be measured in the cave due to field conditions. However, by assuming a 100% recovery rate, the stream discharge is estimated using the following formula:

$$Q=\frac{{M}_{\mathrm{inj}}}{{\int }_{{t}_{\mathrm{i}}}^{{t}_{\mathrm{f}}}C\left(t\right) dt}$$
(1)

where Q is stream discharge (m3/s), Minj = mass of injected tracer (g) and C is tracer concentration (g/m3), integrated through time between ti (time of first detection (i is initial)) and tf (time of the end of recovery (f is final)). The estimated discharge for each BTC is then compared with the stream discharge measured at the gauging station “Resteigne” 4.5 km upstream of the cave. The calculated flow rates are consistent (Table 1), and similar to the measured flow rate of the river (5.3 m3/s the day of the injections) at the gauging station “L5021/Resteigne” upstream the cave (SPW 2023).

Table 1 Tracer tests results
Full size table

In addition to tracer tests, bathymetric surveys of the river were done in some sections of the river by either manually measuring the depth along profiles, or by using a small (43 cm long, 28 cm wide) radio-controlled sonar boat. During the depth measurements, the exact location of the boat was measured using ultra-wideband technology (UWB) for indoor positioning (Van Herbruggen et al. 2019). This was done using UWB anchor nodes that were fixed on the cave walls. Their positions were measured with high precision by the remote sensing technology of T-LiDAR. A compatible UWB tag was also fixed on the boat itself and the distance to each anchor was measured continuously. This allowed for precise positioning of the boat (20 cm accuracy) many times per second during the bathymetric survey. The resulting depth measurements were processed into a point cloud depicting the river bottom. More details of the process and the sonar boat can be found in Van Rentergem and McFarlane (2017). An example of a bathymetric survey using the sonar boat in the Styx room is shown in Fig. 2:

Fig. 2
figure 2

Bathymetric survey of the room “Salle du Styx”, where the fluorometer FS1 is placed. The colored points are measured by the sonar boat while the gray area is the LiDAR scan of the room

Full size image

OTIS model

OTIS is a 1D solute transport model designed for streams and rivers (Runkel 1998). It is mostly used to quantify hydrologic parameters from BTCs. It is often used to model karst rivers that exchange with storage zones (Zhao et al. 2021). Specifically, OTIS was used by Dewaide et al. (2016) in the Han-sur-Lesse River. For this study, the model is used to estimate the dimensions of the conduits, quantify the tracer dispersion and retardation effect for different conduits, and separate a double peak downstream of the confluence based on the contribution of each stream.

Background and conceptual model

The OTIS model is based on a transient storage scheme, referring to the temporary storage of a portion of solute in a “storage zone”. The storage zone can refer to eddies, low-velocity zones or dead zones on the side of the stream; therefore, the model can be used to attempt the quantification of the retardation effect generated by the temporary detainment of solutes in storage zones along the stream. The model uses the advection–dispersion equations in which terms are added to consider for these transient storage processes. The tracer is considered conservative (no degradation or sorption) due to the short distances of the system (<2 km) and the low travel times (<8 h), and this is supported by the calculated flow rates which are in the range of the measured discharge (see section ‘Results’). For the model, the river is divided into different sections, called “reach”, which are then divided into a fixed number of cells. Each cell is attributed with the parameters of the reach it belongs to. The parameters and the main governing equations are as follows (Runkel 1998):

$$\frac{\partial C}{\partial t}=-\frac{Q}{A}\frac{\partial C}{\partial x}+\frac{1}{A}\frac{\partial }{\partial x}\left(AD\frac{\partial C}{\partial x}\right)+\frac{{Q}_{\mathrm{LIN}}}{A}\left({C}_{\mathrm{L}}-C\right)+\alpha ({C}_{\mathrm{S}}-C)$$
(2)
$$\frac{d{C}_{\mathrm{S}}}{dt}=\alpha \frac{A}{{A}_{\mathrm{S}}}(C-{C}_{\mathrm{S}})$$
(3)

where A is the main channel cross-sectional area (m2), AS is the storage zone cross-sectional area (m2), C is the main channel solute concentration (g/m3), CL is the lateral inflow solute concentration (g/m3), CS is the storage zone solute concentration (g/m3), D is the dispersion coefficient (m2/s), Q is the volumetric flow rate (m3/s), QLIN is the lateral inflow rate (m3 s−1 m−1), t is the time (s), x is the distance (m) and α is the storage zone exchange coefficient (s−1). This can be visualized in a simplified conceptual model (Fig. 3).

Fig. 3
figure 3

Conceptual model of OTIS and associated parameters (modified from Runkel 1998)

Full size image

In this study, the model for the diffluence-confluence system of Han-sur-Lesse Cave is divided into multiple reaches with their own length and number of cells. The length of the reaches (Fig. 1b) is mostly based on the placement of fluorometers, field observations and previous knowledge of the system (Bonniver 2011; Dewaide et al. 2016; Quinif and Hallet 2018). The number of cells is determined by the length of the reach (2 cells per meter). For the Dérivation stream, a first reach RD1 is established between the diffluence and the first fluorometer FD1 (Fig. 1b). Two more reaches RD2 and RD3 are then established between FD1 and FD2 fluorometers with an additional 30 m reach (60 cells) until the confluence. The separation between RD2 and RD3 is 30 m downstream of FD1 to allow the input of 0.2 m3/s from the “Réseau Sud” network. This interruption is in the form of a one-cell reach for which a lateral inflow QLIN = +0.2 m3/s is set. For the Styx stream, a first reach RS1 is established between the diffluence and the first fluorometer FS1 (Fig. 1b). A second reach RS2 is established between FS1 and FS2 with an additional 10 m reach (20 cells) to the confluence. For both streams, a one-cell reach is then added at the confluence to allow for the input of the flow rate of the opposite stream. Thus, for the Dérivation stream, a lateral inflow QLIN = +1.4 m3/s is set to represent the confluence with the Styx stream, and for the Styx stream, a lateral inflow QLIN = +4.1 m3/s is set to represent the confluence with the Dérivation stream. A final 110 m reach ROUT is established from the confluence to the outlet fluorometer for both streams. The two streams are then computed separately, while the last reach ROUT keeps its shared parameters for both streams. Observation points are added in both models at the location of the field fluorometers, which allows the observation of simulated BTC and the comparison with the BTC from the field fluorometers. This means that the observation point of the outlet (FOUT) has a signal for both models, which can then be superposed by simply summing up both curves to obtain one simulated BTC for the outlet, allowing for the separate observation of the contribution of both left and right streams to the outlet fluorometers. The inlet flow rates are 3.9 m3/s for the Dérivation stream and 1.4 m3/s for the Styx stream as described in section ‘Materials and methods’. The other parameters are then calibrated as described in the following.

Calibration and validation

For each reach and for both streams, the parameters A, AS, D and α must be calibrated. The calibration was done for injections 1 and 2 by manually changing the parameters of a reach and visualizing the resulting BTC at the fluorometer downstream of this reach. For each attempt, the simulated BTC was represented on the same graph as the field results to observe the match between the two BTCs. At each attempt, one of the four parameters was changed to affect the first time of arrival, the peak time, and the curve tail. This process was repeated until a visually satisfactory match between observed and simulated BTCs was found. The calibration implies the assumption that the tracer cloud is laterally homogenized across the section at the fluorometer locations. This is not granted for the fluorometer FS1 as the injections 1, 2 and 3 are done ~60 m upstream, which implies possible errors in the parameter values for the reach RS1. For the Dérivation stream, the high velocity (>1 m/s) and turbulence (Reynolds number > 3 million) of the first section (RD1) reduce the risk of lateral heterogeneities. After the calibration was done by obtaining a satisfactory match between observed and simulated BTCs for injections 1 and 2, a simulation of injections 3 and 4 was done, which allowed one to validate the calibration process by checking the validity of the parameters for both the “short” BTC (injection 3) and “long” BTC (injection 4).

Modified OTIS model

In this work, the results from field fluorometer BTCs combined with the OTIS curves are mainly considered in a first approach focused on the Han-sur-Lesse system itself. A second approach, focusing on attempting to derive some rules on DC systems in general, is done by changing some key parameters of the OTIS simulations (Fig. 4).

Fig. 4
figure 4

Schematic description of the modified parameters of the OTIS model. a The distance of injection is changed, b the distance of sampling is changed, c the streams’ lengths are changed, d the dispersion is changed

Full size image

Variation of the distance of injection

The distance between the injection point and the diffluence is changed (Fig. 4a). To do this, a new reach is added between the injection point “Gouffre de Belvaux” and the diffluence (IN). This reach is added in both OTIS models of the Dérivation and Styx streams and will share the same parameters, similar to the already computed reach from the confluence to the OUT fluorometer. To ease the process of the parameters’ calibration, an attempt at automatic calibration was made using the implemented OTIS-P algorithm (Runkel 1998; Scott et al. 2003). This tool involves the iterative estimation of optimal parameters based on a Nonlinear Least Squares (NLS) approach. The optimal parameters are reached when a set of convergence criteria is met. The length and number of cells are fixed as 1,520 and 3,040 m respectively, while the other parameters are automatically calibrated.

Variation of the distance of sampling point

To study the influence of the distance between the confluence and the sampling point at the outlet, the BTCs at various distances were simulated using OTIS (Fig. 4b). The methodology involves an increase in the length of the final reach (after the confluence) up to 1,000 m and the addition of multiple sampling points at various distances in the model. The simulated BTC from each point is then calculated by the mathematical superposition of both signals from the two models associated with each stream.

Variation of the system’s dimensions

To investigate the effects of a variation in travel time in both streams, OTIS simulations of the same system with modified geometry are done (Fig. 4c). For this, the dimensions of each stream are modified. First, the Dérivation stream’s length is shortened by 50%, by dividing the length and number of segments of the three reaches by 2. The Dérivation stream then becomes 162.5 m (70 m + 15 m + 77.5 m) in length instead of 325 m (140 m + 30 m + 155 m). Then, the Styx stream is modified in the exact same way by modifying the length of the two reaches. The downstream reach between the confluence and the outlet is unaltered for both cases. These changes are done for both injections 3 and 4.

Variation of transport properties

To investigate the importance of the transport properties of each stream on the final BTC, a modification of the intensity of the longitudinal dispersion is made on both streams for both injections 3 and 4 (Fig. 4d). First, a 50% reduction of the dispersion is made in the Dérivation stream, by simply dividing the initial value by 2 (see Fig. 5; 1.75 and 0.15 m2/s instead of 3.5 and 0.3) for all three reaches. Then, the same reduction is applied to the Styx stream in the exact same way for its two reaches.

Fig. 5
figure 5

Observed and simulated breakthrough curves of injection 1, and Dérivation stream model parameters

Full size image

Results

First, the results from the field tracer tests campaign in the Han-sur-Lesse Cave are presented. The measured and simulated BTCs are shown. Then, the results of the modified OTIS simulations are presented for each changed parameter.

BTCs from the Han-sur-Lesse Cave

The BTCs from the field fluorometers are presented and analyzed in terms of time of first detection, peak time, modal concentration, calculated velocities and flow rate. Then, the BTCs of the OTIS model for each of the fluorometers’ locations are presented and their similarities with the field data are discussed. The fluorescein concentrations are normalized to the injected mass to allow for a comparison between the injections. The results are presented for each of the four injections separately and are also available in a synthetized form in Table 1, where the time of first detection, the peak time, the modal concentration, the calculated velocities and flow rates are shown. The BTCs from the field fluorometers and OTIS simulations are shown in Figs. 4, 5, 6 and 7.

Fig. 6
figure 6

Observed and simulated breakthrough curves of injection 2, and Styx stream model parameters

Full size image
Fig. 7
figure 7

Observed and simulated breakthrough curves of injection 3. (The combined signal is the mathematical superposition of the OUT signals from the OTIS models of both the Styx and Dérivation streams)

Full size image

Injection 1: Dérivation stream

For injection 1, tracer is detected at both fluorometers located in the Dérivation stream (FD1 and FD2) and at the outlet (OUT) (Fig. 5). The time of first detection increases from 1 min 10 s at FD1 to 20 min 30 s at the outlet, and the modal concentration decreases from 2.1 ppb/g at FD1 to 0.15 ppb/g at the outlet. All three BTCs are asymmetrical due to retardation effects along the stream. The calculated velocities are high for the first section (1.9 m/s from the injection point to FD1), and lower for the next two sections (0.63 m/s between FD1 and FD2 and 0.35 m/s between FD2 and OUT).

For the OTIS model, the curves are very similar to the observed data from the field fluorometers (Fig. 5), which can be explained by the calibration of the model which was done based on injections 1 and 2. The parameters of the model were adapted for each reach so that the simulated curves match at best the observed data. The resulting parameters show a strong difference in terms of cross-sectional area; the first reach has a small cross-section area (3 + 2 m2) while the next two reaches have a much wider section area (13 + 7 m2). The last reach, downstream of the confluence, has a large section area of 70 m2. The longitudinal dispersion is much stronger in the first reach (3.5 m2/s) compared to the other reaches (0.3 and 0.05 m2/s).

Injection 2: Styx stream

For injection 2, tracer is detected at both fluorometers located in the Styx stream (FS1 and FS2) and at the outlet (OUT) (Fig. 6). The time of first detection increases from 4 min 5 s at FS1 to 57 min at the outlet, and the modal concentration decreases from 2.6 ppb/g at FS1 to 0.065 ppb/g at the outlet. All three BTCs are asymmetrical due to retardation effects along the stream. The calculated velocities are much lower than for the Dérivation stream (0.24 m/s from the injection point to the FS1 and ~0.13 m/s for the other two sections from the FS1 to the outlet).

For the OTIS model, the curves are very similar to the observed data from the field fluorometers for the same reason as for injection 1 (Fig. 6). The first two reaches have a much lower cross-sectional area (12 and 20.5 m2) than the last reach downstream the confluence (70 m2). The dispersion is stronger for the second reach (0.4 m2/s) than for the first and last reaches (0.05 m2/s). This could be explained by a possible underestimation of the stream length due to the occurrence of the sumps. As the length is measured on the available cave maps and the depths of the sumps are not precisely known, it is possible that the total length of this reach is underestimated, and a higher longitudinal dispersion is required to fit the field BTC.

For all reaches and for both injections 1 and 2, the calibrated cross-sectional area A + AS satisfactorily matches the measured cross-sectional areas from on-site bathymetric surveys. For example, cross-sections from the bathymetric surveys in the Styx room (Fig. 2) showed an average area of 13 m2, which is in the range of the calibrated cross-sectional area A + AS for the Styx stream (Fig. 6).

Injection 3: simultaneous (proximal injection)

For injection 3, tracer is detected at fluorometers located in both Dérivation and Styx streams (FD1, FD2, FS1 and FS2) and at the outlet (OUT) (Fig. 7). In the Dérivation stream, the time of first detection increases from 1 min 20 s at FD1 to 24 min 30 s at the outlet OUT, and the modal concentration decreases from 3.1 ppb/g at FD1 to 0.155 ppb/g at the outlet OUT. In the Styx stream, the time of first detection increases from 3 min 40 s at FS1 to ~60 min at the outlet OUT, and the modal concentration decreases from 2.4 ppb/g at FS1 to 0.053 ppb/g at the outlet OUT. At the outlet, a double-peaked BTC is observed. The first peak can be associated with the tracer cloud coming from the Dérivation stream, while the second peak is associated with the Styx stream. All BTCs are asymmetrical due to retardation effects along the stream. The calculated velocities are very similar to the velocities estimated from injections 1 and 2 (Table 1).

For the OTIS model, the curves calculated with the calibrated parameters seem to correctly match the observed BTCs (Fig. 7), thus making the validation of the model effective. This is especially true for the fluorometers just upstream of the confluence (FD2 and FS2), whereas for the fluorometers downstream of the diffluence (FD1 and FS1), the match is satisfactory but shows slight differences with field data, which might be due to the strong proximity to the injection. For such short distances (60 and 140 m), the curve shape is probably strongly sensitive to the transport parameters, which makes the calibration not easily reproducible. It is also possible that lateral homogenization of the tracer cloud is not yet fully done, which makes the observed BTC from field fluorometers less reliable. This is especially true for the FS1 fluorometer of the Styx stream, where on-site observations have shown that the cloud is not always fully dispersed longitudinally. At the outlet of the system, the double peak is well reproduced by the model. Both peaks satisfactorily match the observed data, with a slightly too high modal concentration for the second peak.

Injection 4: Gouffre Belvaux (distal injection)

For injection 4, tracer is detected at all fluorometers (Fig. 8). First, the tracer cloud reaches the diffluence (IN) after 2 h 30 min with a peak concentration of 0.038 ppb/g. The curve is quite symmetrical, suggesting that retardation effects in the portion “Lesse souterraine” are negligible compared to longitudinal dispersion. At the diffluence, the tracer is separated into two portions for each stream. In the Dérivation stream, the time of first detection is 2 h 34 min at the FD1 and 2 h 45 min at FD2. The modal concentration is the same as the BTC from IN (0.038 ppb/g) for both FD1 and FD2. In the Styx stream, the time of first detection is 2 h 38 min at FS1 and 3 h 25 min at FS2. The modal concentration is 0.035 ppb/g for FS1, which is similar to the inlet BTC (0.038 ppb/g), and 0.027 ppb/g for FS2. At the outlet OUT, a single-peak BTC is observed after 3 h 20 min and a modal concentration of 0.018 ppb/g. The calculated velocities range from 0.14 to 0.16 m/s (averaged since the injection point).

Fig. 8
figure 8

Observed and simulated breakthrough curves of injection 4. (The combined signal is the mathematical superposition of the OUT signals from the OTIS models of both the Styx and Dérivation streams)

Full size image

For injection 4, the OTIS curves correctly match the field data for all fluorometers (Fig. 8), thus making the validation of the model effective. Slight differences in modal concentrations are observed for FD2, FS1 and FS2. The peak of FS2 is slightly too early in the OTIS results. The OTIS curves for the outlet fluorometer OUT correctly match the field data, where the two peaks are merged into one single peak.

A peak is also observed at the RS fluorometer for this injection (Table 1). This is due to a small portion of tracer flowing by the “Réseau Sud” network (Fig. 1a). This portion of tracer is negligible as it has been estimated at ~8 g out of the 200 g injected, thanks to a flow rate measurement on-site. The tracer cloud coming from the “Réseau Sud” is observable at the fluorometer FD2 as a late low-intensity peak (0.29 ppb at ~8 h after injection). This is not considered relevant to the objectives of the study, therefore, the BTC is not shown in Fig. 8. This tracer signal has not been computed in the OTIS simulations as well for the same reasons.

Modified OTIS model

The BTCs obtained from the modified OTIS model (Fig. 4) are presented for each modified parameter.

Variation of the distance of injection

To simulate a variation of the distance of injection, a new reach was added between the injection point “Gouffre de Belvaux” and the diffluence (IN). The resulting simulated BTC at the inlet IN fluorometer for injection 4 has a visually satisfactory fit with the field data (Fig. 8). The calibrated parameters are as follows:

  • Main channel cross-sectional area A = 45 m2

  • Storage zone cross-sectional area AS = 5.5 m2

  • Dispersion coefficient D = 0.8 m2/s

  • Storage zone exchange coefficient α = 8.5 × 10–5 s−1

In Fig. 9, BTCs from various simulations are shown, each associated with a specific distance between the injection and the diffluence. As the distance increases, the BTC from the Styx stream slowly merges with the tail of the BTC from the Dérivation. At 152 m, the peak of the second BTC is not visible anymore, but an inflection point remains and suggests the presence of at least two superposed BTCs. However, after 760 m, the inflection point is not visible anymore, thus making the observation of two BTCs impossible.

Fig. 9
figure 9

Evolution of OTIS-simulated BTC at the outlet OUT as the injection is made further upstream from the diffluence. The distances are associated with reach lengths reduced to 0, 10, 20%, … of the original length (1,520 m). As the distance increases, a merging of the two peaks makes the BTC single-peaked, thus making the observation of a DC system impossible for distal injections (more than 760 m)

Full size image

Variation of the distance of the sampling point

In Fig. 10, BTCs from various simulations are shown, each associated with a specific distance between the sampling point and the confluence. Results show that an increase in the distance of the sampling point from the confluence induces a merge of the two peaks for the proximal injection. This merge makes the BTC become single-peaked after 600 m. Then, a point of inflection remains visible until 1,000 m where no clear sign of a DC system is visible.

Fig. 10
figure 10

Effect of the variation of the distance of the sampling point from the confluence (the distance is indicated above each curve) for the proximal injection. As the distance increases, a merging of the two peaks makes the BTC single-peaked

Full size image

Variation of the system’s dimensions

In Fig. 11, the results from the OTIS simulations with modified stream lengths show strong variations for the altered BTC. First, if the length of the Dérivation stream is reduced by 50%, the first of the two peaks at the outlet is shifted to the left on the time axis for both proximal and distal injections. For the proximal injection, this effect induces an even more contrasted double peak at the outlet (Fig. 11a). For the distal injection, the effect is much less visible as it only makes the BTC slightly more spread out along the time axis (Fig. 11b). Then, if the length of the Styx stream is reduced by 50%, the second of the two peaks at the outlet is shifted to the left on the time axis for both injections. For the proximal injection, this induces a merge of the two peaks into a single-peaked BTC (Fig. 11c). For the distal injection, the effect induces a shorter BTC and a higher modal concentration (Fig. 11d).

Fig. 11
figure 11

OTIS simulations of cases of the Han-sur-Lesse system with modified stream lengths. The stream lengths have a strong influence on the occurrence of a double-peaked BTC

Full size image

Variation of transport properties

In Fig. 12, the results from the OTIS simulations with modified transport properties only show slight variations of the outlet BTC for both injections. In fact, only a slight change is visible for the altered curve of the Styx stream for injection 3 (Fig. 12c). The Dérivation stream altered curve is only slightly modified (barely visible in Fig. 12a). The altered curves for the distal injection 4 show no perceptible differences from the original curves (Fig. 12b,d).

Fig. 12
figure 12

OTIS simulations of cases of the Han-sur-Lesse system with modified stream dispersion. The dispersion of each stream has little to no influence on the occurrence of a double-peaked BTC downstream of the confluence

Full size image

Discussion

Based on field results and OTIS simulations, the effects of a karst DC system on a BTC and on solute transport can be assessed. Thanks to the analysis of field results and initial OTIS model, it is possible to infer some rules about the effect of a DC system on a breakthrough curve specifically for the Han-sur-Lesse system. The results from the modified OTIS models can be used to attempt to derive more general rules about the effect of a DC system. This is done by considering changes in some of the key parameters of the model, and by the assessment of their effects on the resulting BTC at the outlet. Those changes in parameters are tentatively related to real-life scenarios such as a variation in the distance of injection, a variation of the DC system’s dimensions, and the variation in transport properties. First, the effects of the difference of travel times between the two streams is discussed. Then, the effects of the distances of upstream injection and downstream sampling is discussed. Finally, the effects of the longitudinal dispersion are discussed. The discussion is based on field results, simulated BTCs and the modified OTIS model results.

Effects of travel times contrast

The first two injections of the tracer tests campaign done in the Han-sur-Lesse Cave system indicate a difference in travel times between the Dérivation and Styx streams after the diffluence. The maximum flow velocity (averaged between the inlet and the outlet) is higher in the Dérivation stream (0.35 m/s) than in the Styx stream (0.12 m/s). This is explained by two morphological characteristics of the system. First, the first portion of the Dérivation stream (before FD1) has a very short travel time due to the small cross-section area of the stream (5 m2 according to the OTIS model). This is supported by field observations, as the flow velocity is generally very high in this portion of the river. This first section allows the tracer to reach the fluorometers FD2 and OUT faster than the tracer cloud from the Styx stream. Secondly, the three sumps in the Styx stream probably significantly increase the travel distance. The distance measured on the map of Fig. 1a is probably underestimated. Historical reports from cave divers show that the sumps can reach down to 50 m depth. This 3D configuration of the stream means that the travel distance in the Styx stream is higher than the Dérivation stream, and not similar as suggested by the map in Fig. 1a. Both those morphological characteristics can explain the differences in travel time from the diffluence to the outlet in both streams.

The differences in travel time are clearly demonstrated by the times of first detection for injections 1 and 2. The tracer cloud from the Dérivation stream arrives first after 20 min, while the cloud from the Styx stream arrives after ~60 min at the outlet OUT. This difference is responsible for the double peak at the outlet for injection 3. The first peak is due to the tracer cloud from the Dérivation stream, while the second peak is from the Styx stream with a higher travel time. This is supported by the similarities between the results of different injections. The results from injection 1 and the first peak of injection 3, both associated with the Dérivation stream, are very similar in terms of transport properties. The same applies to injection 2 and the second peak of injection 3, which are characteristic of the Styx stream. The OTIS simulations support this idea as well. In fact, the simulated first peak of injection 3 is exclusively from the simulation of the Dérivation stream, while the second peak is from the simulation of the Styx stream only. The mathematical superposition of both peaks fits the double peak, obtained by the outlet fluorometer OUT, well. It is important to note that injections 3 and 4 were not used for the calibration of the OTIS model. The good fit obtained for those injections indicates that the model can reproduce the double peak and thus explaining it.

The results from the OTIS simulations with shorter stream lengths show strong variations of the altered BTCs (Fig. 11). For the proximal injection, a reduction of travel distances for both streams, greatly affects the double peak. As the cloud coming from the Dérivation stream was early (time of first detection is 20 min at the OUT fluorometer) compared to the cloud coming from the Styx (57 min), a shortening of the Dérivation makes the double peak more contrasted (Fig. 11a), while a shortening of the Styx stream induces a merge of the two peaks (Fig. 11c). This shows that the travel time of each stream is an essential factor for the existence of a double peak downstream of a DC system; however, for the distal injection, the length reductions affect only slightly the shape of the single-peaked BTC at the outlet (Fig. 11b and d). This is due to the large BTCs entering the system. As the input curve becomes spread out along the time axis, the time variations induced by the shortening of the DC system’s streams are less significant, thus less visible on the output curve. The further the injection is made upstream of the diffluence, the less contrasted both curves become due to the large spreading of both BTCs along the time axis. This is clearly observable here, as a variation of travel time of each stream has much less impact on the final BTC for the distal injection. The effects of the distance of injection are discussed in the next section.

Effects of injection and sampling distances

Results from tracer tests in the Han-sur-Lesse Cave show that the effect of a DC system on a BTC is strongly dependent on the distance of the injection point from the diffluence. While a double peak is clearly observable for a proximal injection (Fig. 7), for the distal injection done from the “Gouffre de Belvaux” 1,300 m upstream from the diffluence, field results at the outlet downstream from the confluence do not show any double peaks (Fig. 8). Without further analysis, the experimenter would not infer the existence of a diffluence-confluence (DC) system from this BTC alone. In the same way, the OTIS simulation results with the mathematical superposition of the BTC from the Dérivation and Styx streams also do not show any double peak. The curve from the Styx stream still arrives late but it is “hidden” in the tail of the BTC of the Dérivation stream. This is possible due to the longitudinal dispersion of the 1,520 m conduit from the injection point to the diffluence. The longitudinal dispersion allows the curve to be more temporally spread out, making the curve larger. This makes the curve from the Styx stream arrive while the tracer cloud from the Dérivation stream is still passing by at the confluence, and thus makes the superposition of both curves appear as a single peak.

This effect is clearly observable in the modified OTIS results. By varying the distance of injection in the model, an attempt to observe the relationship between the merging and the distance of injection is made. As seen in Fig. 9, when the injection is made further upstream of the diffluence, the double peak gradually becomes merged as a single-peaked BTC at the outlet. The effect of the longitudinal dispersion on the merge of the peaks is observable; therefore, for the karstic system of Han-sur-Lesse and in high flow rate conditions (~5.5 m3/s), a double peak induced by the diffluence-confluence system is observable in the BTC downstream from the confluence if the injection is made at a distance less than ~150 m from the diffluence. Then, a point of inflection remains visible until 760 m where no clear sign of a DC system is visible. These distances are probably dependent on the longitudinal dispersion affecting the transported solute in the section from the injection to the diffluence. If the initial section was less dispersive, a higher distance might be necessary to make the two curves merge into a single peak. This merge is especially significant when the second peak has a much lower concentration than the first one, as it can then be masked in the tail of the first peak. This difference in peak concentrations can be favored by the dilution occurring at the confluence. If the two streams have contrasted flow rates, it is possible that one of the two peaks will be less intense downstream of the confluence. This is especially true when the tracer coming from the faster stream is already flushed out. In the case of Han-sur-Lesse, the second peak is associated with the Styx stream which has a lower flow rate (1.4 m3/s) than the Dérivation stream (4.1 m3/s). This makes the tracer cloud coming from the Styx stream sustain heavier dilution and thus it appears less intense on the OUT fluorometer. This can favor the occurrence of a single-peaked BTC even if the Styx arrives late, as the curve from the Styx can be masked in the tail of the first peak, especially if dispersion makes the curves spread out along the time axis.

In a similar manner, the distance between the confluence and the sampling point also influences the occurrence of a double-peaked BTC. Figure 10 shows that an increase in the distance makes the two peaks gradually merge until the BTC eventually becomes single-peaked, due to the longitudinal dispersion affecting the tracer clouds between the confluence and the sampling point. As the distance increases, both curves are spread out along the time axis, eventually leading to a merger of the two peaks. For the karstic system of Han-sur-Lesse and in high flow rate conditions (~5.5 m3/s), a double peak induced by the diffluence-confluence system is observable in the BTC downstream of the confluence if the fluorometer is located at a distance less than ~600 m from the confluence. Then, a point of inflection remains visible until 1,000 m where no clear sign of a DC system is visible. These distances are probably dependent on the longitudinal dispersion affecting the transported solute in the section from the confluence to the fluorometer.

Effects of transport properties

As discussed in the previous section, the longitudinal dispersion plays a significant role in the occurrence of a double-peaked BTC downstream of a DC system. The dispersion occurring before the diffluence can make the tracer cloud separate into the different streams spread out along the time axis. This allows for the two BTCs downstream of the DC to be merged into a single peak, which can be explained by the fact that the contrast of travel times between the two streams becomes “insignificant” compared the timescale of the input BTC. The same principle works for the dispersion occurring between the confluence and the sampling point; however, results from Fig. 12 suggest that the dispersion occurring inside the DC system in the separate streams has little to no effect on the outlet BTC. This strongly suggests that the transport properties of both streams can only slightly affect the downstream signal for both proximal and distal injections. Therefore, the transport properties of both streams are not a significant factor for the existence of a double-peaked BTC downstream of the diffluence. This is at least valid for DC systems like those in the Han-sur-Lesse Cave and can probably be explained by the relatively short lengths of both streams. Furthermore, the rapid flow in the Dérivation stream probably prevents the dispersion from being significant at the timescale of the tracer tests, in contrast to the Styx stream, which shows a slight variation when the dispersion is reduced (Fig. 12c). It is possible that an increase in dispersion could have an effect. Further analysis of the effects of dispersion could be done either in laboratory designs or in the OTIS model.

Conclusions

Multipeaked breakthrough curves (BTCs) are not unusual when performing dye tracing in karstic rivers. Even though such curves are commonly associated with the separation of the tracer cloud into different streams that converge back together downstream (Assunção et al. 2023; Cen et al. 2021; Goldscheider et al. 2008; Perrin 2004; Perrin and Luetscher 2008; Qi et al. 2018; Smart 1988; Wang et al. 2022b; Werner et al. 1997), they can also be due to the occurrence of a large body of water or a transient storage zone (Dewaide et al. 2018; Sun et al. 2021; Tam et al. 2004; Tinet et al. 2019). In the case of a diffluent river, it is possible for a multipeaked BTC to occur when the separated tracer clouds, flowing into each stream at the point(s) of diffluence, reach the sampling points at different times. This behavior appears to be dependent on the respective lengths, flow velocities, transport properties, and therefore the conduit geometries of each stream. The tracer test methodology (location of the injection and sampling points) also has a strong influence on the BTC. In this study, the goal was to attempt a quantification of the influence of those characteristics on a double-peaked BTC. This was done using a tracer tests campaign in the cave of Han-sur-Lesse (Belgium), combined with the use of a 1D solute transport model.

First, the results show a clear double peak when the injection is made at the diffluence. This is explained by the differences of flow velocities between the two streams, and by the presence of three sumps along one of the streams, thus increasing the difference of lengths. However, when the injection site is located further upstream from the diffluence, the two peaks slowly merge to eventually become one single peak. A point of inflection remains visible in the tail of the BTC, until eventually, no evidence of the occurrence of a diffluence-confluence (DC) system is observable on the BTC. For the system of Han-sur-Lesse, the merging of both curves makes the BTC become single-peaked when the injection is made at more than 150 m from the diffluence. A point of inflection is then visible until 760 m, where eventually no clear sign of a DC system is visible on the BTC. This is explained by the longitudinal dispersion affecting the tracer cloud from the injection point to the diffluence, spreading the BTC along the time axis, and thus favoring a merger of the two peaks. This is especially significant when the second peak has a much lower concentration.

Furthermore, the results show that the distance between the confluence and the sampling point is also a crucial factor for the occurrence of a double-peaked BTC. This is due to the longitudinal dispersion affecting the tracer clouds from the confluence to the sampling point. For the system of Han-sur-Lesse, if the sampling point is located further than 600 m from the confluence, the two peaks are merged enough so that the BTC becomes single-peaked. Then, a point of inflection remains visible until 1,000 m where no clear sign of a DC system is visible.

Finally, the results also suggest that the length of each stream is a crucial factor for the occurrence of a double-peaked BTC. A 50% reduction of the length of one of the two streams can induce a complete merger of the two peaks, thus making the observation of a DC system impossible. This is valid for a proximal injection (at the diffluence) but becomes rather insignificant for a distal injection (1,520 m from the diffluence). Finally, the results suggest that the transport properties of each stream only have a minor influence on the occurrence of a double-peaked BTC in the system of Han-sur-Lesse.

To conclude, a DC system can be highlighted by tracer tests, in the form of a double-peaked BTC or a point of inflection in the tail of the BTC if the following conditions are met:

  • The injection is made close enough to the point of diffluence, which is <150 m (double peak) or 760 m (point of inflection) for Han-sur-Lesse.

  • The sampling point is located not too far from the confluence, which is 600 m (double peak) or 1,000 m (point of inflection) for Han-sur-Lesse.

  • The two (or more) streams have contrasted travel times from the diffluence to the confluence, which is ~30 min for Han-sur-Lesse.

As the merging of the peaks is due to the longitudinal dispersion, the conclusions are also valid for other DC systems in karst areas; however, the values are site-specific. If the lengths and dimensions of each conduit are unknown, it is difficult to estimate if the injection is indeed done close enough to the diffluence or if the fluorometer is placed close enough to the confluence. As a double-peaked BTC due to a DC system is strongly dependent on the distances of the injection and sampling points, pointing out a DC system using dye tracing is a real challenge, which could explain why publications reporting on double-peaked BTCs are not so common in the literature. Studies of other DC systems would be of interest for comparison with these results, and to better quantify those effects on the occurrence of multipeaked BTCs in karst rivers.

As most dye tracer users do not know the travel times before performing a tracer test, the best recommendation the authors can give is to place the injection and the sampling points as close to the DC system as possible. A practical example would be for systems where caves give access to the rivers upstream of a resurgence. The most practical solution would be to place the sampling at the resurgence, but as this could lead to a merging of the multiple peaks, the user should investigate the idea of sampling in the cave.