Modelling daily river discharge using stochastic differential equation in Ouémé at Savè Basin, Benin, West Africa

.


INTRODUCTION
There is much evidence nowadays that the planet is warming up, largely as a result of human generated greenhouse gases (Kundzewicz et al., 2014;IPCC 2014a, b). Extreme hydrological phenomena such as floods and droughts are caused by the acceleration of climate variability. Benin, like all West African countries, has been experiencing increasing climate fluctuations since 1970 (Sintondji et al., 2014). Thus, all sectors depending on water availability (such as agriculture, water supply, hydroelectricity, etc) are now highly vulnerable to the effects of climate change. This situation has prompted hydrologists to question the existing rainfall-runoff modelling methods, and several modelling attempts aim to contribute more effectively to the reduction of uncertainties associated with the assessment of the impacts of climate variability and changes on water resources. As a result, we require an approach that can account for both the dynamics and the random nature of the physical phenomenon.
Mathematically, stochastic systems can be modelled by stochastic differential equation (SDE) in the time continuous case (Zorzano et al., 1999). In the following, we will restrict our considerations to SDE with Gaussian white noise. In many problems, Gaussian white noise has been used as an approximation of the stochastic term. (Pauluhn, 1993). The solution of this SDE is the Markovian diffusion process. One of the most important advantages of the SDE is the associated Fokker-Planck equation (FPE), which allows one to directly derive the time-varying probabilities associated with the outflow (Biao et al., 2016). The FPE models the time evolution of the probability distribution in a system under uncertainty. The major problem for the FPE is that this equation does not have an explicit solution. It can be solved numerically or by approximation using orthogonal polynomials. Alamou (2011) and Biao et al. (2016) used the Hermite polynomial expansion to approximate the exact solution of the FPE. To the best of our knowledge, no study has so far investigated the numerical resolution of the FPE in our study area, even in West Africa; none have investigated the time evolution of the density probability function of the river discharges. The present study fills this gap.
The objective of this paper is, therefore, to consider the uncertainties in the rainfall-runoff modelling. The deterministic hydrological model based on the least action principle (ModHyPMA) is compared with its stochastic model in terms of simulations of discharges. Finite differences and finites volumes methods are among the other methods that have been widely used to numerically solve the FPE. A successful approach to overcoming the limitations of simple finite differences is achieved by Wojtkiewiczet et al. (1999) in terms of higher order finite differences. Thus, the time-dependent probability distribution for the resulting discharge is obtained in the form of the numerical methods of finite differences of higher order (order 4) and finite volumes.

Study area and data used
The Ouémé at Savè catchment stretches out from the center to the North of Benin between 7° 58-10° 12 N and 1° 35-3° 05 E (Figure Alamou et al. 111 1). It covers an area of 23,600 km 2 (Le Barbe et al., 1993) which is about 47.2% of the whole Ouémé catchment. On a global view, Benin extends from the Niger River to the Atlantic Ocean, with a relatively flat terrain, small mountains (about 600 m), and low coastal plains with marshlands, lakes, and lagoons. The study area has a unimodal rainfall season (from mid March to October) that peaks in August. The interannual mean rainfall is around 1075 mm; the minimum is 680 mm (in 1983), and the maximum is 1693 mm (in 1963). The Ouémé catchment landscape is characterized by gallery forest, savannah, woodlands, agricultural lands, pastures and mosaics of cropland and bush fallow, plantation with parkia, cashew, and palm trees (Bossa, 2007). In the Ouémé catchment, the main landscape elements are the crest and the upper, middle, and lower slopes (together referred to high peneplains), followed by the valley fringe and colluvial footslopes, and the valley bottoms and terraces (low peneplains and floodplains) (Igué et al, 2000). Data used in this study consist of daily rainfall data, daily potential evapotranspiration (ETP), and daily discharge data. Rainfall and ETP data were provided by Meteo-Benin, while discharge data were provided by the National Directorate of Water (DG-Eau Benin). The period 1961 -2010 has been chosen as the study period (good compromise, taking into account the length of the data available in the different stations). Spatialized regional daily mean rainfall was obtained by kriging (Matheron, 1970) with an exponential variogram.

