A Novel Interannual Rainfall Runoff Equation Derived from Ol’Dekop’s Model Using Artificial Neural Networks

In water resources management, modeling water balance factors is necessary to control dams, agriculture, irrigation, and also to provide water supply for drinking and industries. Generally, conceptual and physical models present challenges to find more hydro-climatic parameters, which show good performance in the assessment of runoff in different climatic regions. Accordingly, a dynamic and reliable model is proposed to estimate inter-annual rainfall-runoff in five climatic regions of northern Algeria. This is a new improvement of Ol’Dekop’s equation, which models the residual values obtained between real and predicted data using artificial neuron networks (ANNs), namely by ANN1 and ANN2 sub-models. In this work, a set of climatic and geographical variables, obtained from 16 basins, which are inter-annual rainfall (IAR), watershed area (S), and watercourse (WC), were used as input data in the first model. Further, the ANN1 output results and De Martonne index (I) were classified, and were then processed by ANN2 to further increase reliability, and make the model more dynamic and unaffected by the climatic characteristic of the area. The final model proved the best performance in the entire region compared to a set of parametric and non-parametric water balance models used in this study, where the R2Adj obtained from each test gave values between 0.9103 and 0.9923.


Introduction
Precipitations are the origin of water resources, which undergo different quantitative and qualitative transformations on the slopes. The losses of rainwater are almost observed as a form of infiltration, retention in the soil, and evaporation. Some part of this quantity can also flow into wadis until it moves into the sea. The estimation of actual evapotranspiration is the component most required for estimating water and energy balance equations [1,2]. This resource can be accessed on the inter-annual scale, according to the history of climatic and hydrometric measurements that could be given in the outlet of some watersheds. For ungauged watersheds, the lack of some information posed a big problem in the estimation of the mean annual flow, which is the main scientific challenge for many hydrologists [3][4][5]. The water balance concept was considered to study the hydrological behavior of watersheds and to describe the relationship between water and thermal components of the earth. This relativity was defined by a mathematical ratio between rainfall (R), rainfall-runoff (RR), and real evapotranspiration (Ea) [6]. The estimation of rainfall-runoff is necessary as the first step to search for the best evaluation of Ea. In literature, the first attempts were started by Schreiber [7] and Ol'Dekop [8]. Then, Budyko [9] proposed an average model of both previous equations to minimize the estimation errors that were given by Schreiber over 16 watersheds numbered from 01 to 17, except basin 13 which represents the southern part of the country ( Figure 1A). In this region, the spatial precipitation distribution is very heterogeneous, characterized by a strong gradient from north to south, and a weak one from east to west [25]. North Algeria has a typically Mediterranean climate, which is hot and dry in summer. Contrariwise, in winter, it is mild and rainy. Annual precipitation data varies between 400 mm and 1500 mm. In the highlands and the Saharan Atlas the climate is generally semi-arid, and the rainfall does not exceed 500 mm annually. The pre-Saharan and Saharan regions have a very arid climate, which is almost devoid of precipitation, and the mean annual rainfall there varies between 50 mm and 200 mm [26,27]. The temperature of the country varies between day and night, as well as between summer and winter [28]. According to the data, the mean daily temperature is observed in January with a value between around 10 • C and 12 • C. In July it is comprised between 25 • C and 27 • C. In the highlands and the Saharan Atlas, the temperatures recorded a thermal amplitude higher than that of the coastal regions, given by an average daily temperature in January which is around 2 • C and 9 • C. However, in July these are between 19 • C and 33 • C. Algeria area, which includes a surface of 480,000 km 2 [24]. It is bordered in the north by the Mediterranean Sea, while the southern part overlooks the Great Sahara. It is located between a longitude of −2.2° and 8.6°, and a latitude of 33° and 37°, which is distributed over 16 watersheds numbered from 01 to 17, except basin 13 which represents the southern part of the country ( Figure 1A). In this region, the spatial precipitation distribution is very heterogeneous, characterized by a strong gradient from north to south, and a weak one from east to west [25]. North Algeria has a typically Mediterranean climate, which is hot and dry in summer. Contrariwise, in winter, it is mild and rainy. Annual precipitation data varies between 400 mm and 1500 mm. In the highlands and the Saharan Atlas the climate is generally semi-arid, and the rainfall does not exceed 500 mm annually. The pre-Saharan and Saharan regions have a very arid climate, which is almost devoid of precipitation, and the mean annual rainfall there varies between 50 mm and 200 mm [26,27]. The temperature of the country varies between day and night, as well as between summer and winter [28]. According to the data, the mean daily temperature is observed in January with a value between around 10 °C and 12 °C. In July it is comprised between 25 °C and 27 °C. In the highlands and the Saharan Atlas, the temperatures recorded a thermal amplitude higher than that of the coastal regions, given by an average daily temperature in January which is around 2 °C and 9 °C. However, in July these are between 19 °C and 33 °C. Spatial variation maps displaying the maximum and the minimum altitudes observed in 16 watersheds in northern Algeria are shown in Figure 1B,C, respectively. The maps show that the northern Algeria morphology is higher on the southern side than on the northern, where the minimum altitudes observed on the marine side reach measurements between 17.9 m and 271 m in the basins numbered 3, 2, and 4. However, the maximum altitudes in this region are between 262.2 m and 1079 m. In the southeastern region, the The hydrometric network of northern Algeria consists of 200 measuring stations, which are distributed over 16 watersheds. The present work aims to estimate the runoff of 102 sub-basins obtained by defining the extreme outlets of the wadis in the entire region. In the middle location of each sub-basin, we chose a meteorological and a hydrological station in the same coordinate to provide sufficient data that were used to assess rainfall flow between 1965 and 2020. The climatic data used in this study are precipitation (IAR) and temperature (IAT). The data was obtained from the available dataset on the climatic stations provided by the National Centers for Environmental Information (NCEI-NOAA). https: //doiwww.ncdc.noaa.gov/ (accessed on 12 May 2021), and Climate Knowledge Portal https://climateknowledgeportal.worldbank.org/ (accessed on 17 May 2021), whereas the hydrological and geographical data, which are the watershed area (S) and the length of the watercourse (WC), were obtained from hydrometric stations provided by the Algerian National Agency of Hydrological Resources (A.N.R.H).
The De Martonne aridity index (I) was used in this study as input data of the dynamic model given by the artificial neuron network (ANN 2 ) and to classify the bioclimatic floor of the northern Algeria area. The I index values were obtained according to the following Equation (1) [29]: The actual measurement of inter-annual rainfall-runoff (IARR R ) was obtained from the A.N.R.H service. This data was used to estimate the computational residues (IARR') obtained by Ol'Dekop's model and to control the performance of the proposed model in the study area [21], that is: where IARR E represents the inter-annual rainfall-runoff data estimated by the Ol'Dekop model, IAEa is the mean annual real evapotranspiration, and IARR R is the inter-annual real rainfall-runoff.

Artificial Neuron Network
An artificial neuron network (ANN) is a data-driven process with a flexible mathematical algorithm capable of solving the complex nonlinear relationships between input and output datasets. In fact, it mimics the biological neuron architecture [30]. It is a family of parallel architectures used to solve the most complex mathematical problems in modeling, optimization, and prediction [31,32]. In practice, using ANN involves taking into consideration three main elements:

•
The interconnection between the input data ensures good results through the process.

•
The transfer function controls the generation of the neural output.

•
The summing function and the statistical parameters describe how the weights of input data are adjusted during the treatment [30].
During computation the ANN receives data from the input layer, then a combination between selected data is performed by the hidden layer using the summing function and a number of statistical control parameters. In general, the summing function formula is represented by: where net j is the mean of weighted input for the jth neuron; W ij is the weight from the ith neuron in the previous layer to the jth neuron in the current layer, and X i is the input from the ith to the jth neuron.
The transfer function Ψ is also used by the hidden layer to generate the final result, which is called by output data Y j , using the function (5). The training stops when the error obtained by the validation test reaches the minimum.
where net j represents the obtained data from the summing function (4), θ j is the external threshold that is obtained from the summing step. In literature, hydro-climatic modeling using ANN typically uses the feed-forward network structure with either one or multiple layers, depending on the objective of the study [33,34]. When using ANN for modeling climatic and hydrological phenomena, a set of fundamental decisions need to be made [33], which are: • The choice to use the appropriate neural network architecture.

•
The best selection of a suitable training algorithm, study periods, and network structure.

•
The way to pre-process and post-process the input and output data, respectively.

ANN Transfer Function
In the experimental part, the transfer function used by the ANN model to estimate the output data was given by the trend equation of the multiple linear regression model (MLR). This latter one was applied to a set of input data (X i *) that could prove a good regression with IARR'*. The test was conducted after linearizing the data series using the Ln equation. In this machine learning, the output data represents the predicted residual values obtained by the Ol'Dekop model. The concept of multiple linear regression used to study the linear relationship between the dependent variable Y and the vector of regressors (X, X 2 , . . . , X k ) is given by the following function [35]: where α is the intercept, β is the slope or the coefficient, k is the number of observations and ε represents the estimation error.

Water Balance Model
A set of nonparametric and parametric models used to compare and analyze the performance of the proposed model are detailed by the Equations (7)-(11), as detailed next.

Schreiber
Schreiber [7] proposed a simple model to estimate inter-annual evapotranspiration (IAEa) in terms of inter-annual precipitation (IAR) and mean annual potential evapotranspiration (IAEo), that is:

Yang
Yang [37] proposed an alternative model (10) to assess mean annual actual evapotranspiration by entering local parameter (n) characteristics on Budyko's hypothesis. The model uses the watershed characteristics for a better estimation, that is where n > 0.

Sharif
This model is an improvement of the Mezentsev-Choudhury-Yang (MCY) proposition, achieved by replacing b, k, and n parameters with 0, 2, and 1, respectively [38].
where by IAR is the inter-annual rainfall; IAE 0 is inter-annual potential evapotranspiration; and IAE a is inter-annual actual evapotranspiration.

Metrics of Performance
The statistical criteria applied in each modeling step, which are also used to compare the performance of the obtained model with non-parametric and parametric models mostly cited in the literature, are well detailed in Table 1. Where by Qs is the estimated runoff; Qo is the observed runoff; N is the total number of ordinates; K is the number of independent variables; and e i is the residual for the time period i. Table 1. Statistical criteria used to evaluate the performance of sub-models proposed in all climatic regions.

Criteria Statistical Formula Reference
Coefficient of determination (R 2 ) Mean squared error (MSE) Root mean square error (RMSE) Mean absolute error (MAE)

Results
In this section, we start our study by spatially analyzing and classifying the input data used by each sub-model. Then we show the ANN design, such as the choice of input data by the selection criteria, the regression equation used as a transfer function, and the output results, which are mainly used to estimate the residual data (IARR'). Finally, the

Data Description and Classification
Data distribution and variability analyses of selected variables used for modeling IARR' are shown in this section, which aims to study the behavior and the regression relationship between input variables and the IARR' dataset. This study is based mainly on qualitative and quantitative tests, where the results of the data distribution are given by P-P plot, Q-Q plot (which allowed us to compare the empirical and theoretical distribution of data and cumulative quartiles, respectively, using normal law), and scattergram graph ( Figure 2). The Q-Q plot and the scattergram graph are obtained after linearizing the used variables by applying the Ln function, to simplify the comparison between data series distributions, which have different measurement units. Moreover, Table 2 gives us a set of statistical parameters to quantify the variability of the data series. The P-P plot and the Q-Q plot show that IAR, S, and WC data series have a similar variability to the IARR's dataset. However, IAR and S data distributions are closer to IARR than the WC data series ( Figure 2). On the other hand, IAT shows a different behavior of data distribution, where the graphs prove that the stochastic model, which defines the IAT variability, is closer to the normal law. The p-value results obtained by the D KS fit test using a set of distribution laws ( Table 2) show that IARR' and IAR follow the Weibull 3 law, where the p-values equal 0.9789 and 0.8572, respectively. On contrary, IAT data have a GEV distribution, where the p-value is equal to 0.8764. Moreover, S and WC show a similar distribution, in which both variables follow the law of gamma 2. According to the fit results, the last three variables accept the Weibull 3 as a second closest fit model, where p-values equal 0.6136, 0.8363, and 0.8263, respectively. In Figure 2, the scattergram shows a descriptive comparison between the variability of data cited above. The graphs show a similar variability between IARR', IAR, and S data series, where the majority of the values of each variable are very close between them and below the average of their series (which equals 49.113, 494.6529, and 719.7106, respectively). On the other hand, Table 2 shows that 25% of the dataset, which is bounded between the third quartile and the maximum value, has very large variability, given by the interval of [64.3229, 284.3221], [607.00, 1107.00], and [1028. 25,4050.00], respectively. The WC variable shows a slight difference between data distribution. Inversely, the variability is more similar to IARR'. In this series, the mean and the median are close, and equal 45.60 and 51.10, respectively. However, 25% of the WC values give a very high variability, which is given by a value range of [65.5250, 179.80]. On the other hand, the IAT dataset shows a very low variability given by a variation coefficient of 0.1180. This last dataset proves a convergence between the median and the mean, which are equal to 15.2792 and 15.4569, respectively. Moreover, the scattergram shows that the upper and the lower IAT values, compared to the average, have a similar distance, given by the interval of [−2, 0] and [0, 2], respectively. In Table 2, this similarity is given by [11.5833, 15.2792] and [15.2792, 21.7583], respectively. According to this analysis, the greatest variability was obtained from IARR' and S datasets, whereby the variation coefficients equal 1.1039 and 1.0379, respectively.
The spatial inter-annual rainfall distribution and the De Martonne index obtained by applying Equation (1) are mapped using the data of 102 meteorological stations to determine the bioclimatic floor of each watershed of northern Algeria between 1965 and 2020. This study helps to provide the application areas used to control the performance of the model proposed in this section. Climate classification analysis is shown in Figure 3. A very large variability of rainfall from north to south is shown in Map (a), wherein the southern part, the IAR reaches up to 200 mm. However, in the north, the rainfall reaches values higher than 800 mm.   The spatial inter-annual rainfall distribution and the De Martonne index obtained by applying Equation (1) are mapped using the data of 102 meteorological stations to determine the bioclimatic floor of each watershed of northern Algeria between 1965 and 2020. This study helps to provide the application areas used to control the performance of the model proposed in this section. Climate classification analysis is shown in Figure 3. A very large variability of rainfall from north to south is shown in Map (a), wherein the southern part, the IAR reaches up to 200 mm. However, in the north, the rainfall reaches values higher than 800 mm. According to the bioclimatic classification given by [45], northern Algeria has a climatic diversity spread over five floors, from very humid to very dry ( Figure 3B). The figure shows that watershed 02, which is located in the northern part and overlooking the Mediterranean Sea has a very humid climate, characterized by IAR values varying between 700 mm and >800 mm ( Figure 3A). On the other hand, watersheds 09, 15, and 10 According to the bioclimatic classification given by [45], northern Algeria has a climatic diversity spread over five floors, from very humid to very dry ( Figure 3B). The figure shows that watershed 02, which is located in the northern part and overlooking the Mediterranean Sea has a very humid climate, characterized by IAR values varying between 700 mm and >800 mm ( Figure 3A). On the other hand, watersheds 09, 15, and 10 have the greatest climatic diversity, varying between very humid, humid, semi-humid, Mediterranean, and semi-dry. Where in this area, the rainfall varies between 400 mm and >800 mm. This diversity depends mainly on the geographical and geological characteristics of the region. In the northeastern part, watershed 03 is characterized by a very humid, humid, and semihumid climate, where the IAR varies between 500 mm and >800 mm. The catchment area 01, 14, and 12 have a climatic diversity which is between humid, semi-humid, Mediterranean, and semi-dry. In this area, the minimum rainfall was observed in watershed 14, which arrives at 300 mm. Differently, in watersheds 01 and 12, the minimum rainfall values reach up to 200 mm in the southern areas.

Proposed ANNs
The different steps followed to obtain the best ANNs model for estimating the computational errors of IARR' given by the Ol'Dekop model are represented in this section. This study allows us to develop a new form of water balance model based on a set of climatic and geomorphological variables, which can be applied in a different area, without being conditioned by the aridity state of the watershed. The statistical tests used in this study are: residual analysis curve, R 2 , R 2 Adj , MSE, RMSE, and Durbin-Watson (DW). These latter tests were used to analyze the performance of each transfer equation used by the ANN model, and also to quantify the reliability degree of the proposed model compared to a set of parametrical and non-parametrical water balance models, which are what is mostly used in the literature. Our model is classified into two steps, given as ANN 1 and ANN 2 (Figure 4), which show that initially a local model (ANN 1 ) was given to estimate IARR in all northern Algeria watersheds. Then an improvement was made to increase the reliability of the previous model in each basin climatic area and to make it more dynamic and applicable in different regions (ANN 2 ). Figure 4 shows that the two previous sub-models are the type of feed-forward network with the architecture of (3-2-1-1) and (10-5-1-1), respectively. In the first attempts of this modeling, we used IAR, S, and WC variables as input layers in the ANN 1 model. After that, we classified the estimation results that were obtained by this model in groups, according to each bioclimatic level. In this step, we also used the aridity index I of each hydrological station as input data in the ANN 2 model to determine, in each climatic area, its transfer function, which allows us to deduce the final rainfall-runoff model. Figure 4, shows that the intermediate nodes (denoted by IRR 1 , IRR 2 , IRR 3 , IRR 4 , IRR 5 , IRR 6 , and IRR 7, respectively) make it possible to apply sub-processing, using a summing and transfer function on the input data "output layer" to estimate the output data in each step of the ANN model. In our case, the summing function combines the input variables two by two to have the best modeling results. Moreover, the choice of combination between variables was justified by the results of the correlation test, which were applied to the linearized variables compared to the linearized IARR' (IARR'*) ( Table 3). The table shows the degree of correlation between a set of candidate variables that were cited in the previous section of this study and IARR'*, using two correlation forms. A direct correlation was found between the IARR'* and (S*, WC*, IAT*, IAR*, and I*), then an indirect correlation between IARR'* and (WC*, S*), by studying the relationship between IARR'* and (( IARR * WC * ), ( IARR * S * )) then between (( IARR * WC * ), ( IARR * S * )) and (Wc*, S*), respectively.
ANN model, and also to quantify the reliability degree of the proposed model compared to a set of parametrical and non-parametrical water balance models, which are what is mostly used in the literature. Our model is classified into two steps, given as ANN1 and ANN2 (Figure 4), which show that initially a local model (ANN1) was given to estimate IARR in all northern Algeria watersheds. Then an improvement was made to increase the reliability of the previous model in each basin climatic area and to make it more dynamic and applicable in different regions (ANN2). Step ANN1: Basic model of inter-annual rainfallrunoff (local model).

2nd
Step ANN2: Spatial model of interannual rainfall-runoff (dynamique model).  Inter-annual rainfall-runoff residual data of Ol'Dekop (IARR'), inter-annual rainfall (IAR), inter-annual temperature (IAT), watershed area (S), watercourse (WC), De martone index (I), * input variable linearized by the Ln function, as follow: X i * = Ln X i . Table 3 shows that IARR'* has the best correlation with IAR*, which equals 0.8513. It is also strongly correlated with climatic data obtained from the I* index. Contrariwise, the IAT* shows a weak correlation. However, the geo-hydrologic variables such as WC* and S* are slightly correlated with the response variable (IARR'*), which equals 0.5141 and 0.5023, respectively. The results highlight that the relation proposed by ( IARR * WC * ) and ( IARR * S * ) shows a very good correlation with IARR'*, given by 0.8124 and 0.7322, respectively. Furthermore, the new variables also show a strong correlation with WC* and S*, respectively, given by correlation coefficients of 0.6168 and 0.7923, respectively.

