Predicting Sediment Discharge at Water Treatment Plant Under Different Land Use Scenarios Coupling Expert-Based GIS Model and Deep Neural Network

Excessive sediment discharge at karstic springs and thus, water treatment plants, can be highly disruptive. It is essential for catchment stakeholders and drinking water supplier to reduce the impact of sediment on potable water supply, but their strategic choices must be based on simulations, integrating surface and groundwater transfers, and taking into account possible changes in land use. Karstic environments are particularly challenging as they face a lack of accurate physical description for the modelling process, and they can be seen as a black-box due to the non-linearity of the processes generating 15 sediment discharge. The aim of the study was to assess the sediment discharge variability at a water treatment plant according to multiple realistic land use scenarios. To reach that goal, we developed a new coupled modelling approach with an erosionrunoff GIS model (WaterSed) and a deep neural network. The model was used in the Radicatel catchment (106 km2 in Normandy, France) where karstic spring water is extracted to a water treatment plant. The sediment discharge was simulated for five designed storm projects under current land use and compared to three land use scenarios (baseline, ploughing up of 20 grassland, eco-engineering, best farming practices). Daily rainfall time series and WaterSed modelling outputs extracted at connected sinkholes were used as input data for the deep neural network model. The model structure was found by a classical trial and error procedure, and the model was trained on two significant hydrologic years. Evaluation on a test set showed a good performance of the model (NSE = 0.82), and the application of a monthly-backward chaining nested cross validation revealed that the model is able to generalize on new datasets. Simulations made for the three land use scenarios suggested that 25 ploughing up 33 % of grasslands would not increase significantly sediment discharge at the water treatment plant (5 % in average). In the opposite, eco-engineering and best farming practices will significantly reduce sediment discharge at the water treatment plant (respectively in the range of 10-44 and 24-61 %). The coupling of these two strategies is the most efficient since it affects the hydro-sedimentary production and transfer processes (decreasing sediment discharge from 40 to 80 %). The coupled modelling approach developed in this study offers interesting opportunities for sediment discharge prediction at karstic 30 springs or water treatment plant under multiple land use scenarios. It also provides robust decision-making tools for land use planning and drinking water suppliers. https://doi.org/10.5194/hess-2020-363 Preprint. Discussion started: 22 September 2020 c © Author(s) 2020. CC BY 4.0 License.


Methodology
The proposed coupled modelling for the simulation of SD at WTP incorporates five different components that can be used separately or as part of an integrated modelling framework. These components are described in detail in the following sections. 110

Data handling
Appropriate data handling can help address various concerns in DNN modelling, such as their ability of generalization beyond their training limits (Kourgialas et al., 2015). Moreover, coupling GIS and DNN modelling can require important computational effort; therefore, the input data must be carefully selected. In this study, turbidity data distribution was explored using scatter plot (Fig. A1). The input turbidity data was chosen so that it includes the best dispersion of values (0-370 NTU) 115 and a strong recurrence of extreme turbidity values (> 150 NTU). Thus, two significant hydrologic years (H12-H13; from 10/1/1998 to 9/30/2000) were selected, accounting for 731 daily events. Thus, these two hydrologic years were selected as the reference period for the calibration and training of the model. SD (kg d -1 ) was assessed considering turbidity time series, mean pumping flow rate at the Radicatel WTP and a mean sediment concentration resulting from field investigations (~1 mg L -1 ; Hanin, 2011). In accordance with previous studies by Masséi et al. (2006) and Hanin (2011) in karstic environment in 120 Normandy, a lag of 1 day was applied to the SD time series to properly match the rainfall input. Erosion and runoff were calculated with WaterSed model for 269 events to reduce computational efforts, considering that erosion and runoff occurs only for significant rainfall events (P > 2.5 mm; threshold below which no effective rainfall is generated; and runoff and sediment discharge entering connected sinkholes was set to 0). https://doi.org/10.5194/hess-2020-363 Preprint. Discussion started: 22 September 2020 c Author(s) 2020. CC BY 4.0 License.

