2.1. Lumped Parameter Modeling
Given the lack of spatial information on the structure of the karst reservoir, we propose a conceptual non-linear reservoir model. This model allows the simulation of flows in a karst environment while considering surface infiltrations and withdrawals from the different compartments of the watershed
We aim to propose a hydrogeological model based on the KarstMod modeling platform [
10], which can be modulated to simulate discharge at the outlet and to quantify the internal flows between the different compartments of the model and can be applied on various karst catchments [
11,
12,
13,
14]. KarstMod is developed by the Service National d’Observation du Karst (SNO Karst—INSU/CNRS). The proposed model comprises three compartments organized in two levels. The upper level corresponds to reservoir
E (Epikarst). It represents the unsaturated part of the system, which is also the location of a transient aquifer. This reservoir is connected with the two reservoirs of the lower level:
C (Conduit) and
M (Matrix). The
C conduit reservoir corresponds to the flows in the main karst conduits of the system and more generally to the rapid flows observed in parts of the system with higher hydraulic conductivity. Reservoir
M (matrix) accounts for flows in the fractured environments of the system, albeit in relation to the conduits. However, the flows here are slower. Reservoir
E controls the recharge of the system and so is subject to both precipitation and evapotranspiration. Reservoirs
M and
C are in connection with reservoir
E and may be subject to external forcing such as pumping (groundwater abstraction) or injections (river losses).
From a general point of view and in application of the laws of conservation of mass, the model is characterized by the following differential system to be solved at each time step:
with
E,
M and
C (L) water levels in reservoirs
E,
M and
C, respectively,
P (L T
−1): Precipitation rate,
ET (L T
−1): Evapotranspiration rate,
,
(L T
−1): Withdrawals from reservoirs
M and
C respectively, per unit area and finally
,
,
,
,
(L T
−1): Internal flow rates per unit area,
(
): Catchment area, Qs (L
3.T
−1): Watershed outlet flow.
In general, for a given reservoir
A (and therefore characterized by a level also noted as
A in mm) and its outlet noted as
S, the non-linear head loss law is expressed as follows:
where
(L T
−1) is the pressure drop coefficient between the two reservoirs,
A (L) is the water level in the reservoir, and
is an exponent that introduces a potential non-linearity in the model. Finally, we propose to introduce a flow between the conduit and the matrix (which can be positive if the flow is from reservoir
M to reservoir
C and it can be negative if the flow is from reservoir
C to reservoir
M). This flow integrates the interactions observed physically between these two compartments, which are particularly important during floods or during low water support [
12,
15]. Here, this internal flow noted
is expressed by the following expression:
where
(L T
−1) is pressure drop coefficient between the two reservoirs,
M and
C (L) are the respective water levels in reservoir
M and reservoir
C and
is an exponent that introduces a potential non-linearity of the model.
2.2. The Touvre Karst System
The Touvre karst system consists of a binary system (
Figure 1). The infiltration consists of a delayed infiltration of effective rainfall on the karstic recharge area, but also of direct infiltration of surface water of the Tardoire, Bandiat and Bonnieure rivers [
16,
17].
Figure 1.
Localisation of the Touvre karst aquifer and of the main hydrological, piezometric and meteorological stations. The percentages of different use of water over the four main sub-catchments (Touvre, Bandiat, Bonnieure and Tardoire) are also mentioned (stations’ acronyms: AGR: Agris, AMA: Amant, ANG: Angoulème-Brie-Champoniers, BRI: Brie, CHA: Chazelles, COUR: La Couronne, MART: Marthon, MON: Montignac, MONT: Montembœuf, RIV: Rivières and STF: Saint-Front).
Figure 1.
Localisation of the Touvre karst aquifer and of the main hydrological, piezometric and meteorological stations. The percentages of different use of water over the four main sub-catchments (Touvre, Bandiat, Bonnieure and Tardoire) are also mentioned (stations’ acronyms: AGR: Agris, AMA: Amant, ANG: Angoulème-Brie-Champoniers, BRI: Brie, CHA: Chazelles, COUR: La Couronne, MART: Marthon, MON: Montignac, MONT: Montembœuf, RIV: Rivières and STF: Saint-Front).
All the geological formations of the massif to the east of the Bajocian to Kimmeridgian springs are likely to develop karst networks. Moreover, the aquifer is well delineated geologically with a base corresponding to the impermeable marls of the Toarcian-Alenian and with a roof outcropping or sub-outcropping under the Quaternary alluvial formations, an Tertiary or Secondary (Cenomanian) continental detritic formations. While the karst is well described in its superficial part, it is almost unknown in its deep parts, even with deep drilling between 200 and 300 m that has revealed transmissive networks.
The limestone massif constituting the aquifer is well delimited. The eastern limit corresponds to the transition from the Toarcian-Alenian marl to the Bajocian limestone. The western boundary corresponds to the limit of the Upper Jurassic reef facies beyond which the formations become marlier, and therefore less likely to develop karstic networks. To the south-west, the limit corresponds to the Echelle fault. The northern and southern limits correspond to piezometric ridges.
The springs of the Touvre have three main outlets (the Bouillant, the Dormant and the Font de Lussac) and a secondary outlet (the Lèche) which is located at about 500 m from the first. They drain a massif of limestone of the Middle and Upper Jurassic, which extends from the northeast to the southeast of the emergences on the margins of the Massif Central. These outlets are located 7 km east of Angoulême.
The Bandiat and the Tardoire rivers have their source outside the sedimentary zone on impermeable crystalline terrain in the Massif Central. As soon as the rivers leave the impermeable terrain and enter the fractured limestone regions, losses occur. At present, the water infiltration occurs to such a degree that the waters of the two rivers are lost in their entirety for up to nine months a year. The Tardoire only joins the Charente river in periods of high water. The losses start at Feuillade on the Bandiat and at Montbron on the Tardoire.
The influence of the Bonnieure is less clear. In fact, in high water conditions, a piezometric peak can be observed between the Bonnieure and the Touvre watershed, probably reflecting the contribution of water from the surface water table contained in the sandy-clay formations. At low water, however, this ridge no longer exists. During these periods, there is probably a process of feeding the karst water table by the Bonnieure. Previous artificial tracing operations did not clearly demonstrate any relationship of the waters of the Bonnieure to the north with the sources of the Touvre, but this connection still seems a reliable hypothesis.
The Touvre watershed has been the subject of studies aimed at understanding its functioning. Larocque et al. [
18] proposed a systematic study using the tools of correlative and spectral analysis of the various hydrological, climatic and hydrogeochemical signals available in the study area. This analysis highlighted the strong hydraulic spatial heterogeneity of the La Rochefoucauld karst system. The La Rochefoucauld karst system highlights a relatively low level of karstification with a high regulation time of around 75 days.
The model presented by Le Moine et al. [
19] was not specifically designed for a karst watershed but aimed to consider potential losses to other catchments in the GR4J models. In particular, the authors applied this model to consider karst losses from the Bandiat and Tardoire catchments into the Touvre karst system. In addition, an extension of the study proposed a non-linear reservoir model of the Touvre integrating these losses as input. The results are encouraging, but the influence of withdrawals on low-water discharge in particular could not be quantified.
A neural network approach developed by Kurtulus and Razack [
20] for the period 1980–1986 made it possible to propose an original non-linear simulation of measure discharge based on precipitation and potential evapotranspiration data. However, beyond the learning phase of the artificial neural networks model, the prediction phases of the model indicated a degradation of the correlation between measured and simulated discharge, in particular for low discharge. However, this approach did not make it possible to quantify the impact of withdrawals on the watershed. Finally, Laroque et al. [
21,
22] also proposed a distributed numerical modelling approach for this watershed by exploring the transient aspects in particular.
Therefore, the Touvre watershed does not currently have a predictive model over long periods integrating a predominant karstic component between the different components of the massif and allowing the consideration of the withdrawals taking place on the whole of the karstic and non-karst impluvium of the watershed. This contribution aims to provide a first order estimation of the impact of the withdrawals based on a conceptual model in order to provide a first estimation to all the stakeholders of water supply over the region.
2.3. Data
According to long-term weekly measurements made by a hydrological station at Ruelle (within 1 km from Touvre resurgence) between 1895 and 1996, the minimum interannual weekly discharge is 1.4 m3/s, the average interannual discharge is 11.1 m3/s and the maximum interannual discharge is 40 m3/s. In order to properly address this issue, the daily time step is the most appropriate. In the context of this study, we relied on daily data measured a little upstream from the Touvre resurgence at Foulpougne.
The Touvre station at Foulpougne resurgence has been monitored since 1980. In addition, to consider the infiltrations on the karstic impluvium, we used the daily discharge data measured at the stations of Montbron for the Tardoire, Feuillade for the Bandiat and Saint Ciers for the Bonnieure, also extracted from the HYDRO database (
http://www.hydro.eaufrance.fr/ from 1 January 2016 up to 31 December 2018). Finally, in order to be able to estimate the losses of the Bandiat and the Tardoire during their passage on the limestone impluvium, the daily discharge data measured at the Coulgens station were introduced into the model. Indeed, the discharges measured at Coulgens are characterized during the summer period by a total drying up and therefore an infiltration of the rivers in the karst. During the high-water period, the infiltration via the losses is also active and the discharge increases at Coulgens station.
In order to estimate the effective rainfall, an estimation of daily rainfall was carried out based on rainfall records of up to eight meteorological stations managed by Meteo-France (
Figure 1). However, during the study period, daily rainfall data were not available for the whole period 2006–2018 and we therefore proposed an estimation with Thiessen polygons of variable size. A statistical analysis showed that this estimation did not lead to any change in the mean or variance of the series. For the estimation of daily ETP, we used Oudin’s formula [
23] for the available periods and then applied the Thiessen polygon method. As cultivated fields represent 25% of the watershed and local agriculture is mainly based on the culture of corn, a cultural coefficient of 1.1 was applied to produce an estimate of daily AET over the period from the 1st of June up to the 15th of September [
24]. Again, statistical analysis showed that this estimation with a variable number of stations did not lead to any change in the mean or variance of the series.
Concerning the estimation of anthropogenic abstractions, the annual abstraction data on the Touvre catchment area and the three other surface catchment areas, covering the period 2006–2018, were provided by the Adour Garonne Water Agency. In 2015, there were more than 175 declared abstraction points in all four watersheds (including 84 in the karstic impluvium of the Touvre). In the Touvre watershed, 4.6 Mm3 of water were withdrawn for agricultural needs and 1.1 Mm3 for domestic drinking water needs, mainly via boreholes. In the other three watersheds, 2.5 Mm3 were withdrawn for agricultural needs and 3.29 Mm3 for domestic drinking water needs, mainly via river intakes or alluvial groundwater drilling.
The total annual withdrawals on the watershed roughly correspond to 5% of the total water available at the outlet of the aquifer. This is quite low compared to the 50% of annual withdrawals considering the Lez spring for Montpellier city water supply [
25] but comparable with the 15% of annual withdrawals put in evidence by Sivelle et al. [
26] considering the Oeillal karst aquifer used for Narbonne city water supply.
At first glance, the 20 Mm3 of annual withdrawals may seem minimal compared to the 400 Mm3 supplied by the springs. However, these withdrawals are essentially concentrated in the summer months, when the flow of the springs is also low. Focusing on these three summer months, these withdrawals correspond to approximately 11.5 Mm3 for a volume of 50 Mm3. The main question is to evaluate the influence of these withdrawals on the flows available at the watershed outlet downstream where a strong conchicolous activity is present.
The application of the parsimonious model leads to several hypotheses concerning the daily internal forcing and in order to make both direct infiltration from surface rivers and groundwater abstractions. The daily internal forcing includes here the daily agricultural, industrial and domestic withdrawals from the conduit but also part of the discharge of the Bandiat and Tardoire rivers, which are located at the level of the contacts with the limestones of the Touvre watershed. With regard to the estimation of daily infiltration from the Bandiat and Tardoire rivers, field observation at Coulgens led us to propose an infiltration rate towards the conduit reservoir of 90% and an infiltration rate towards the matrix reservoir set at 10% on a permanent basis.
With regard to withdrawals, it was necessary at this stage to propose a model for the distribution of annual agricultural withdrawals on a monthly and then daily basis [
27]. It was necessary to specify over such a vast area, the level of abstraction on the scale of a catchment area as the accumulation of numerous individual decisions by each irrigator. We make the hypothesis that the agronomic optimum corresponds to irrigation management at the MTE (maximum evapotranspiration). The MTE is determined by the crop coefficient which evolves according to the development stage of the crop: MTE = crop coefficient × PTE. Following Allen et al. [
24], we set the crop coefficient at 1.1 for the period from June up to mid-September and its variation is assumed constant from one year to the next one. Irrigation is supposed to provide the plants with water to complement the rainfall and the soil reserves. Thus, the theoretical monthly water demand of the crop is evaluated by agro-climatic zone using a water balance on the RFU fixed at 60 mm. Water supply by irrigation is capped at 5 mm/day in relation to the available agricultural equipment. Daily withdrawals are estimated by dividing by the number of days in the month making the hypothesis of an overall stationarity of irrigation at the daily scale during low water periods. With regard to domestic and industrial withdrawals, daily withdrawals are simply taken as a uniform distribution of the annual withdrawal because no indication of under or over consumption has been detected in this watershed.
To quantify the rapid infiltrations in the model, we relied again on experimental field data, namely the artificial tracings carried out in the watershed at the level of the Bandiat and Tardoire losses. The various tracings carried out made it possible to highlight different transit times in high and low water. The speeds measured vary from 50 to 100 m/h depending on the place and time of injection, from 50 to 60 m/h in low water periods and from 100 to 120 m/h in high water periods (i.e., an overall transfer time of five to 10 days). From the estimated residence time distribution functions on these tracings, we have proposed infiltration transfer functions for the periods of high and low water. In the high-water phase from November to May, the transfer function includes a delay of five days and the sampling is then distributed over the four following days following a weighting {0.1, 0.2, 0.4, 0.2, 0.1}. In the low water phase from June to October, the transfer function includes a delay of 10 days and the sample is then distributed over the next four days using a weighting {0.5, 0.3, 0.15, 0.05}.
The daily internal forcing
includes the daily agricultural, industrial and domestic withdrawals from the matrix reservoir and a term corresponding to a potential diffuse infiltration of the Bonnieure river located to the north of the Touvre system. As far as the withdrawals are concerned, it is clear that the influence is not immediate for the system and even slower and more diffuse than for the conduits. While the injection points of the tracings do not allow us to deduce a clear connection between the Bonnieure and the Touvre, we have estimated, in view of the work of Laroque et al. [
18] in particular, that a supply to the karstic system of the Touvre was possible. In order not to add parameters and on the basis of the transmissivities estimated by Laroque et al. [
21], we therefore proposed to include a transfer function including a delay of 10 days followed by a constant infiltration over 10 days. In order not to overcomplicate the model, this transfer function was also applied to the withdrawals identified on the matrix reservoir.
Classically, the calibration criterion chosen in hydrological modelling is the Nash coefficient [
28], which corresponds to the mean square deviation between observed and simulated discharge. We tried many calibration criteria available in the literature (Pearson Product Moment Correlation, Akaike information criterion, Bayesian information criterion, Kling Gupta efficiency) but our choice finally fell on the non-parametric RKG criterion developed by Pool et al. [
29], defined by:
The rest of the study will corroborate the existing literature that shows that this criterion allows a good compromise in the optimization procedure and most often leads to the best results.