IARR Modeling
Describing different modeling steps used to propose a new equation of IARR estimation is based mainly on the computational error analysis of the Ol'Dekop model in a set of bioclimatic floors. A non-linear regression relationship between IARR' and selected input variables provided in (Table 3) are shown graphically in Figure 5 as the first step of this analysis. In this regard, a direct regression is applied between the response variable (IARR') and the input variables (IAR, I). Then, intermediate variables were used to express the indirect relationship between (S, Wc) and IARR'. The results show a similar regression for each pair of the dataset (IARR ', IAR) and (IARR', I), given by an R 2 which equals 0.8608 and 0.8134, respectively.   The graphs shown in Figure 5A,B prove some trends of values, which are shown in the range of the maximum value. Regression models cited above are defined by Equations (12) and (13): On the other hand, the figure also shows a good nonlinear regression between IARR' and IARR S , and also between IARR' and IARR Wc according to R 2 results, which equal 0.7501 and 0.8209, respectively, given by Figure 5C,D. The regression relationship between both input variables are defined by Equations (14) and (15), respectively: Moreover, both ratios ( IARR S ) and ( IARR Wc ) proved a good regression trend with S and WC variables, respectively. Where R 2 equals 0.8103 and 0.6314, respectively. Both statistical relationships between input and response variables are defined by Equations (16) and (17): In this study, we found that all the cases of regression cited in Figure 5 followed the power model trend. Figure 6 represents the linear regression graphs, which express the degree of correlation between the real values of IARR'* and the estimated values that were obtained by the IRR 1 and IRR 2 models. In this step, we proved the correlation's degree and the reliability of the obtained model. We have well explained the different multiple regression models, which are applied to the pair of variables (IAR*, S*) and (IAR*, WC*), using the intermediate variable ( IARR * S * ) and ( IARR * WC * ), respectively, which showed a good nonlinear regression with S and WC data, and also with IARR' (Figure 5). Moreover, Table 4 shows a set of statistical parameters relating to this modeling, which gives information about the reliability and the trend analysis of the sub-models that are noted by (A, B, C, and D), compared to the linearized real data (IARR'*). According to the results, we found that the regression model obtained from (IAR*, IARR * S * ) and (IAR*, IARR * WC * ) have a very good regression, proved by an (R 2 , R 2 adj ) results which equal (0.9301, 0.9286) and (0.9405, 0.9393), respectively. Moreover, the errors given by MSE, RMSE, and DW show that all models did not prove a large trend deviation compared to IARR'* real data. The computational steps of IRR 1 and IRR 2 models are well detailed by the following Equations (18) 1 linear model ((a 1 X 1 + a 2 X 2 + B)), 2 final equation represented by power model (B × X 1 a 1 × X 2 a 2 ), no data (-), panels represented in Figure 1A-E, output response (Y), input variable (X), determination coefficient (R 2 ), adjusted coefficient of determination (R 2 Adj ), mean squared error (MSE), root mean square error (RMSE), Durbin-Watson coefficient (DW).
To deduce the previous equations as functions of fundamental variables IAR, S, and WC, we start with Equation (18) by replacing ( IARR S ) with (S) using Equation (16). So we obtain: On the other hand, we use Equation (17) Figure 6C,D shows the reliability of the results obtained by the IRR 1 and IRR 2 models that are given by Equations (18) and (19), respectively. The corresponding graphs show a good fit of regression between the actual and the estimated values of IARR'*. This reliability was shown by R 2 , and R 2 Adj statistical parameters, which equal (0.9036, 0.9006) and (0.9105, 0.909), respectively (Table 4). We apply the Exp function in Equations (20) and (21), to obtain the IRR model used by ANN 1 .
We have, In this step, we found two reliable equations to estimate IARR'. Figure 6E shows that the best regression can be obtained as a function of the three variables (IAR, S, and WC), which is given by Equation (24).
According to Equations (22) and (23), the general model IRR is defined as follows: Table 4 shows that the IRR model gives a good estimation, proven by a set of statistical parameters. The results obtained by this model give an R 2 and R 2 Adj , which are equal to 0.9518 and 0.9508, respectively. Moreover, the computational errors given by MSE, RMSE, and DW parameters show values of 208.539, 11.4409, and 0.7094, respectively. The ANN 1 model shows that the use of geomorphological parameters increases the reliability of estimation when compared with the simple regression model obtained by IAR only ( Figure 5A).
In the second step of this study, we improved the proposed model, which is defined by Equation (24) to be more dynamic and applicable to each bioclimatic region. We applied multiple linear regression to aridity index series (I) and the predicted data obtained previously by the IRR model. This technique was applied separately to each bioclimatic stage in order to find the corresponding estimation model of each area. Table 5 shows all the statistical parameters relating to each local regression model. The results show that the model obtained in a very humid climate area gives a more reliable estimation compared to the performance models of other bioclimatic floors, where the R 2 and R 2 Adj given for the IRR 7 model equal 0.9072 and 0.9001, respectively ( Figure 4 and Table 5). The table shows that the obtained models have different reliability in each climate region, wherein the Mediterranean area, the R 2, and R 2 Adj were proven to perform well, and equal 0.7804 and 0.7647, respectively.  6 , and IRR 7 , which were obtained by modeling IARR' in semi-dry, Mediterranean, semi-humid, humid, and very humid climate floors, respectively, are defined by Equations (25) We apply the Exp function on Equations (25)- (29)to generate the final estimation model (IRR F ) in terms of IAR, S, WC, and I ( We apply the different models obtained in each bioclimatic area using all the datasets of northern Algeria to show the trend that can be caused by the static models. In this step, we want to propose a dynamic model, by eliminating all constants and finding a mathematical relationship with variables that can prove a good correlation. Table 6 shows statistic results, obtained from IRR 3 , IRR 4 , IRR 5 , IRR 6 , IRR 7 , and IRR F models. The table shows that all previous models cited above proved an R 2 and R 2 Adj greater than 0.80. Moreover, the error trends increase in these models, which are proven by MSE, RMSE, and DW parameters. On the other hand, the IRR F model has the best reliability, given by an R2, which equals 0.9841. This model proved a low tendency, where the MSE, RMSE, and DW equal 62.5948, 5.9117, and 0.5250, respectively. The final model (IRR F ) was obtained by applying the weighted average using the R 2 values that are obtained in Table 6 as weighting coefficients to estimate the a1 and a2 coefficients of this model. Where the IRR F equation is defined as follows:  (3), Mediterranean (4), Semi-Humid (5), Humid (6), Very Humid (7), final inter-annual rainfall-runoff residuals model (IRR F ), regional interannual rainfall-runoff residuals model (IRR), determination coefficient (R 2 ), adjusted coefficient of determination (R 2 Adj ), mean squared error (MSE), root mean square error (RMSE), and Durbin-Watson coefficient (DW). We have, We have also, Using the results obtained by Equations (36) and (37) in Equation (35), we find: To make the obtained Equation (38) more dynamic, the constant (Cte) is defined in terms of IAR data, which is given by Equation (42). The Cte must be obtained by each watershed to take into consideration the variability of each climate area. For this, we propose in Equation (39) the hypothesis that the predicted and real data of the Ol'dekop residuals (IARR') are almost equal.
We replace EIRR F in Equation (39) by the formula defined in Equation (38), by doing so, we obtain IARR ∼ = Cte × EIRR 0.27808 × I 1.2364 (40) Cte ∼ = IARR EIRR 0.27808 × I 1.2364 Figure 7 shows that Cte data has a very good correlation with IAR, where R 2 equals 0.7235. When we use the trend equation obtained from the regression model to represent the Cte variable as a function of IAR, we obtain: (41) 4. Discussion Figure 8 shows regression graphs obtained from real and predicted IARR data that are estimated by a set of non-parametric and parametric water balance models most cited in the literature review. These are Schreiber, Ol'Dekop, Budyko, Yang (n = 2), and Sharif, as given by Figure 8B-F, respectively. These models were applied to the dataset series used in this modeling to compare the results with the predicted IARR values obtained by the new model, given in Figure 8A. The results demonstrate that the best IARR estimation are obtained by the proposed model, which proved a very good match between the estimated and the actual data ( Figure 8A), where the model proved no trend and a very good regression, given by an R 2 , which equals 0.9924. On the other hand, Schreiber and Yang's models give similar results in which R 2 equals 0.9211 and 0.9334, respectively ( Figure  8B,E). Budyko's model also gives a good result, and performs less than the previous models, given by an R 2 equal to 0.9041. On the other hand, Ol'Dekop and Sharif's models proved a similar performance given by an R 2 equal to 0.8525 and 0.8739, respectively. When we use Equations (22) and (39) in Equation (36), we find:  (42) In Equation (43), we used Equations (8) and (42) to define the new form of the Ol'Dekope model used for the IARR estimation.
4. Discussion Figure 8 shows regression graphs obtained from real and predicted IARR data that are estimated by a set of non-parametric and parametric water balance models most cited in the literature review. These are Schreiber, Ol'Dekop, Budyko, Yang (n = 2), and Sharif, as given by Figure 8B-F, respectively. These models were applied to the dataset series used in this modeling to compare the results with the predicted IARR values obtained by the new model, given in Figure 8A. The results demonstrate that the best IARR estimation are obtained by the proposed model, which proved a very good match between the estimated and the actual data ( Figure 8A), where the model proved no trend and a very good regression, given by an R 2 , which equals 0.9924. On the other hand, Schreiber and Yang's models give similar results in which R 2 equals 0.9211 and 0.9334, respectively ( Figure 8B,E). Budyko's model also gives a good result, and performs less than the previous models, given by an R 2 equal to 0.9041. On the other hand, Ol'Dekop and Sharif's models proved a similar performance given by an R 2 equal to 0.8525 and 0.8739, respectively. Sensors 2022, 22, x FOR PEER REVIEW 21 of 28 Moreover, in Figure 8C,F, the predicted data obtained by both models show a high tendency. Trend analysis of the new equation compared to the models cited above is given in Figure 9 and Table 7, whereby the figure represents residual curves obtained between measured and predicted data, and the table shows a set of statistical parameters used to analyze the variability and the performance of each model. Figure 9 shows that the new proposition gives the lowest error and no trend of residuals was obtained in Figure 9A when compared to data processed by all runoff models in the 102 sub-basin. Furthermore, the parametric model given by Yang (n = 2) is also reliable and can be taken as the second Moreover, in Figure 8C,F, the predicted data obtained by both models show a high tendency. Trend analysis of the new equation compared to the models cited above is given in Figure 9 and Table 7, whereby the figure represents residual curves obtained between measured and predicted data, and the table shows a set of statistical parameters used to analyze the variability and the performance of each model. Figure 9 shows that the new proposition gives the lowest error and no trend of residuals was obtained in Figure 9A when compared to data processed by all runoff models in the 102 sub-basin. Furthermore, the parametric model given by Yang (n = 2) is also reliable and can be taken as the second choice to estimate IARR in the region that has a great climatic diversity-i.e., the northern Algeria area. choice to estimate IARR in the region that has a great climatic diversity-i.e., the northern Algeria area.  On the other hand, Sharif's model shows good results compared to Schreiber, Ol'dekop, and Budyko, where the mean residuals obtained by this last one are symmetrical, distributed between negative and positive value ranges for the axis (X = 0). In this study, the Ol'Dekop model shows a very big tendency for residual data compared to other models, reaching up to 300 mm. In addition, Schreiber and Budyko produce a significant trend of data compared to real observations in watersheds, ranked from 61 to 102, which are located in semi-humid, humid, and hyper-humid climate regions (Figure 3). Table 7 shows that the IARR data series, which is obtained by the proposed model has a similar variability to the real data, given by a coefficient of variation, which equals 9292.465 and 9797.2623, respectively. Moreover, the mean values of both series are 81.4977 and 81.1249, respectively.
The results prove that the variability of data obtained by Schreiber's model is closer to real data compared to the data obtained by other models, where the coefficient of variation given by this method equals 5772.4952. On the other hand, the statistical analysis of data distribution shows that Yang's model provides data closer to the real dataset compared to data obtained by Schreiber's model; whereby the first quartile, the median, the mean, and the third quartile results equal 31.2638, 46.9607, 116.9946, and 86.4257, respectively (Table 7). In this analysis, the new model shows the best performance compared to all models used in this study proven by R 2 , R 2 Adj , RMSE, and MAE parameters, which equal 0.9924, 0.9923, 8.5073, and 5.2053, respectively. Figure 10 shows a set of performance tests applied on data series obtained by all models used in this study in five climate regions of northern Algeria.
The results show that the new model has the best reliability and more efficiency in estimating IARR in all watersheds. Where R 2 and R 2 Adj show values between 0.9 and 0.99, which was proved to be the best performance in the five climatic regions. Furthermore, the model's performance increases according to the rainfall condition of the basin, leading it to perform more in humid regions than in the arid. The new model is even able to assess flows produced by little precipitation, which is proven in the figure by an R 2 Adj greater than 0.9 in the semi-arid region. On the contrary, the other models are regional and they can be reliable in more arid areas than in humid areas. The R 2 and the R 2 Adj show that Yang (n = 2) and Sharif models perform better than non-parametric models, such as Schreiber, Ol'Dekop, and Budyko models in the semi-humid, humid, and very humid areas. Contrariwise, the performance of these models is better than the Yang and Sharif models in the arid region. The histogram of RMSE and MAE shows that the proposed model has the lowest amount of errors, which are closer to the axis (X = 0) during all observations. However, the other models have more computational residuals, which increase more in humid areas. Ol'Dekop's model has the highest values of RMSE and MAE observed in semi-humid, humid, and very humid climate areas. The results deduced by this analysis show the power of the proposed model to estimate runoff in different regions of the world. The use of the S and Wc input variables and the calibration applied to the model parameters using the De Martonne series is proof of the dynamicity of the new model and its ability to work with the varying climatic and geomorphological conditions of the watershed, even if it is ungauged. This advantage was justified by the steady performance and lack of a significant trend residual estimation during all the processes carried out in the 102 sub-basins. The results show that the new model has the best reliability and more efficiency in estimating IARR in all watersheds. Where R 2 and R 2 Adj show values between 0.9 and 0.99, which was proved to be the best performance in the five climatic regions. Furthermore, the model's performance increases according to the rainfall condition of the basin, leading

Conclusions and Future Work
This research work aims to propose a new equation derived from Ol'Dekop's model to estimate the inter-annual rainfall runoff in a large region, which is characterized by climatic diversity. The proposed model was obtained by applying new ANNs to estimate the residual values (IARR') given when comparing actual and predicted runoff assessed by Ol'Dekop. In this regard, multiple regressions were used as transfer functions and applied to a set of climatic and geomorphological input variables such as IAR, Wc, S, and I. As the first step, the sub-model, namely ANN 1 , proposed to estimate IARR' in the entire region. Then an improvement was made by classifying output data in each climatic region using the De Martonne index (I), which was also used to calibrate the coefficients of the final model (IRR F ). The performance analysis shows that the use of the IRR sub equation which is obtained by ANN 1 as a function of IAR, S, and Wc proved good reliability (R 2 Adj = 0.9518) with a significant error (RMSE = 208.539 mm) when estimating IARR' in the entire region. Contrariwise, the IRR F sub equation given by ANN 2 shows an improvement in the estimation model performance, where the R 2 Adj and RMSE equal 0.9789 and 62.5948 mm, respectively. The new equation proposed in this model shows the best performance compared to all parametric and non-parametric water balance models used in this study. This performance is given by an R 2 Adj equal to 0.9923 obtained for the entire northern region of Algeria, with the absence of significant errors proven by residual analyses. Inversely, the compared models are less performed. They are set by the R 2 Adj test between 0.8525 and 0.9333, with a tendency of residuals observed in some watersheds that are located in humid and very humid regions. The comparative performance analysis, which was applied in five climatic regions, shows the best and steady performance of the proposed model in all areas, given by an R 2 Adj more than 0.92, with RMSEs less than 15 mm. However, the other models are regional and perform better in arid regions than in humid. A large trend of residuals is observed by non-parametric models such as Schreiber, Ol'Dekop, and Budyko in the very humid region, where the RMSEs are greater than 100 mm.
Our future work aims to further improve this model designed to estimate rainfallrunoff data in ungauged watersheds using different time scales. We also want to propose a new model using machine learning for estimating rainfall-runoff in sandy basins.
Author Contributions: All authors of this paper have directly participated in the planning, execution, and analysis of this study. O.M. worked on data collection, analysis, and mapping. A.A. worked on the ANN design, statistical modeling, and validation. A.L. worked in statistical modeling, investigation, and validation. K.M. participated in validation and methodology followed in the redaction. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.