Erosion and runoff modelling 125
The WaterSed model uses a raster-based distributed approach to model the spatial distribution of runoff and soil erosion within a catchment for a given rainfall event. The WaterSed model is an upgrade of the STREAM model (Souchère et al., 1998;Cerdan et al., 2001) simulating hydrological processes by conceptualizing each raster grid cell as a reservoir whose properties are calculated at the event-scale and by routing both water and sediment according the surface flow network. The WaterSed model requires six datasets to compute runoff and erosion for any location in the catchment: 130  Digital Elevation Model (DEM) to extract slope and the runoff flow network,  Stream network (to burn the stream network into DEM) and river width,  Land cover and soil texture map to associate each land cover in the catchment with the appropriate soil surface characteristics,  Decision tables, adapted for the local conditions, to associate each soil surface characteristic observed in the land 135 cover or field with a steady-state infiltration rate, a Manning's roughness coefficient, a single potential sediment concentration and a soil erodibility value, and  Rainfall events including total precipitation, antecedent moisture (48h) and effective rainfall duration.
The DEM (5m, BD Alti®) depressions were filled according the algorithm developed by Wang and Liu (2006). The stream network location and width were provided by the BD TOPO®. A land cover map was developed for 2016 by combining two 140 national databases: the French Land Parcel Identification System (RPG) and the Soil Observatory of Upper Normandy (OSCOM). A soil texture map with three classes: clay, silt, and sand were derived from the Regional Soil Referential (RRP) with a precision of 1:250000. The parametrization consists of a characterization of the main parameters influencing runoff and infiltration in the studied context: soil surface crusting, surface roughness and crop cover. These soil surface characteristics were defined for each month and for each crop class according to the cropping calendar developed by Evrard et al. (2010) and 145 Delmas et al. (2012). Steady-state infiltration and potential sediment concentration were assigned to each soil surface characteristics according to Cerdan et al. (2001) and Cerdan et al. (2002b). The WaterSed model requires Manning's roughness coefficient in order to compute flow velocity. Based on soil surface characteristics, Manning's values were derived from surface roughness (Morgan, 2013) and the percentage of crop cover (Gilley et al., 1991). Last, soil surface characteristics were also used to define the soil erodibility factor, varying in the [0-1] interval, by adapting table developed by Souchère et al. 150 (2003a). Calibration parameters were extracted from a previous study (Landemaine, 2016;Landemaine et al., 2020b) conducted in the Austreberthe catchment, located 30 km east of the Radicatel catchment.

DNN Configuration
In this study, a multi-layer feed forward network (DNN) was built under Python version 3.6 using the high-level API Keras (Chollet, 2015) and Tensorflow (Abadi et al., 2016) as a background engine. We used the rainfall time series (i.e. daily 155 cumulative rainfall (P); and 48h-antecedent rainfall (P48)) retrieved from SAFRAN database, and WaterSed modelling outputs (runoff (RWS) and sediment discharge (SDWS)) at connected sinkholes as input data for the DNN model. SD time series at the Radicatel WTP was selected as output considering the two significant hydrologic years (H12-H13). Data (time series length n = 731) were split as a training set (70 %; n = 511) and a test set (30 %; n = 220). The training set was also split into training and validation set to fine-tune the hyperparameters (70-30 %; n = 358-153). The optimum structure and configuration of the 160 network was found by a classical trial and error procedure (training-evaluation process through optimization of errors; Ortiz- Rodriguez et al., 2013). When performing the DNN training, the number of iterations epoch was set to 15, and the batch size to 1. Stochastic gradient descent (batch size = 1) was used to calculate error and update the model for each example in the training dataset. The number of hidden layers and neurons depends on the characteristics of the input data, and there is no specific rule for selecting these parameters (Le et al., 2019). The final structure of the DNN was composed of one input/output 165 layer and three hidden layers. We set 41 neurons in the hidden layers as follow: 30-10-1. Before training, input variables were rescaled between 0 and 1 using normalization and re-transformed for the predictions using the normalization parameters. We used the rectified linear (ReLU) as activation function for the first two nodes and linear activation function was set in the last hidden layer. Mean square error (MSE) was chosen as loss function and the Adaptive moment estimation (Adam) was adopted as the model optimization algorithm. 170

