Framework for assessing lateral flows and fluxes during floods in a conduit-flow dominated karst system using an inverse diffusive model

The aim of this study is to present a framework giving new keys to characterize the spatio-temporal variability of lateral exchanges for flows and fluxes in a karst conduit network during flood events. An inverse model using an analytical solution of the diffusive wave model is applied on data from two successive gauging stations to simulate exchange dynamics after recharge. The study site is the karst conduit network of the Fourbanne aquifer in the French Jura Mountains, which includes two reaches of 5-10 km characterizing the network from sinkhole to cave stream, and to the spring. The model is 5 applied after separation of the base and the flood components on discharge and total dissolved solids (TDS) in order to assess lateral flows and mass-fluxes and compare them to help identify water origin. Our results showed various lateral contributions in space between the two reaches located in the unsaturated (R1), and in both unsaturated and saturated zone (R2) as well as in time, according to hydrological conditions. Globally, the two reaches show a distinct response to flood routing, with important lateral inflows on R1 and large outflows on R2. By combining these results with mass flux exchanges and the 10 analysis of flood routing parameters distribution, we showed that lateral inflows on R1 are the addition of diffuse infiltration (observed whatever the hydrological conditions) and localized infiltration in the secondary conduit network (tributaries) in the unsaturated zone, except in extreme dry periods. On R2, despite inflows on the base component, lateral outflows are observed during floods. This pattern was attributed to the concept of reversal flows of conduit/matrix exchanges, inducing a complex water mixing effect in the saturated zone. From our results we build the functional scheme of the karst system. It demonstrates 15 the impact of the saturated zone on matrix/conduit exchanges in this shallow phreatic aquifer, and highlights the important role of the unsaturated zone on storage and transfer functions of the system.