ModHyPMA Model description
The hydrological model based on the least action principle (ModHyPMA) is described by Afouda et al. (2004), Afouda and Alamou (2010), Alamou (2011) and Biao et al. (2016) as follows: (1) (2) Where: Z describes the variation of the initial state of the catchment, ψ describes the model input, Y describes the model structure, λ and μ are the physical parameters of the model, and Q stands for the river discharge.
Equation 1 describes explicitly the production process (that is, the action of the unsaturated zone that accounts for evaporation and evapotranspiration, and divides the resulting rainfall event into two components: overland and underground) and Equation 2 describes the transformation process (that is, the process by which the rainfall volumes for the overland component and the underground component are transformed into runoff) (Alamou, 2011). The scheme and main equations of the deterministic ModHyPMA are given in Figure 2.

Stochastic differential equation describing the River Basin and its associated Fokker-Planck equation
To account for the different types of uncertainties (that is, uncertainties related to our imperfect knowledge of the physical phenomenon and uncertainties related to the quantitative evaluation of the parameters of the environment, particularly rainfall), the stochastic formulation of ModHyPMA is derived from its foundation in the SDE, which has already been successfully applied in a wide range of hydrological applications. From Equation 2 we can derive Equation 3 as : ( 3) Let us focus on Equation 2, is the variation of discharge Q with respect to time. In the case of daily discharge, this variation can be written in the form of ∆ which represents the variation of discharge from day t to day + ∆t, with ∆t = 1day. Let us now model the uncertainties by a Gaussian white noise process which can be added to the structure of the deterministic model Y. The stochastic process represents the variation which is therefore given by = ̂+ with ̂= − 2 −1 + 1 , . We will now continue the reasoning with the process, which is also a Gaussian white noise, with r̅ = [ ] and σ r representing the intensity of the noise. For simplicity, let us set = r̅ and G = σ r . Equation 4 reads, By replacing the expression of ̂ and r t , becomes: The question which arises now is to find the values of the constants R and G. To this end, we can deduce that, This process can be seen as the residuals from the approximation of by the deterministic model ModHyPMA.
Therefore, can be obtained from the observed data using: Thus, the estimators of R and G are given in the form: Let us set, Thus, Equation12 can be written in the form, The Gaussian white noise ε t can be described as the formal derivative in time of a Wiener process, W(t), i.e. dW(t) = ε(t)dt. Therefore, the SDE that describes the River Basin is given by Equation 14.
The FPE that is associated with the SDE describing the River Basin is (Risken, 1989): (15) with the initial conditions,

Numerical approximation of the solution of the stochastic differential equation (SDE): Euler Scheme
The numerical Euler Scheme to the SDE given in Equation 14) is: with +∆ − = ∆~ 0, 1 and ∆ = 1 .
The scheme and the main equations of our stochastic approach are summarized in Figure 3.

Model performance criteria
To evaluate the model performance for calibration, the following criteria were taken into account: the Nash-Sutcliffe Efficiency (NSE) (Nash and Sutcliff, 1970), the coefficient of determination (R 2 ), the root mean square error (RMSE) and the mean absolute error (MAE). The Nash-Sutcliffe Efficiency is defined as: where , , and ̅ stand, respectively for the observed discharge, simulated discharge, and the arithmetic mean of for all events, i = 1 to n. The NSE can attain values from -1 to 1. The value of 1 indicates the total agreement of observed and simulated discharges. An efficiency of 0 (NSE = 0) indicates that the model predictions are as accurate as the mean of observed data; whereas an efficiency less than zero (NSE < 0) occurs when the observed mean is a better prediction than the model; in other words, when the residual variance is larger than the data variance. Essentially, the closer the model is to 1, the more accurate the model is. The coefficient of determination is a number that indicates how well a data fits a statistical model, sometimes simply a line or curve. The coefficient of determination ranges from 0 to 1. An R 2 of 1 indicates that the regression line perfectly fits the data. The RMSE represents the average distance between the simulated and observed data, Numerical solution of the Fokker-Planck equation using the finite difference of high order 4 and finite volume methods