Performance evaluation
The choice of the training/test set remains arbitrary and may need to be evaluated to compute a robust estimate of model error.
With time series data, care must be taken when splitting the data in order to prevent data leakage. To address this issue, we adapted a month-backward chaining nested cross-validation procedure, which provides an almost unbiased estimate of the true error (Varma and Simon, 2006). The procedure contains an inner loop for parameter tuning, and an outer loop for error 175 estimation (see Fig. B1). The parameters that minimize error are chosen on the inner loop and we add an outer loop, which splits the data into multiple training/test sets. Then, the error is averaged on each split to evaluate the overall performance of the model. The performance of the model was evaluated through a range of statistical error measures including the coefficient of determination (R²), the Nash-Sutcliffe Efficiency (NSE), and the root mean square error (RMSE). The coefficient of determination (R²) ranges from 0 to 1 and describes the amount of observed variance explained by the model. A value of 1 180 suggests that the model can explain all the observed variance, whereas a value of 0 implies no correlation.
Where is the observed variable at time , ′ is the predicted variable at time , ̅ and ̅ ′ the mean values of observed and predicted variable at time , and the number of observations.
The root mean square error (RMSE) corresponds to the standard deviation of the residuals (prediction errors). The residuals 185 are a measure of the distance from the data points of the regression line. It evaluates how the predictions match to the observations, and values may range from no fit (+∞) to perfect fit (0) based on the relative change of the data.
The Nash-Sutcliffe Efficiency (Nash and Sutcliffe, 1970), indicates the model's ability to predict variables different from the mean, and gives the proportion of the initial variance accounted by the model. NSE values vary between -∞ (poor model) and 190 1, indicating a perfect fit between observed and predicted values.

Designed storm projects and land use scenarios
Five designed storm projects were constructed for the modelling approach based on expert knowledge and depth-durationfrequency curves of the French Meteorological Survey (Météo-France) available from 1996-2006 on the studied area (Table  195 1). We considered winter events for low return period rainfall (0.5 and 2 year-return period). In the studied region, low daily cumulative rainfall can lead to severe erosion by water runoff, due to saturated soils induced by heavy cumulated rainfall on the antecedent days (Le Bissonnais et al., 2002;Evrard et al., 2010). For higher return period (i.e. 10, 50, and 100 year-return period) we considered spring events that are characterized by stronger rainfall intensities and can be particularly damaging in this region (Evrard et al., 2007). 200 were incorporated in the model in order to predict SD variability at the WTP and evaluate their impacts. All scenarios are  the catchment (increasing crop cover or delaying the formation of the slaking crust). 50 % of plots were randomly selected and we applied a 15 % increase in infiltration capacity, respecting the actual proportions of winter and spring crops on the catchment. The applied value was set based on the study of Maetens et al. (2012), who synthetized the reduction in erosion and runoff following different agricultural practices across Europe. According to their results, the 15 % increase in infiltration capacity is a low case assumption easily achieved through simplified cultural 225 techniques (e.g. minimum tillage, no till, direct seeding, crop cover, etc.).

Modelling the inputs
In order to fed the DNN, erosion and runoff were simulated with the WaterSed model for 269 rainfall events over the entire Radicatel catchment. We used the WaterSed parametrization that was carried out on the adjacent catchment (La Lézarde; 230 Landemaine et al., 2020a) and that already proved to be valid on the same area of the Pays de Caux (Landemaine, 2016;Landemaine et al., 2020b) and Belgium (Baartman et al., 2020). For each event, WaterSed outputs (i.e. runoff (RWS; m 3 ) and sediment discharge (SDWS; kg d -1 ) values) were extracted and summed over all positively connected sinkholes.