Introduction
Hydraulic transfers and solute transport processes in karst aquifers are known very complex due to the organization of underground void structures leading to preferential drainage axis through a conduit network embedded in a less permeable fissured matrix (Kiraly, 2003;Ford and Williams, 2013).Flow processes are driven by the spatial variability of the geometrical elements that constitute the conduit network as full pipes, open channels, or pipe constrictions (Covington et al., 2009), leading exchanges in karst conduits in both unsaturated ("river stretches") and saturated zone ("flooded area").This study furthermore showed that the parametrization of solute transport in the system was strongly impacted by interactions within the saturated zone.However, all these works were designed for low flow periods and are therefore not suitable to investigate the temporal evolution of lateral flows during flood events in conduit-flow dominated karst systems with large quick flows.
To assess hydrodynamic processes, Saint-Venant equations (SVE) describe unsteady flow in partially filled conduits (Saint-Venant, 1871) and are generally used to simulate discharge in conduit-flow dominated karst aquifers.As the conduit flow in the unsaturated zone is mainly controlled by free-surface conditions, Manning's equations are favoured, although they underestimate head losses due to turbulent flows (Jeannin and Marechal, 1995).Equations to be used for the saturated zone should be adapted for pressurized conditions, such as the Preissmann slot model in ModBraC (Reimann et al., 2011), or the Darcy-Weisbach equation in pipe-flow models like the Storm Water Management Model (SWMM) (Jeannin, 2001;Campbell and Sullivan, 2002;Peterson and Wicks, 2006;Chen and Goldscheider, 2014).However, application of these models is questionnable, because detailed informations-often not available -on hydraulic parameters are required, particularly concerning the location and geometry of conduits.As a consequence, such physically-based models are often applied in a degraded mode (Le Moine et al., 2008), and are sometimes over-parameterized even in case of relatively well-known conduit networks.An alternative approach to model hydraulic processes would be to test the ability of simplified SVE with parsimonious parameters.
For that, following Moussa and Bocquillon (1996b), the diffusive wave equation can be considered as a simplification of the full SVE, and is even a higher order approximation than the uniform formulae (i.e.Manning's formula).This approach was used successfully by Charlier et al. (2015) to assess lateral flows in karst rivers with numerous lateral in-and outflows.It is therefore a promising alternative to the more complex approaches cited above.
One of the methods used to simplify the SVE, which consists of the continuity equation and the momentum equation, is the diffusion-wave method (DW) that neglects the acceleration terms (Moussa and Bocquillon, 1996b).Combined with Manning's equation or Chezy's equation, the DW equations can be simplified to one single equation (Moussa, 1996;Fan and Li, 2006;Wang et al., 2014).This equation can be applied to a wide range of phenomena in different fields as exposed by Singh (2002) for the kinematic wave.In fact, DW is largely used for flood routing model but also for solute transport, known as the advection-diffusion model (Cimorelli et al., 2014).With the assumption of constant values for wave celerity and diffusion coefficient (Hayami (1951) hypothesis) and in the particular case of uniformly distributed lateral flow, an analytical solution unconditionally stable of the Hayami model exists Moussa (1996).The inverse model needs as input the inflow hydrograph, the outflow hydrograph and a predetermination of the celerity and diffusivity; the inverse model enables to identify the temporal distribution of the lateral flow which can be either positive (lateral inflow) or negative (lateral outflow).The applied direct and inverse model of this analytical solution belong to the hydrological model MHYDAS (Distributed Hydrological Modelling of AgroSystems) and was described in Moussa et al. (2002).By using the inverse model, it is possible to simulate lateral inflows or outflows between two gauging stations during a flood event and to identify more accurately exchanges from lateral hydrogeological units along the main conduit of a karst system.
To model conservative solute transport along a 1D flow path, the advection-dispersion equation (ADE) is largely used in hydrology (Runkel, 1996;Baeumer et al., 2001) and karst hydrology (Hauns et al., 2001;Luhmann et al., 2012).Under steadystate flow conditions, and in agreement with the mass conservation law, solute mass fluxes can be expressed by ADE.However, ADE is more challenging to implement in case of unsteady-state flow conditions, especially when lateral exchanges occur.
In the present study, following (Singh, 2002), the diffusive wave model was applied for flux transfer, the flux being either discharge or mass flux.Under the hypothesis of uniformly distributed lateral fluxes, the analytical Hayami solution and the corresponding inverse model proposed by (Moussa, 1996) were used.Thereby, we present a new framework combining the diffusive wave model not only as a flood routing model, but also as a transfer function for mass propagation.Inverse modelling is used to identify lateral flows and estimate their concentrations coupling both, water flow and mass fluxes.The framework is tested on discharge and total dissolved solids data from two reaches of a karst conduit for several flood events under various hydrological conditions, with the aim to characterize lateral variability of flow and mass-flux exchanges.(Moussa and Bocquillon, 1996b): where Between two stations I (inflow) and O (outflow), the model is applied on the flood component (Q Based on the Hayami assumptions (Hayami, 1951), considering C Q and D Q as constant parameters over time along a channel network of length L, and without downstream boundary condition, the direct model of the diffusive wave equation without lateral exchange can be written as follows (Moussa, 1996): The symbol * representing the mathematical convolution function and the Hayami kernel function K Q (t), which is expressed as follows: (5) A first estimation of C Q and D Q can then be obtained comparing Q I,routed (t) to Q O,f lood (t).

Diffusive wave model with lateral flows
By considering the existence of lateral flow exchanges along a channel reach, we obtain (Moussa et al, 1996): where q(x, t) [L 2 .T −1 ] is the lateral flow rate per unit length as a function of distance along the channel reach x.The expression q(x, t) may be positive or negative depending on the occurrence of lateral inflow or outflow, respectively (Fig. 1).
An analytical resolution of the diffusive wave equation with lateral flows is proposed by Moussa (1996) based on the Hayami assumptions, accounting for uniformly distributed lateral flow between two gauging stations I (inflow) and O (outflow).
where dλ is the variable along which the integral function is calculated.
The lateral flood flows Q A,f lood (t) along the reach, under the hypothesis of uniformly distributed lateral flow, is expressed as: The inverse model enables the identification of the temporal distribution of the lateral inflows or outflows q(x, t) over the channel reach.According to Moussa (1996), by knowing Eq. 7 gives then: The resolution of equations ( 10) and (11) requires first the identification of K Q (t) using equation ( 5), and consequently a predetermination of the two parameters C Q and D Q , to calculate afterwards lateral flow Q A (t) as follows:

Assessing lateral mass fluxes
Classical 1D advection-dispersion equation (ADE) for steady-state flow conditions is analogous to the DW equation (Eq.1), replacing discharge by solute concentration, and celerity and diffusivity parameters by advective velocity and dispersion parameters, respectively.In unsteady-state flow conditions, and when lateral fluxes occurred, the application of ADE is not so obvious.We propose here to assess lateral mass fluxes (defined by Eq. 15) during flood, applying the DW model accounting for lateral exchanges, as a transfer function.We assume in that case that the mass conservation law is maintained and is equivalent to the volume conservation law in the diffusive wave model.Thus, the analytical solution of Moussa (1996) is used to resolve the conservative solute mass transfer accounting for uniformly distributed lateral fluxes.
Then, following the method described above for water flows, lateral solute mass exchange M A (t) is calculated by adapting Eqs. 5, 7, 10 and 11 as follows: The kernel function for the solute mass flux, K M (t), is expressed as: Finally, this modelling framework combining the calculation of both lateral water flows Q A (t) and solute mass fluxes M A (t), allows the assessment of the solute concentration of the lateral flows S A (t) in the case of positive or negative Q A (t) and

Framework
We propose in this section a step by step structure to help readers using our framework to investigate the exchange dynamics (2) Calculation of the total mass fluxes M I,tot and M O,tot using Eq. ( 15).(4) Calculation of the lateral base component for water flow (Q A,base (t)) and mass fluxes (M A,base (t)), using Eq. ( 14).
(5) Modelling of the lateral flood exchanges, using: A. first, the direct model (Eq.( 4)) to parametrize with a trial-and-error optimization, calculating Q I,routed and M I,routed to get the best fit with Q O,f lood and M O,f lood , respectively.The model is parsimonious with only the two parameters (celerity and diffusivity) to optimize.These two parameters are independent, hence the parametrization procedure using a trial-and-error optimization give similar results than a automatic optimization and converge rapidly to a unique and a stable solution.Nevertheless, if the user need to optimize more parameters, an automatic optimization procedure is necessary.
B. Afterwards, the inverse model is applied to simulate lateral flood exchanges Q A,f lood (t) and M A,f lood (t), using Eqs.
(7) Calculation using Eq. ( 21) of the solute concentrations of lateral base and flood flows S A,base (t) and S A,f lood (t), respectively.And determination, if possible, of the total solute concentration S A,tot (t).
In our study, the application of this framework is done separately on various selected flood event.Hence, the two parameters ) are calibrated for each event.Then, the relationships between the parameters and the variations of water flow and mass flux are analysed.It can be expected for example that C Q & D Q should increase with a floodpeak increase.Finally, the compilation of all results lead us to define a functional scheme of the studied karst system.
3 Study site

Field situation
The study site is located in the Doubs river valley at the northern limit of the Jura Mountains in Eastern France, near the village of Fourbanne (47  (Charmoille, 2005).The recharge area of the site covers about 30 km 2 .The site was selected for The climate is temperate with both oceanic and mountainous influence.Rainfall averages 1200 mm/year, occurring mainly in autumn and winter, but with slightly higher intensities in summer (Vermot-Desroches, 2015).Based on long-term records in the Jura Mountains, the number of rainy days per year is 140 on average (Charlier et al., 2012), corresponding for the most part to low-intensity rainfall events: 50% of rainy days had less than 3 mm of rainfall, whereas days between 15 and 30 mm of rainfall represented only 10% of rainy days.From the available time series, 7 flood events with complete data sets for all stations and various rainfall intensities were selected.A synthetic characterization of the 7 events is given in Fig. 4, where they are sorted from event 1 to 7 in function of decreasing baseflow at system outlet (station s3).The event duration varied from 24 to 52 hours with total precipitation amounts between 3.4 to 46.6 mm.A progressive decrease of the hydrogeological response with decreasing baseflow was found for events 1 to 5, whereas events 6 and 7 from an extremely dry period in June 2015 behaved differently with an important inflow from the station s1.

Field monitoring and data processing
The EC is directly related to the Total Dissolved Solids (TDS) assuming that TDS represent mainly conductive ionic compounds.The TDS values were therefore calculated directly from EC by using a constant factor of 0.64 (100 S.cm −1 = 64 g.l −1 ), which is commonly used by OTT CTD probes and consistent with the literature (Lloyd and Heathcote, 1985) according to the water mineralization range of the data set.

Model application to the study site
To illustrate the model behaviour described theoretically in Section 2 and to help to define a parametrization strategy, this section presents a sensitivity analysis on a representative event.We assume that the analysis carried out on a benchmark flood event is representative of the model parameters sensitivity on the other events.

Sensitivity analyses
This analysis was carried out on the celerity C Q and the diffusivity D Q for the direct and the inverse model, both applied on water flows.The celerity C Q and the diffusivity D Q describe the propagation and the flattening of the flood peak, respectively.
The model applied on fluxes using C M & D M is assumed to show an identical behaviour.(Moussa and Bocquillon, 1996b;Yu et al., 2000;Chahinian et al., 2006;Charlier et al., 2009) and confirm that the lower C Q and the higher D Q , the lower are the peak flow intensity and the transfer velocity.Figs.5a'-b' illustrate the simulation of lateral flows using the inverse model knowing input and output signals.It shows that C Q and D Q influence the lateral flow rates as well as the exchange direction (inflows or outflows) which could be inversed during a same flood event.Low C Q and D Q values lead to lateral inflows followed by outflows, whereas high C Q and D Q values lead to outflows followed by inflows.

Parametrization strategy
From the sensitivity analysis, a parametrization strategy was defined for the application of the direct model for the four parameters: C Q and D Q for discharge, and C M and D M for the mass fluxes.The parameters were optimized by the trial-and-error method based on simulations of the routed input signals and by using the following criteria: i) C was determined first from the peak delay between two succeeding gauging stations (celerity = reach length / delay), ii), whereas D was adjusted in order to get a best fit of the shape of the hydrograph or mass-chemograph with the observed output signal.