Finite difference method
In numerical analysis, finite-difference methods are a class of numerical techniques for solving differential equations by approximating derivatives with finite differences. Today, finitedifference methods are one of the most common approaches to the numerical solution of partial derivative equations, along with finite element methods (Grossmann et al., 2007). The higher the order of the finite difference method, the greater the accuracy of the numerical solution. This is the reason for chosing the finite differences of order 4 for solving the FPE. The procedure used in solving boundary problems by finite difference method is as follows: ii. Numerical scheme of the finite differences of the partial derivative equation: We used the finite difference scheme of order 4 for the first and second derivatives of P(Q, t) with respect to Q and the backward scheme of order 1 for the first derivative P(Q, t) with respect to t. The choice of the backward scheme for the first derivative of P(Q, t) with respect to t is due to the fact that it leads to an implicit scheme that is always stable.
iii. Finite difference equation at mesh points: At each fixed date ≥ 1, the index i must be varied between 1 and m -1 in order to write the finite difference equations at m -1 mesh points. We have to solve, for each given j, a system of m -1 equations with m -1 unknowns.
iv. System of discrete algebraic equations: Before finding the = 1 − ( ) =1 = 0 + ∆ ∈ ℕ, ∆ = 1 = 0 + ∈ ℕ discrete solutions to the equation at any date j, we have to find all the solutions up to j -1 since we need information up to this order. Therefore, the numerical solution of our problem must begin with the resolution of the system obtained at j = 1. The system of equations is given by: The system to be solved for j = 1 is therefore { 1 }: 1 1 = 0 . The system to be solved at date j is therefore the following: v. Solving the system of equations: At each date j, there is a system of equations (Sj) to be solved. For the resolution of the system, direct Gaussian method has been applied.

Finite volume method
The finite volume method is used to represent and evaluate partial differential equations in the form of algebraic equations (LeVeque, 2002). In the finite volume method, volume integrals in a partial differential equation that contains a divergence term are converted to surface integrals using the divergence theorem. Finite volume methods can be compared and contrasted with the finite difference methods, which approximate derivatives using nodal values, or finite element methods, which create local approximations of a solution using local data, and a global approximation constructed by stitching them together. In contrast, a finite volume method evaluates exact expressions for the average value of the solution over some volume and uses these data to construct approximations of the solution within cells (Fallah et al., 2000). At a date t = n, the numerical solution of the problem is obtained by solving the following: It is therefore necessary to start from the numerical resolution of the Alamou et al. 115 system obtained at n = 1 to find the discrete solutions of the equation at any date n.

Simulation of discharge with deterministic ModHyPMA
The hydrological model has been calibrated over the period 2003 -2007 and validated over the period 2009 -2010. Figure 4 shows the result of the simulated hydrograph compared with the observed discharge for the calibration period. Figure 5 presents the same data for the validation period. The difference between the observed and simulated results can be seen by a simple visual control, and the numerical values for NSE, R 2 , RMSE and MAE as presented in Table 1. The recession curve is quite well simulated. However, the uncertainties associated with the peaks are greater than those associated with low flow. These findings are in line with Biao et al. (2016). The fact that the discharge peaks are not well simulated can be attributed to data errors (Andréassian et al., 2010;Kuczera et al., 2010). The improper representation of uncertainty is an intrinsic drawback of the deterministic hydrological models since they do not include components that enable the preservation of the associated statistical characteristics of the observed data. For both the calibration and the validation periods, the NSE and R 2 are greater than 0.70 (Table 1), while RMSE of 123 and MAE of 51 for the calibration period and RMSE of 184 and MAE of 90 for the validation period were achieved. These results indicate that ModHyPMA is suitable for the simulation of river discharge in the Ouémé River Basin. However, there are still many sources of uncertainty not being taken into account by ModHyPMA.

Stochastic model
The simulation with the stochastic model has been performed over the period 2009-2010. Figure 6 shows in the same graph the observed discharge, the simulated discharge with deterministic ModHyPMA, and the simulated discharge with the stochastic model. It can be seen from this figure that the simulated peaks with the stochastic model are much close to the peaks from observed data than the simulated peaks with the deterministic model. The NSE and R 2 are greater than 0.89, while RMSE of 113 and MAE of 76 were achieved.