DNN: Calibration and Generalization
We used the runoff and sediment discharge predicted by the WaterSed model for the two selected hydrologic years. Predicted values of runoff and sediment discharge were extracted over the connected sinkholes and summed to be used as inputs for the 245 DNN. The rainfall time series (P, P48) available on the catchment were also considered as inputs. Turbidity time series recorded at the WTP was transformed into sediment discharge knowing the mean pumping flow rate at the Radicatel WTP, and considering a mean sediment concentration of 1 mg L -1 from field investigations made by Hanin (2011). Thus, sediment discharge at the WTP was considered as output of the DNN. The DNN was trained from October 1998 to February 2000 (n = 511 daily events). The remaining period from March 2000 to September 2000 was used as the test set (n = 220 daily events). 250 Modelling results (Fig. 4) suggested a good agreement between observed and predicted SD at the WTP over the test set, with a NSE criterion reaching 0.82 and a RMSE around 383 kg. The temporal variability was also well reproduced. The results over the training set were slightly less efficient and reached a NSE of 0.6 and a RMSE of 420 kg. Overall results were above the classical threshold values for satisfactory model, defined between 0.5 < NSE < 0.65 in hydrological studies (Moriasi et al., 2007;Ritter et al., 2013). Over the investigated period, the cumulative SD observed at the WTP reached 158611 kg whereas 255 the predicted cumulative SD reached 129253 kg (underestimation of 19 %). Comparing SDWS and observed SD at the WTP attested a 60 % recovery rate from WaterSed outputs which was consistent with previous hydrogeological tracing results on connected sinkholes. Hanin (2011) suggested the existence of fast karstic connections between sinkholes and springs on Radicatel catchment, with a 62 % recovery rate. The recovery rate dropped to 50 % for the predicted SD. For the 126 events for which erosion and runoff occurred on the catchment, the cumulative amount for both observed and predicted SD was 260 estimated to 88114 kg and 78870 kg respectively, suggesting that direct transfer accounts for 55-61 %. These results were https://doi.org/10.5194/hess-2020-363 Preprint. Discussion started: 22 September 2020 c Author(s) 2020. CC BY 4.0 License. consistent with previous published results of Masséi et al. (2003) in same karstic environment (Norville catchment, Normandy, France), which evaluated the part of direct transfer to 73 % during erosive events.

265
The monthly-backward chaining nested cross-validation procedure made it possible to assess the generalization capacity of the DNN on a new data set. This procedure was repeated twelve times on all the data, while keeping at least one full hydrologic year as input for the modelling. The modelling results were less efficient than for the complete dataset but overall satisfactory (Fig. 5). The median NSE values for the training and the test sets was above the lower limit of the threshold (Moriasi et al., 2007;Ritter et al., 2013). The median RMSE values on the training and the test sets was in the same order of magnitude as for 270 the complete dataset (300 -500 kg). Despite honorable performances, the difficulties of the model to generalize suggested a classical problem in machine learning, the bias-variance dilemma (Geman et al., 1992). Ideally, the model should accurately reflect irregularities in the training data while generalizing to data testing, but it is not possible to do both at the same time.
The higher NSE and lower RMSE values on training sets suggest a complex model with high variance and a reliable bias that represents the training data fairly well but present a learning risk on the test data. The model can thus represent part of the 275 random noise of the learning dataset.