Simulation of lateral exchanges
Figure 6 illustrates the application of the framework (depicted in Fig. 2) for reaches R1 and R2 of the study site for flood event no. 1 with 21 mm of rainfall during high flow conditions.
The lateral base exchanges calculated for reach R1 demonstrate that the output signal observed at station s2 for discharge and mass flux cannot be entirely explained by the contribution of the input signal, but that it was essentially derived from lateral inflows between s1 and s2.The mass flux model shows in addition that the lateral inflows were strongly mineralized, with higher TDS values than for stations s1 and s2.The mass fluxes of the base component at station S2 were thus essentially derived from lateral inflows along reach 1.This pattern was different for reach R2, where lateral base inflows were 2-fold lower and with similar TDS values.Regarding the flood components, the simulated lateral exchange indicates important lateral inflow and high mass influx along R1.The concentration estimations indicate a similar evolution than for station s2, which is characterized by a TDS dilution during the flood.In contrast, the dynamic was totally different along R2, where outflows are simulated.The concentration estimations of the lateral component deduced from outflows and mass outfluxes are very low compared to the measurements at stations s2 and s3, suggesting the presence of more complex processes than simple outflows, as will be discussed later on in Section 5.2.

Variability of lateral exchanges
Following the example detailed for flood event no. 1 in the previous section, this section aims at summarizing results of the model application on all selected flood events in order to get informations on the general hydrological functioning of the field site.