Numerical solutions of the Fokker-Planck equation
The numerical solutions of the Fokker-Planck equation are the distributions probability of the river discharge. These distributions express the lack of confidence (that     numerical methods: the finite differences of order 4 and the finite volume methods. The minimum value of discharge Q is equal to 0 m 3 /s and the maximum value Qmax = 5000 m 3 /s, a value that cannot be exceeded by the discharge.
Let us start first with m = 10,000 meshes to compare the two methods over 100 days. Figures 7, 8 and 9 show, respectively the numerical solutions of the two methods for days 0 + 2, 0 + 20 and 0 + 100. We have gone to the order of 10 -4 before it can be noticed that the density probability of the discharge from the finite difference method is hidden by the density probability obtained with the finite volumes approach. From Figure 10, we notice a difference of at most 2.595×10 -4 -2.59×10 -4 = 5 ×10 -7 between the two peaks. The above results lead us to conclude that the two solutions are convergent for 10,000 meshes. We now reduce the number of meshes to 1000 at 0 + 100 days, and then to 100 to compare the two numerical methods. The results obtained are presented in Figures 11 and 12.
From Figures 11 and 12, the slight difference between the two solutions appears from m = 100 meshes. All  these results are satisfactory and allow us to conclude that the two numerical solutions are convergent. We can therefore confirm that the derived density probability of discharge is close to the exact solutions of the FPE. Let us now continue with only the finite difference method and we consider 5,000 meshes. Figure 13 shows the density probability of discharge at some dates ( 0 + 2, 0 + 10, 0 + 50 0 + 100 days). It can be seen that the density probability of discharge flattens out and covers more discharges as time evolves. This can be justified by the fact that there are more uncertainties when moving from 0 0 + 10. The further Figure 10. Comparison between the density probability of discharge derived from the finite difference and finite volumes methods at 0 + 100 days. away from 0 , the more there will be uncertainties in the discharge concerning Q 0 . Indeed, the solution of the FPE is the density of probability knowing Q (t 0 ) = Q 0 : this is a transition density probability. We also notice that the probability density is symmetric like a normal distribution when it covers an interval that does not contain the boundaries (as in the case 0 + 2); this is due to the Gaussian white noise hypothesis we made at the beginning. However, the densities probability become asymmetric when they cover an interval containing one of  the limits; this is due to the limit conditions. Figure 14 shows the evolution of density probability of discharge in the function of time and discharge. This confirms also the findings described in Figure 14. The above results are also in line with the findings of Biao et al. (2016) who used the Hermite orthogonal polynomial to approximate the exact solution of the FPE.
Using the initial discharge 0 = 553.73 3 / at date 22/07/2011 ( 0 ), one can derive the mean and the confidence interval (CI) at 90% around the mean discharge at day 0 + 10 ( Figure 15). We can derive that: CI = [144 and 637 m 3 /s] and the mean = 410 m 3 /s. So, there is 90% chance that the observed discharge of the day 0 + 10 is between 144 and 637 m 3 /s. The observed  discharge at the date 0 + 10 is indeed within the confidence interval, close to the mean discharge. We can also calculate the probability that the discharge Q belongs to a given interval. For example 300 < < 550 = 0.602. This means that there is a 60.2% chance that the discharge of the day belongs to the interval (300 and 550 m 3 /s).

Conclusion
The main contribution of this paper was to use a stochastic approach to better account for both the dynamic and stochastic character of the hydrological phenomenon. The achievement of the objective of this paper stemmed from the combination of two modelling approaches: first, deterministic modelling of the hydrological system by using ModHyPMA; and second, the stochastic formulation of ModHyPMA in terms of a stochastic differential equation (SDE). In comparison to the deterministic model, we discovered that the stochastic model improves simulations of discharge in the Ouémé at Savè Basin. The resolution of the Fokker Planck equation that is associated with the SDE was done using two different numerical methods: finite differences and finite volumes. The convergence of the two numerical methods investigated allows us to provide a solution that is close to the exact solution of the Fokker-Planck equation. This numerical solution can be used to calculate the uncertainty in the discharges in the Ouémé at Savè Basin. Although the use of the SDE and the associated FPE as proposed in this paper may become more complicated, the potential benefits in decision making, data collection, and information value are promising.