Scenarios Predictions 280
Once the calibration was completed and the architecture of the DNN defined, we applied the model for all predefined scenarios.
As a first step, and using the WaterSed model, we quantified SDWS and RWS at the connected sinkholes for the three additional scenarios and the five DSP, comparing them to the baseline scenarios ( Fig. 6a-b). For the first scenario (S_grass), 33 % of grasslands were ploughed up, which led to an increase of the spatial extent of runoff generation on the catchment. SDS_grass ranged from 8468 to 188264 kg for P0.5 to P100, which led to an average increase of 4.74 % compared to the results on the 285 baseline scenario. The SD increase was maximal (+7.45 %) for the 10-years DSP. The increase in runoff was higher, with an average increase of 8.4 % on the five DSP, ranging from 2211 to 67097 m 3 . The second scenario (S_engi) considered the implementation of erosion control measures by 2050 (i.e. 181 fascines and 13.1 ha of grass strips), which led to a global decrease in sediment fluxes rates. Simulated SDS_engi ranged from 1659 to 133458 kg to connected sinkholes, resulting in an average decrease of 44 %. This scenario was more effective on small return periods (i.e. 0.5 and 2-years), with a SD decrease 290 between 79 and 55 % respectively. Erosion control measures were less effective on the three other DSP, with a decrease ranging from 25 to 34 %. The effects on runoff was not significant, with an average decrease of 7 %. These results were not surprising considering that this management plan only included 13.1 ha of grass strips, and the main interest of the fascines lies in their ability to trap suspended sediments, especially the coarsest elements (AREAS, 2012). The last scenario (S_farm) considered the adoption of good farming practices (i.e. +15 % infiltration capacity) on 50 % of the plots. It can be observed 295 that this scenario significantly reduced sediment production on the hillslopes (Fig. 7a-b). The simulated SDS_farm ranged from 741 to 134001 kg for an average decrease of 49 % compared to the baseline scenario results. The simulated values on water runoff ranged from 720 to 47442 m 3 . This scenario was more effective in reducing SD and runoff, in average, than the ecoengineering scenario. The SD decrease induced by good farming practices ranged from 25 to 90 % at the connected sinkholes.   DSP characteristics and predicted SD/runoff by the WaterSed model at connected sinkholes, and for all scenarios, were injected as inputs in the DNN model. On the baseline scenario, predicted SD at the WTP ranged from 576 to 12200 kg (see section 4.2), and the Fig. 8 shows the global trends of the scenarios modelling results. For the S_grass scenario, predicted SD at the 310 WTP were slightly higher than the baseline scenario, reaching 12894 kg for the 100-years DSP. This scenario resulted in an average increase of 4.5 % for the SD at the WTP. The SD increase was less significant on small return periods (+0.42-3.29 %) than longer return periods (+5.68-6.82 %). On the third scenario (S_engi), the modelling results suggested a significant SD decrease at the WTP. The predicted SD ranged from 515 to 9810 kg that led to an average decrease of 25.4 %. The SD reduction was particularly effective on the 2-years DSP reaching a decrease of 44.7 %. The smallest SD reduction was observed for the 315 0.5 years DSP (10 %) and ranged from 19 to 30 % on longer return periods. These results were consistent with those found by Fournier et al. (2008), who showed by a parametric interpolation model that erosion control measures at connected sinkholes on the Caumont aquifer led to a 36 % decrease of the turbidity in the Varras WTP (Normandy, France). The last scenario https://doi.org/10.5194/hess-2020-363 Preprint. Discussion started: 22 September 2020 c Author(s) 2020. CC BY 4.0 License.
(S_farm) was the more efficient with a predicted SD ranging from 223 to 9165 kg and a high average decrease of 43 %. The SD reduction ranged from 28.9 to 61.3 % and was the more efficient for the 0.5 and 2-years return period (55.6-61.3 %). 320

Coupled modelling approach: strength and uncertainties
The original modelling approach developed in this study allowed the coupling of surface and subsurface sediment transfer processes. As karstic processes are complex and difficult to model due to the lack of knowledge of their geometry, the coupling of a process-based GIS model and a data-driven model (i.e. DNN) proved to be a powerful tool. This approach makes it easy to assess the impact of different land-use scenarios on drinking water quality at drinking water treatment plants. Despite very 330 encouraging results, one may noticed that the coupled model was trained on only two years of complete selected measurements, and therefore may not have captured the full distribution of possible cumulative daily rainfall on the catchment and/or turbidity values at the drinking water treatment plant. To ensure that the model may not have a weakness on the representativeness of extremes events, we performed an additional evaluation using the Generalized Extreme Value Distribution (GEV) that was largely applied to extreme events such as rainfall or river discharges (Carreau et al., 2013;De Michele and Avanzi, 2018). We 335 applied the GEV distribution to SD time series at the WTP of Radicatel. The maximum annual SD observed at the WTP was extracted for all complete hydrologic years (i.e. annual maximal blocks; n = 22) and the frequency distribution was assessed using the package 'extRemes' under R software (R Development Core Team, 2008;Gilleland and Katz, 2016; see Fig. C1).
We compared the SD simulated at the WTP by the DNN to the distribution values predicted by the GEV. Firstly, we simulated overall SD for the baseline scenario at the connected sinkholes and for the five DSP, using the WaterSed model. Then, we 340 used these modelling results as inputs for the DNN model and applied it to the five DSP. Secondly, we compared the results with the calculated GEV distribution (Fig. 9). The SD predicted at the WTP by the DNN model was in good agreement with https://doi.org/10.5194/hess-2020-363 Preprint. Discussion started: 22 September 2020 c Author(s) 2020. CC BY 4.0 License. the distribution calculated by GEV even with a high dispersion for highest return period. Even if it is well known that deep learning-based methods may results in weak performance for extreme events (Zhang et al., 2019), the results obtained here give confidence in the model's ability to simulate extreme events thanks to a careful selection of input data. The degree of 345 confidence in the model's output could be improved with longer time series. One limitation of the presented coupled approach is the data availability and its ability to be validated at other study sites, as the turbidity data backup are not always properly done in all WTP. From now on, the drinking water supplier should keep turbidity records to allow the application of these data-driven models, with the guarantee of a quality control to avoid interruptions in the time series.