Lateral exchanges for the base component
Figures 7a-c present lateral exchanges for base water flow and base mass fluxes for all events of reaches R1 (orange labels) and R2 (purple labels).When comparing the modelled mean lateral base flow exchange in function of the measured mean base input, two distinct linear relationships can be observed (Fig. 7a).For both reaches, lateral exchanges were positive, indicating inflows that increased linearly with average base input.For reach R1, lateral water inflow was 2.6 fold higher than mean input flow from station s1, whereas for reach R2 the mean lateral water inflow represented only 0.4 times the mean inflow from station s2.For mass fluxes (Fig. 7c) a very similar relationship was found.However, for R1 the slope of the correlation was much steeper than for water flow (4.3 against 2.6), indicating that the lateral inflow water was more mineralized than the input flow from station s1 already present in the system.In contrast, for reach R2 a lower slope was found for mass flux than for water flow (0.3 against 0.4), meaning that lateral inflow was less mineralized than input flow from station s2.In this Section we present the distribution of the model parameters -celerity and diffusivity -in function of maximum input flood flow intensities, because these data furnish information on transport dynamics (Fig. 8).
The C Q water flow celerity parameter increased linearly with input peakflow for events 1 to 5 for reaches R1 and R2 (Fig. 8a, graph in semi-log scale).Reach R1, entirely located in the unsaturated zone, showed lower celerities compared to reach R2, which was localized in both unsaturated and saturated zone.Events 6 and 7, corresponding to the extreme dry period, did not fit this trend for both reaches, and were characterized by lower C Q values.This behaviour may be related to the low degree of water saturation of the system at the beginning of the flood.The C M mass flux celerity parameter showed similarly a relationship with the flood flow intensities for events 1 to 5, but only for reach R1 and again with a different behaviour for events 6 and 7 (Fig 8c).On the contrary, we did not found any relationships for reach R2, suggesting a more complex behaviour for this part of the system, most probably related to the presence of the saturated zone.All events from reaches R1 and R2 had very low D Q diffusivities.Only events 1-4 on reach R2 showed increasing D Q values with increasing input peakflow (Fig 8b).We note that these 4 data points were all from a wet period and from reach R2, which localized in the non-saturated and the saturated zone.The D M mass flux diffusivity parameter showed a very similar distribution as D Q , i.e. again with the only elevated values for events 1-4 on reach R2 (Fig 8d).

Assessment of the saturated level of the conduit network
The analysis of the parameter distributions showed distinct patterns for R1 and R2, which can be attributed to the presence of the saturated zone in the lower parts of reach R2 (see Figure 3 for hydrogeological scheme of the main conduit).In fact, flood routing in the unsaturated zone is related to flow in open conduits, whereas in the saturated zone it is controlled by pressure transfer leading to an almost instantaneous propagation.Consequently, the difference of flood routing between unsaturated and saturated conduits along the system should be observed in the C Q parametrization of the 2 reaches R1 and R2.In fact, by considering a mean constant slope of the main karst conduit, we hypothesize a constant C Q value along the unsaturated zone in R1 and R2.The relation between C Q values of R1 and R2 for each event should therefore furnish information on the localization of the limit between the saturated and the unsaturated zone.Based on these hypotheses, the Eq.( 22) is used to estimate the percentage of conduit length of reach R2 located within the unsaturated zone, defined as "U" value. .100 This limit is supposed to fluctuate in function of the baseflow condition.Figure 9a represents the U values of all events in function of the mean baseflow at the system outlet in station s3 (Fourbanne spring).A linear relationship is observed for events 1 to 6 showing -as expected -an increasing U with a decreasing mean baseflow, which is coherent with the saturation level of the system measured by discharge.With an estimated total length of 5.4 km for R2, the saturated conduit length is approximatively 1.5 to 2.2 km.Event 7 shows that the system behaved differently for flood events during extreme dry periods.This point will be discussed in the next Section.
5 Discussion and conclusion

Modelling framework
Our study intends to present a new framework to quantify the temporal evolution of lateral flows and fluxes during floods in a well-developed karst conduit networks.It uses the diffusive wave (DW) model, which is a physically-based approach, parsimonious and easy-to-use.We used the model with an inverse approach allowing simulations of both lateral flows and mass fluxes under unsteady-state conditions following the assumption of uniformly distributed lateral exchanges.
Previous modelling studies focusing on the quantification of lateral exchanges in karst aquifers were conducted during base flow at steady-state flow conditions (Dewaide et al., 2016).In contrast, the framework of the present study furnishes new keys to identify such exchanges during flood events.It allows furthermore to quantify them, to estimate their concentration, and to characterize their origin(e.g.strongly mineralized lateral inflows from the rock matrix).An important advantage of the DW model over other approaches is its a parsimonious parametrization using only two parameters (celerity C and diffusivity D) for both flow and mass fluxes processes.The Hayami model also offers an unconditionally stable method with very low computing time, unlike numerical methods inducing numerical instabilities (Moussa and Bocquillon, 1996a, b).Finally, the short time- The proposed framework requires to decompose the base and the flood components, being consistent with the duality of flow processes dynamics.As proposed by many authors (Atkinson, 1977, among others), the base component typifies slow flow in the system, which is strongly influenced by interaction with the rock matrix (or low permeable volumes) with high storage capacities, whereas the flood component typifies quick flow within the conduit network with low storage capacities, which is strongly linked with flood intensity.Thus, considering that baseflow is used to estimate base exchanges and the diffusive wave flood routing model to assess flood exchanges, our framework gives a deconvolution of the two components of the lateral exchanges helping to interpret the involved processes.
Our methodology quantifies total lateral water flow and mass flux exchange for a given reach, but does not allow to identify simultaneously occurring local lateral flows and fluxes as already highlighted by several authors (Payn et al., 2009;Szeftel et al., 2011;Charlier et al., 2015).The results for the reach R2 furthermore denote the difficulty to decompose exchanges that occurred in part in the unsaturated and in part in the saturated zone.Lateral flood outflows for R2 were mainly combined with low flood outfluxes with total mineralisations lower than those observed for the input and output stations.This may be the result of complex exchanges occurring between conduit and matrix (or low permeable volume) water within the saturated zone, with out-and inflows occurring concomitantly during the same flood and having contrasted concentrations.The DW model thus furnishes numerous information about lateral exchanges in karst aquifers.Its limitation is that it considers a reach of a conduit system as a "homogeneous" unit with uniformly distributed lateral exchanges.It is a diagnostic tool which naturally cannot be used to decipher too much complexity in the case of too remote monitoring stations with highly variable lateral contributions.
As reported by Szeftel et al. (2011), we point out the limits of our framework to identify a single, unambiguous model structure representing transport processes when boundary conditions are poorly described.