Implications for Future Land Use Strategies
The modelling results in this study suggest that two different land use strategies (i.e. increase of capacity infiltration and eco-355 engineering) can significantly reduce the SD incoming to connected sinkholes during extreme rainfall events (up to 90 %), and therefore, decrease the SD at the WTP (from 10 to 61%). The first strategy leads to a decrease in sediment production through simplified cultural techniques, while the second affects transfer processes. The adoption of better farming practices (e.g increase of crop cover, no-tillage, reduced tillage) inducing an increase of 15 % of the capacity infiltration on 50 % of the agricultural plots on the catchment, appears to be slightly more efficient than the implementation of eco-engineering 360 infrastructures (181 fascines + 13.1 ha of grass strips) for reducing SD at the WTP. These results thus show that it is more interesting to adopt a land use strategy aimed at reducing hydro-sedimentary production directly on cultivated plots than transfer processes on the catchment, but these strategies could be combined. Additional simulations integrating both effects of best farming practices on cultivated plots and the implementation of eco-engineering infrastructures show a tenfold effect. As illustrated in Fig. 10, the coupling of the two strategies makes it possible to reduce significantly the SD at the WTP. The 365 combined effect will not add up, but we can expect an improved decrease in sediment discharge at the WTP from 40 % to 80 https://doi.org/10.5194/hess-2020-363 Preprint. Discussion started: 22 September 2020 c Author(s) 2020. CC BY 4.0 License.
%. Both land use strategies can also provide interesting effects on biodiversity and ecosystem services (Posthumus et al., 2015) and/or can improve soil resilience and promote sustainable agriculture (Lal, 2014). Public policies leading to the implementation of erosion control measures can be economically viable and efficient (Patault et al., 2020). Even if some simplified cultural techniques may imply negative returns for farmers, they may be eligible for agri-environment payments 370 (Posthumus et al., 2015). The implementation of these land use strategies may also require specific maintenance to keep their initial performance (Frankl et al., 2018) or specific training and machines for farmers. According to the simulations, the ploughing up of 33 % of grasslands for the benefit of agricultural plots on the catchment by 2050, will not increase significantly the SD at the WTP (less than 5 % in average). Our simulations just extended the current trend observed in the studied region. Despite all, some precautions must be taken regarding the results of this scenario which could be higher or lower, depending on the localization of the ploughed-up plots. Same observations were made by Souchère 380 et al. (2005)

Conclusions
In this study, a new coupled modelling approach was developed in order to help decision-makers choosing an adapted erosion 390 and runoff management strategy to reduce the impact of sediment discharge on drinking water supply. Expert-based GIS model (WaterSed) was used to simulate erosion and runoff at connected sinkholes on the Radicatel catchment (Normandy, France) and was coupled to a data-driven model (i.e. DNN) to simulate karstic transfer. This new approach can be easily implemented, does not require a knowledge of the geometry of the karstic system studied, and demonstrated the values of understanding hillslopes erosion and runoff processes to model hydro-sedimentary in karstic systems. Our modelling results suggest that the 395 coupled model performed well during the calibrating phase. The coupled model was able to generalize on unknown datasets through an adapted monthly backward chaining nested cross-validation procedure, but the results are just slightly over the threshold limit and suggest that more research should be done for the application on other catchments. The results also suggest that a land use change scenario considering the adoption of simplified cultural techniques can significantly reduce sediment discharge at the water treatment plant (up to 61%). This scenario also performed the one who considered the implementation 400 of eco-engineering control measures reducing erosion and runoff on the catchment (up to 45% at the water treatment plant).
The coupling of the two cited land use strategies has proven to be the most effective since it operates both on the hydrosedimentary production and transfer processes (decreasing SD at WTP from 40 to 80 %). Finally, ploughing up actual grasslands on the catchment will not increase significantly sediment discharge at the water treatment plant. In the framework of this study, we can consider that the results might be different considering a different spatial organization of grasslands on 405 the catchment and we suggest conducting a specific study on this hypothesis. https://doi.org/10.5194/hess-2020-363 Preprint. Discussion started: 22 September 2020 c Author(s) 2020. CC BY 4.0 License.