Functional scheme
The modelling framework is proposed as a diagnostic tool to assess the dynamics of lateral exchanges between the main conduit and the neighbouring compartments of the karst aquifers during floods.The model application on several contrasted events on the Fourbanne karst network gives several points of discussions to evaluate exchanges in various hydrological conditions and allow to dress a general hydrological functional scheme of the site.This scheme, presented in Fig. 10, highlights the lateral (a) For major precipitations during high flow periods (events no. 1 and 2; fig10a) we showed the importance of lateral exchanges for both base and flood components.For the baseflow component both reaches were fed by higher mineralised lateral inflows.However, for the flood component the model yielded different results for the two reaches.Reach R1, located entirely in the unsaturated zone, obtained lateral inflows which were more mineralized than the sinking stream at the reach input.We relate these higher mineralized inputs to the arrival of tributaries from adjacent sinkholes (see their location on Fig. 3).On the contrary, reach R2 shows less mineralized base inflows and mainly slightly mineralized flood outflows, indicating high losses, which we relate to mixing processes in the saturated part of the conduit.
(b) For minor precipitations during low flow periods (events no. 3, 4 and 5; fig10b) reach R1 showed similarly mineralised inflows from the base and the flood components, but with a lower inflow amount in coherence with the lowest rainfall intensity of these events.On the contrary, the baseflow component of reach R2 was mainly characterized by inflows, whereas weak outflows were found for the flood component.This shows that the system was mainly influenced by the drainage of water from the rock matrix even during low flow periods.The weak outflows found for the flood component of reach R2 probably occurred within the saturated zone.
(c) For major precipitations during extremely dry periods (events no.6 and 7; fig 10c) a quite different behaviour was observed.Reach R1 presented very low base inflows with similar mineralizations as for the input flow, probably derived from sinking stream tributaries joining the unsaturated conduit network.However, important outflows were observed for the flood component, indicating high losses towards the rock matrix in the unsaturated zone.On the contrary, reach R2 presented very low baseflow inputs of highly mineralized water.The flood component was characterized mainly by lateral outflows followed by lowest inflows .It is interesting to denote this reversal of the lateral flood flows from outto inflows, indicating an evolution of the exchanges during the flood event, that we may be related to conduit/matrix interactions in the saturated zone.
Even if we observed increasing lateral inflows with increasing baseflow conditions, a slow inflow component remained always present all along the network (i.e. in the unsaturated as well as in the saturated zone), whatever the hydrological conditions (cf.section 4.2.1.).These constant inflows could probably originate mainly from diffuse infiltration in the unsaturated zone and from lateral drainage systems in the saturated zone.However, the mineralization of the infiltrating water was higher along R1 compared to R2, reflecting probably different recharge mechanisms.It seems that reach 1 collected additional inflows from strongly mineralized secondary tributaries mainly localized in the upstream part of the aquifer (R1), in accordance with the presence of sinkholes in the north-eastern part of the study area (Fig. 3).
wave model without lateral flows The diffusive wave equation is an approximation of the Saint-Venant equations used to model 1D unsteady flow in open channels ) Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-565,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 11 November 2016 c Author(s) 2016.CC-BY 3.0 License.where M (t) is the solute mass flux [MT −1 ], S(t) is the solute concentration [ML −3 ] and Q(t) is the discharge [L 3 T −1 ] As previously described for water flows, the model involves first the determination of the flood (M I,f lood (t) and M O,f lood (t)) and base (M I,base (t) and M O,base (t)) components of the total fluxes (M I,tot (t) and M O,tot (t)) based on mass-chemograph separation: of water flows and conservative solute mass fluxes along a channel reach between two gauging stations.To simulate lateral exchange flows Q A and fluxes M A during floods, required data are discharge and solute time series from both stations covering a complete flood event.The inverse modelling of the lateral flood component include four parameters corresponding to C Q and D Q for water flow and C M and D M for mass fluxes.A graphic representation of this framework is given in Fig. 2: (1) Collection of discharge (Q I,tot & Q O,tot ) and concentration (S I,tot & S O,tot ) data from input and output stations.
Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-565,2016   Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 11 November 2016 c Author(s) 2016.CC-BY 3.0 License.(3) Determination of the base and flood components by hydrograph (Q I,tot , Q O,tot ) separation defined by Eqs.(2) & (3) and mass-chemograph (M I,tot and M O,tot ) separation defined by Eqs (16) & (17).The base and flood components are separated with the constant slope method (McCuen, 2004) by using the inflection point on the hydrograph recession.The inflection point determined on the hydrograph is used to separate base and flood components for both hydrograph and mass-chemograph.
Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-565,2016   Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 11 November 2016 c Author(s) 2016.CC-BY 3.0 License.its dominantly allogenic recharge and its well-developed, partially accessible conduit network.The upstream recharge area corresponds to a surface watershed underlain by impervious lower Jurassic marls.The conduit network is fed by swallow holes located at normal faults, bringing into contact lower Jurassic shales and karstified middle Jurassic limestones (Fig.3b).The Verne swallow hole constitutes the main infiltration point of the Fourbanne karst system and corresponds to monitoring station s1.The allochtonous recharge joins the well-developed cave stream of the "En-Versennes" cave, which is completely explored over a distance of about 8 km.An artificial well near the village of Fontenotte gives direct access to the cave stream, where monitoring station s2 was installed, about 5 km downstream of station s1.The Fontenotte cave stream can be followed for another 2 km downstream of station s2, where it disappears into an inaccessible conduit network.It joins finally the Fourbanne spring, which is fed by a saturated siphon of 25 m depth and explored by cave-divers over the last 500 m.Monitoring station s3 was installed at the Fourbanne spring.
Discharge (Q) and electrical conductivity (EC) were monitored continuously at 15 minute intervals from January to June 2015 at the three stations (Fig.3b): the Verne swallow hole (station s1), the Fontenotte cave stream (station s2,) and the Fourbanne spring (station s3).Stations s1 and s2 were equipped with OTT CTD probes for conductivity, temperature, and water level monitoring, whereas a OTT-Hydrolab DS5X multiparameter probe was used at station s3 for conductivity and water temperature, and a OTT Orpheus mini probe for water level.The three stations divided the main conduit into two reaches: R1 from s1 to s2 (3.1 km) and R2 from s2 to s3(5.4 km, Fig 3c).Hourly precipitation were recorded by a fully automatic Campbell BWS200 weather station installed next to monitoring station s2 in the Fontenotte village.

Figs
Figs. 5a-a' represent the routed input hydrograph and the simulated lateral flow, respectively, varying C Q from 0.1 to 0.3 m.s −1 with a fixed D Q of 0.1 m 2 .s−1 (i.e.low D Q value).Similarly, Figs.5b-b' represent the same graphs, but with a fixed C Q of 0.15 m.s −1 (i.e.low C Q value) and varying D Q from 0.1 to 10 m 2 .s−1 .Figs. 5a-b illustrate the application of the direct model simulating the propagation of the input signal in order to fit the output signal with the two parameters set (C Q & D Q ) by considering no lateral exchange.The graphs clearly demonstrate that the routed input signal is much more sensitive to celerity than diffusivity.Varying C Q by a factor of 3 has a stronger impact on the results than varying D Q by a factor of 100.As an example, a 10% variation of C Q (with constant D Q ) modifies the maximum of Q I,f lood * K Q of 9%, while 10% variation of D Q (with constant C Q ) changes to vary the maximum of Q I,f lood * K Q of 0.6% only.Moreover, this sensitivity test illustrates the impact of both parameters on the propagation velocity and the shape of the flood peaks: increasing C Q values yield more rapidly propagating and higher flood peaks, whereas increasing D Q values lead to flattened peaks.These observations are in agreement with the literature(Moussa and Bocquillon, 1996b;Yu et al., Figs. 7b-d illustrate lateral water flows and mass fluxes calculated for flood flow.To summarize the dynamics of the lateral exchanges during the flood, minimum and maximum values are presented rather than average values in order to characterize intensities of both lateral gains and losses which may occur during the same event.Along reach R1, two distinct groups are observed depending on peakflow from input station s1 (Fig. 7b): (i) Events 1 to 5 with low input values are characterized by lateral inflows, with high maxima and minima close to zero.(ii) Events 6 and 7 with high input values from strong rain events during an extremely dry period, show with maxima close to zero and strongly negative minima indicating important lateral losses.Comparing the slopes of the linear relationships for water flows (Fig. 7b) and mass fluxes (Fig. 7d), it appears that the water inputs of the first group were strongly mineralized, whereas the water losses of the second group were characterized by low mineralization.For reach R2, all maxima of water flow and mass flux are close to zero, whereas the minima are negatively correlated with input values, indicating increasing lateral losses with increasing peakflow.Comparing the slopes between water flood flow (fig 7b) and mass flood flux (fig 7d), it appears that lateral losses were systematically less mineralized than input water from station s2.4.3 Transport dynamics towards the conduit network 4.3.1 Distribution of model parameters Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-565,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 11 November 2016 c Author(s) 2016.CC-BY 3.0 License.step analysis of the results show that the modelling framework may be used as a diagnostic tool to investigate the dynamics of lateral exchanges from flood analysis.Regarding parametrization, the analysis of the distribution of water flow parameters (celerity C Q and diffusivity D Q ) of 7 flood events allowed to characterize the variability of flood routing variability in more detail for the reach R1 located in the unsaturated zone and the reach R2 covering both unsaturated and saturated zones.Both reaches showed linearly increasing C Q values with increasing flood intensities (Section 4.3.1.).The difference of C Q parametrization between the two reaches was then used in Section 4.3.2 to quantify the fluctuation of the unsaturated/saturated boundary within the karst conduit.Mass flux parameters (C M and D M ) are less obvious to interpret.Reach R1 showed a linear relationship of C M in function of mass flux intensity, whereas a comparable trend was absent for reach R2 attesting the difficulty to characterize transport processes within the saturated conduit network from the available monitoring design.
Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-565,2016   Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 11 November 2016 c Author(s) 2016.CC-BY 3.0 License.contributions for the base and the flood components.In fact, we assume that -as discussed before -base and flood components are mainly related to specific processes as conduit/matrix exchanges distinguished by slow and fast flow, respectively.
Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-565,2016   Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 11 November 2016 c Author(s) 2016.CC-BY 3.0 License.unsaturated and saturated zone.One of the main points was the ability of our framework to propose a deconvolution of the output hydrograph as well as mass-chemograph allowing to quantify the lateral contributions in terms of flows and mineralization.Finally, as illustrated by the building of a conceptual model of the study case, our framework provided opportunities to experiment and investigate lateral flows and fluxes of karst conduit networks, known to have large amount of concomitant inand outflows.

Figure 1 .
Figure 1.Diffusive wave equation to model lateral flood flow exchanges along a channel reach of length L. The black curve Q I,f lood depicts the evolution of the flow rate with time at the beginning of the channel (input), and the blue curve Q O,f lood at the end (output).The dashed green curve Q A,f lood corresponds to the lateral flow exchanges which are positive for lateral inflows or negative for lateral outflows.

Figure 2 .
Figure 2. Framework to investigate lateral exchange dynamics of water flows and solute mass fluxes along a channel reach.

Figure 5 .
Figure 5. Sensitivity analysis of the model parametrization.Graphs (a-b) illustrate the direct model which simulate the routed input Q I,routed (dashed black lines) without lateral exchange, while graphs (a'-b') illustrate the inverse model which simulate lateral flows Q A,f lood (dashed green lines).Graphs (a-a') and (b-b') correspond to sensitivity test of CQ and DQ , respectively.

Figure 6 .
Figure 6.Framework application on the event no. 1 in high flow condition along the 2 reaches R1 (top) and R2 (bottom).Total flow and fluxes, base component and flood component, are represented in the first, second, and third column, respectively.

Figure 7 .
Figure 7. Base and flood analyses of the selected event set.Orange and purple symbols correspond to the lateral exchanges modelling along R1 and R2, respectively.Base analyses is performed on mean values ( ) exclusively whereas flood analyses take in account minimum ( ) and maximum ( ) values of the calculated lateral exchanges (Q A,f lood & M A,f lood ).Base (a) and flood (b) lateral flow (Q A,base & Q A,f lood ) are compared to the mean base and flood flow input (Q I,base & Q I,f lood ) and base (c) and flood (d) lateral fluxes (M A,base & M A,f lood ) are compared to the mean base and flood fluxes input (M I,base & M I,f lood ).

Figure 8 .
Figure 8. Parametrization analyses of the selected event set.Orange and purple symbols correspond to R1 and R2, respectively.Each parameter of the model is compared with the maximum intensity of the flood flow input signal (logarithmic scale).

Figure 9 .
Figure 9. Representation of the limit estimation between unsaturated and saturated conduit along R2 in function of the mean spring base flow.U calculation corresponding to the percentage of conduit length of R2 located within the unsaturated zone.