## Abstract

In this study, Artificial Intelligence (AI) models along with ensemble techniques were employed for predicting suspended sediment load (SSL) via single station and multi-station scenarios. Feed Forward Neural Networks (FFNNs), Adaptive Neuro-Fuzzy Inference System (ANFIS), and Support Vector Regression (SVR) were the employed AI models, and simple averaging (SA), weighted averaging (WA), and neural averaging (NA) were the ensemble techniques developed for combining the outputs of the individual AI models to gain more accurate estimations of the SSL. For this purpose, twenty-year observed streamflow and SSL data of three gauging stations, located in the Missouri and Upper Mississippi regions were utilized at both daily and monthly scales. The obtained results of both scenarios indicated the supremacy of ensemble techniques to single AI models. The neural ensemble demonstrated more reliable performance compared to other ensemble techniques. For instance, in the first scenario, the ensemble technique increased the predicted results up to 20% in the verification phase of the daily and monthly modeling and up to 5 and 8%, respectively, in the verification step of the second scenario.

## HIGHLIGHTS

Three AI models were applied for SSL modeling of a river.

Ensemble of output of models was performed by combining three methods.

Modeling was done using single and multi-station strategies.

Ensemble technique increased the predicted results up to 20% in verification phase.

## INTRODUCTION

The whole amount of sediment output, carried from a watershed by rivers and streams, is described as sediment load (SL) and can be categorized into two groups of bed load (BL) and suspended sediment load (SSL). Fine sediments fall within SSL and could be moved long distances before deposition, even under base and low-flow conditions (Salih *et al.* 2019). SL prediction is a crucial task in hydro-environmental and water resources management as it can provide beneficial information about watershed erodibility and the deposition of sediment caused by the scouring phenomenon. Moreover, it affects water quality through: (i) creating visual pollution by increasing the turbidity of the water which also increases the drinking and domestic water's purifying expenses; (ii) carrying contaminants like pesticides, nutrients, and other chemicals; and (iii) reducing the dissolved oxygen in the water and consequently limiting the water life. It also has a significant effect on the design and maintenance of water constructions such as the dead volume of a dam, stable channels and river-bed equilibrium (Pandey *et al.* 2020). Moreover, one of the destructive phenomena caused by sediment transport is scouring and it is often dominant around bridge piers (Singh *et al.* 2020) and in culvert outlets (Najafzadeh & Kargar 2019; Pandey & Azamathulla 2021).

Considering all these impacts, the accurate prediction of SL is crucial and is a highly challenging task since it has a vast number of temporal and spatial variabilities. Owing to this, several approaches have been introduced to model and address this issue (Verstraeten & Poesen 2001; Ward *Et Al.* 2009) which are usually categorized into three main groups: physical-based, conceptual, and black-box models (Nourani & Mano 2007). Conceptual and physical models deal with equations and rules involved in the phenomenon and could be beneficial in gaining knowledge of the process but they have their limitations and complexities (Li 1979; Asselman 1999; Xu 2002; John *et al.* 2021; Pu *et al.* 2021). The main obstacle for physical-based and conceptual models is that these models require high-resolution and high-quality data of flow and sediment which are not often accessible (Sivakumar 2006) and in case of direct measurement of these data, aside from its high expense, any subtle error may affect the modeling results. This also could be time-consuming which encourages researchers to employ data-driven (black-box) models in cases where accurate forecasting is much more important than physical illumination. In recent decades, artificial intelligence (AI) methods have received great attention and been utilized in a wide spectrum of hydrologic processes due to their practical advantages and reliable results (ASCE 2000; Sahoo *et al.* 2006; Solgi *et al.* 2014). Since sedimentation and its transport process involve high sophistication, dynamism, and nonlinearity in both spatial and temporal scales, employing AI approaches often results in dependable outcomes (Jain 2001; Kisi 2004; Zhu *et al.* 2007; Kisi *et al.* 2008; Partal & Cigizoglu 2008; Nourani 2009a; Mirbagheri *et al.* 2010). Nagy *et al.* (2002) predicted suspended sediment concentration (SSC) employing an artificial neural network (ANN) via several stream data collected from reliable sources. The network was developed based on several input variables including Reynolds number, Froude number, Stream Width Ratio, and Mobility number, and the SSC as the output parameters. Comparing the ANN model and various commonly used sediment discharge formulas performed on observed data, the ANN demonstrated better results in model testing than the other commonly used sediment discharge formulas. Cigizoglu (2004) developed a multi-layer perceptron (MLP) network for estimation of daily SSL and the results showed that MLPs capture the complex non-linear behavior of the SSL series much better than the conventional models. Azamathulla *et al.* (2013) employed an extension of genetic programming (GP) called gene expression programming (GEP) along with an adaptive neuro-fuzzy inference system (ANFIS) and regression models to predict the SSL relation of three Malaysian rivers. The results indicated the outperformance of the GEP approach. Nourani (2014) reviewed AI-based models employed for modeling SSL and investigated the advantages and disadvantages of both AI-based and hybrid models like GP and indicated that an ANN is more reliable than GP. Nourani & Andalib (2015) used a wavelet-based least-squares support vector machine (WLSSVM) along with a wavelet-based ANN (WANN) to predict SSL at daily and monthly scales. Comparison of the results for both approaches revealed that WLSVM is more reliable and robust than WANN for estimation of sediment loads.

Considering the uncertainties involved in the sedimentation process on one hand and the ability of a fuzzy concept to handle the uncertainties on the other, different neuro-fuzzy (NF) and ANFIS models have been utilized for sediment load estimation. Rajaee *et al.* (2009) modeled SSL using ANNs, NF, Multi Linear Regression (MLR), and conventional sediment rating curve (SRC). Obtained results demonstrated that NF performs better than the others due to its nonlinear nature and its ability to handle uncertainties. To predict the functional relationship of sediment transport in sewer pipe systems, Azamathulla *et al.* (2012) used ANFIS and proved its robustness for practitioners. Support vector machine (SVM), as another kind of AI model, is an quite recently developed model that demonstrates reliable performance and can be utilized as an alternative method of ANN since SVM, instead of minimizing error, considers operational risk as the objective function and minimizes it. Nourani *et al.* (2016) performed spatio-temporal modeling of monthly SSL employing an SVM model and determined the non-linear relationships among the SSL for three stations located on the Ajichay River. Zounemat-Kermani *et al.* (2020) modeled both SSL and BL using conventional methods including SRC, MLR, and also AI models including ANFIS, support vector regression (SVR), and their integrated version with a genetic algorithm (GA-ANFIS, GA-SVR). The major findings of the research pointed out the outperformance of AI models over conventional models regardless of input combinations and the prediction results improved by utilizing integrated versions of the AI models.

Although AI-based models have been successfully applied in numerous studies to model SL and the results proved the efficiency of AI models, it is a common problem that employing different models to solve a unique issue may lead to distinct outcomes. Overall, as can be understood from the presented literature review of the applied AI models in SL prediction, there is almost no consensus among researchers to present a unique model as a leading one in estimation of SSL since these models are case sensitive and their performance may depend on the time of the employed data. For instance, one model may efficiently capture maximum values while another model may well present lower amounts (i.e. different performances of a single model for different seasons of the year) in a time series forecasting task. Hence, various features of the underlying pattern of a phenomenon could be obtained more exactly by compounding different models through the ensemble technique. The ensemble approach has already been discussed and applied in different fields of engineering and SL is not an exception (Shamseldin *et al.* 1997; Zhang 2003; Sharghi *et al.* 2018; Nourani *et al.* 2019, 2021; Sharafati *et al.* 2020). As an example, least-square ensemble models as well as an ANN model was used with different wavelet families linked to it (Wavelet-ANN models) to estimate the SL of the Toutle River in Washington state in which the ensemble models outperformed the other single models. However, as far as the authors are aware, the sediment process has not been investigated through the multi-station approach together with the ensemble concept. So the aim of this study was to utilize ensemble techniques to predict SSL in current time steps through single station and multi-station scenarios. For this purpose, FFNN, SVR, and ANFIS were employed to model the SSL of three stations located in the Missouri and Upper Mississippi regions at daily and monthly scales. Although there are numerous AI models to employ for modeling SSL, the motivations behind choosing these three models are that: FFNN is one of the fundamental models used in this field, SVR tries to minimize the operational risk and ANFIS employs fuzzy logic which is essential to handle the uncertainties of the modeling results. The reason for including a multi-station scenario in this study was that it could handle the uncertainty and non-linearity of SSL modeling via exploring the interaction of the stations’ information by recognizing the spatial and temporal variabilities, and consequently could provide more robust results. Moreover, the multi-station scenario made it possible to estimate the SSL of the downstream station using upstream data, and therefore measuring SSL values in downstream stations becomes unnecessary in the future when using this model. Afterwards, to enhance the prediction efficiency, the outputs of the mentioned AI models for each scenario were used as inputs to develop ensemble techniques. In this regard, two linear ensemble techniques, including simple averaging (SA) and weighted averaging (WA), and one nonlinear neural averaging (NA) were used for each scenario. So the objectives of this study can be summarized as: (i) AI-based modeling of SSL; (ii) improving the AI-based modeling via ensemble approach; and (iii) discovering multi-station relationships.

## MATERIALS AND METHODS

### Hydrometric stations and data sets

Three stations on the Missouri and Mississippi rivers located in three distinct states (i.e. Nebraska, Illinois, Missouri) were considered to carry out this research. The stations belong to three discrete basins and sub-basins of the Missouri and Upper Mississippi regions. The first station is located in Nebraska City and refers to the Missouri-Nishnabotna basin. The other two stations belong to the Upper Mississippi-Salt and Upper Mississippi-Meramec basins and are respectively located below Grafton after the confluence of the Mississippi and Illinois rivers, and in St. Louis after the confluence of the Mississippi and Missouri rivers.

To have an overview of all the characteristics of the basins, Figure 1 shows the spatial locations of the three stations and Table 1 presents some general information about the stations and land use statistics of the related basins.

. | Missouri River at Nebraska city . | Mississippi River Below Grafton . | Mississippi River at St. louis . |
---|---|---|---|

Denoted by | A | B | C |

Station USGS ID | 06807000 | 05587455 | 07010000 |

Hydrologic unit code (HUC) number | 10240001 | 07110009 | 07140101 |

State | Nebraska | Illinois | Missouri |

Hydro-region name | Missouri Region | Upper Mississippi Region | Upper Mississippi Region |

Sub-basin name | Keg-Weeping Water | Peruque-Piasa | Cahokia-Joachim |

Drainage area (km^{2}) | 1,061,896 | 443,665 | 1,805,223 |

Latitude | 40°40′55.01″N | 38°57′4.17″N | 38°37′44.40″N |

Longitude | 95°50'49.01″W | 90°22'16.41″W | 90°10'47.20″W |

State soil | Holdrege | Silty Clay Loam | Menfro |

Percentage of forest land | 8.34% | 17.58% | 11.41% |

Percentage of agricultural land | 24.71% | 60.67% | 39.31% |

Percentage of urban land | 2.46% | 8.38% | 4.55% |

. | Missouri River at Nebraska city . | Mississippi River Below Grafton . | Mississippi River at St. louis . |
---|---|---|---|

Denoted by | A | B | C |

Station USGS ID | 06807000 | 05587455 | 07010000 |

Hydrologic unit code (HUC) number | 10240001 | 07110009 | 07140101 |

State | Nebraska | Illinois | Missouri |

Hydro-region name | Missouri Region | Upper Mississippi Region | Upper Mississippi Region |

Sub-basin name | Keg-Weeping Water | Peruque-Piasa | Cahokia-Joachim |

Drainage area (km^{2}) | 1,061,896 | 443,665 | 1,805,223 |

Latitude | 40°40′55.01″N | 38°57′4.17″N | 38°37′44.40″N |

Longitude | 95°50'49.01″W | 90°22'16.41″W | 90°10'47.20″W |

State soil | Holdrege | Silty Clay Loam | Menfro |

Percentage of forest land | 8.34% | 17.58% | 11.41% |

Percentage of agricultural land | 24.71% | 60.67% | 39.31% |

Percentage of urban land | 2.46% | 8.38% | 4.55% |

According to Figure 1 and Table 1, stations A and B are respectively in the upstream of the Missouri and Mississippi rivers and station C is located in the downstream of the Mississippi River. Concerning the drainage area, station C has the largest drainage area with station B having the smallest. Although the state soils of all three stations are subjected to severe erosion, as demonstrated in Table 1, the major land uses of the related watershed of station B are agriculture and forest lands whereas station A has fewer green areas than the other two stations. This may cause the basin of station A to be subjected to severe soil erosion and consequently bring about a larger amount of SSL compared to station B. Accordingly, it is expected that station A has more interactions with station C due to its less suitable land cover and land use.

For better understanding of catchment flow discharge and SSL conditions, Figure 2 provides time series of SSL and discharge of station C on a daily scale and statistical analysis of the used data on a daily scale are also summarized in Table 2.

Station . | Statistical parameter . | All data . | Training data . | Verification data . | |||
---|---|---|---|---|---|---|---|

Q(/s) . | SSL(tons/d) . | Q(/s) . | SSL(tons/d) . | Q(/s) . | SSL(tons/d) . | ||

A | 1,194 | 55,760 | 1,196 | 57,263 | 1,186 | 51,252 | |

6,258 | 701,968 | 6,258 | 701,968 | 3,398 | 525,022 | ||

245 | 2,460 | 245 | 3,920 | 442 | 2,460 | ||

Standard deviation () | 675 | 68,835 | 746 | 71,285 | 395 | 60,689 | |

Coefficient of Variation^{a} (CV) | 0.56 | 1.23 | 0.62 | 1.24 | 0.33 | 1.18 | |

B | 3,870 | 63,877 | 3,632 | 61,129 | 4,582 | 72,122 | |

13,734 | 754,000 | 12,431 | 599,000 | 13,734 | 754,000 | ||

379 | 489 | 379 | 755 | 535 | 489 | ||

Standard deviation () | 2,572 | 83,168 | 2,400 | 73,504 | 2,918 | 106,653 | |

Coefficient of Variation^{a} (CV) | 0.66 | 1.30 | 0.66 | 1.20 | 0.63 | 1.47 | |

C | 6,340 | 216,127 | 6,173 | 214,353 | 6,838 | 221,451 | |

23,389 | 2,093,129 | 20,275 | 2,093,129 | 23,389 | 1,826,514 | ||

1,529 | 4,230 | 1,529 | 6,550 | 1,625 | 4,230 | ||

Standard deviation () | 3,871 | 296,763 | 3,777 | 292,785 | 4,100 | 308,410 | |

Coefficient of Variation^{a} (CV) | 0.61 | 1.37 | 0.61 | 1.36 | 0.60 | 1.39 |

Station . | Statistical parameter . | All data . | Training data . | Verification data . | |||
---|---|---|---|---|---|---|---|

Q(/s) . | SSL(tons/d) . | Q(/s) . | SSL(tons/d) . | Q(/s) . | SSL(tons/d) . | ||

A | 1,194 | 55,760 | 1,196 | 57,263 | 1,186 | 51,252 | |

6,258 | 701,968 | 6,258 | 701,968 | 3,398 | 525,022 | ||

245 | 2,460 | 245 | 3,920 | 442 | 2,460 | ||

Standard deviation () | 675 | 68,835 | 746 | 71,285 | 395 | 60,689 | |

Coefficient of Variation^{a} (CV) | 0.56 | 1.23 | 0.62 | 1.24 | 0.33 | 1.18 | |

B | 3,870 | 63,877 | 3,632 | 61,129 | 4,582 | 72,122 | |

13,734 | 754,000 | 12,431 | 599,000 | 13,734 | 754,000 | ||

379 | 489 | 379 | 755 | 535 | 489 | ||

Standard deviation () | 2,572 | 83,168 | 2,400 | 73,504 | 2,918 | 106,653 | |

Coefficient of Variation^{a} (CV) | 0.66 | 1.30 | 0.66 | 1.20 | 0.63 | 1.47 | |

C | 6,340 | 216,127 | 6,173 | 214,353 | 6,838 | 221,451 | |

23,389 | 2,093,129 | 20,275 | 2,093,129 | 23,389 | 1,826,514 | ||

1,529 | 4,230 | 1,529 | 6,550 | 1,625 | 4,230 | ||

Standard deviation () | 3,871 | 296,763 | 3,777 | 292,785 | 4,100 | 308,410 | |

Coefficient of Variation^{a} (CV) | 0.61 | 1.37 | 0.61 | 1.36 | 0.60 | 1.39 |

^{a}Dimensionless.

Checking the data set for trend and randomness and regarding the SSL and discharge time series of the used data set (Figure 2), it can be understood that flow and erosion conditions have not encountered any significant transformation from 1997 to 2017 and also there is no trend and randomness in the used data set as checked by the Mann-Kendall test.

As Table 2 reports, stations A and B have similar mean SSL values while the mean flow discharge of station B is almost three times bigger than the mean flow discharge of station A, indicating much sedimentation happening in the Missouri region. Not surprisingly, both flow discharge and SSL for station C are remarkably larger than for the other two stations since it is located after the confluence of the Missouri and Mississippi rivers. Furthermore, higher values of coefficient of variation (CV) for SSL compared to the flow discharge data show the higher deviation of SSL data from the mean and describe SSL as a more sporadic phenomenon than the flow rate series itself. It is worth noting that when some uncertainties of data are added into the dataset, the ANFIS model could handle some of the uncertainties due to its employment of fuzzy concepts and the ensemble technique could reduce the uncertainties of the applied models.

*X*and

*Y*can be written in the form of (Yang

*et al.*2000):where

*x*and

*y*are the probability distributions of variables

*X*and

*Y*; H(

*x*) and H(

*y*) show respectively the entropies of distributions

*x*and

*y*, and H(

*x*,

*y*) is their joint entropy as:where p

_{XY}(

*x*,

*y*) is the joint distribution.

In this study, to assess the possibility of applying a multi station scenario, MI was used to explore the correlation of the discharge and SSL of stations A and B and discharge of station C with the SSL time series of station C (see Table 3).

Time series . | . | . | . | . | . |
---|---|---|---|---|---|

1.19 | 1.27 | 1.50 | 1.32 | 1.52 |

Time series . | . | . | . | . | . |
---|---|---|---|---|---|

1.19 | 1.27 | 1.50 | 1.32 | 1.52 |

denotes the SSL of Station *i* and denotes the flow discharge of station *i*, all in current time step. As can be seen from Table 3 both stations A and B have strong interdependency with the SSL of station C. However, the discharge of station B has a slightly larger MI value which could be due to the closer distance of stations B and C.

### Proposed methodology

In this study, three different AI-based black-box models including ANN (a commonly used AI method), ANFIS (an AI method which deals with uncertainties involved in the process), and SVR were developed to predict the SSL of the study stations. Then, to improve accuracy and decrease the uncertainty of the predictive models, the outputs of these single models were fed into the ensemble unit comprised of three ensemble techniques including SA, WA, and NA methods. The modeling was carried out through two scenarios of single station and multi-station modeling.

- (i)
First scenario

*i*stands for a modeling station and SSL of the

*i*th station at time step

*t*is considered as a function (

*f*) of

*i*th station's SSL value at previous time steps (

*t*

*−*

*1, t*

*−*

*2, …*) up to lag time

*n*, and discharge value in time

*t*and previous time steps (

*t*

*−*

*1, t*

*−*

*2, …*) up to lag time

*m*. It is worth noting that the most efficient lags (n, m) and input could be determined through a trial and error procedure.

(ii) Second scenario

It is worth noting that, as can be seen in Figure 1, in addition to the main rivers of the study area (i.e. Missouri and Mississippi) there are several other tributaries that may have some effect on the SSL amount of the studied rivers. The effect of these tributaries have been included by developing black-box models which employ the historical data and contain related information about the upstream conditions.

### Feed forward neural network (FFNN)

Originally inspired by biological neural networks, an ANN is a family of statistical learning algorithms that can be used to estimate or approximate nonlinear functions with an arbitrary number of inputs. Feed-forward neural networks (FFNNs) with a back-propagation (BP) learning algorithm are well-known utilized strategies in dealing with various engineering issues. An FFNN consists of interconnected processing elements known as nodes with unique characteristics of information processing such as learning, nonlinearity, noise tolerance, and generalization capability. FFNN structure contains three layers, namely, the input, hidden, and output layers (see Figure 3). In an FFNN, the inputs presented to the input layer's neuron are propagated in a forward direction and a nonlinear function known as activation function is used to compute the output vector. It has been already shown that a FFNN trained by a BP algorithm with only three layers of input, hidden and target layers can be decently employed for different subjects of water resources engineering (ASCE 2000; Nourani 2017).

### Adaptive neural fuzzy inference system (ANFIS)

The ANFIS model includes the merits of both neural network and fuzzy control systems, by unifying ideas from fuzzy control and neural networks. The fuzzy structure has a strong inference system and has no learning ability while the neural network has powerful learning ability. ANFIS presents these two desired features in one single model. A fuzzy database, fuzzifier, and defuzzifier are the main components of every fuzzy system (Nourani & Komasi 2013). A fuzzy database is comprised of an inference engine and fuzzy rule base. As illustrated by Jang *et al.* (1997), the fuzzy rule base involves rules that are based on fuzzy propositions and fuzzy inference applies operation analysis. The ANFIS model architecture consists of five layers with layer 1 representing the input layer; layer 2 representing the input membership function (MF); layer 3 representing rules; layer 4 representing the output MFs; and layer 5 representing the output configured as illustrated in Figure 4.

*f*, for instance with two input vectors of

*x*and

*y*, the first order Sugeno inference engine may be applied to two fuzzy if–then rules as (Aqil

*et al.*2007):in which

*A*,

_{1}*A*and

_{2}*B*,

_{1}*B*show the MFs of inputs (i.e.,

_{2}*x*and

*y*),

*p*,

_{1}*q*,

_{1}*r*and

_{1}*p*,

_{2}*q*,

_{2}*r*are the target function parameters. The operation of an ANFIS can be briefly described as follows:

_{2}*i*th node in layer

*k*is denoted as

*Q*. Assuming a generalized bell function (gbellmf) as the MF, the output

^{k}_{i}*Q*can be computed as (Jang

^{1}_{i}*et al.*1997):where {

*a*,

_{i}*b*,

_{i}*c*} are adaptable premise variables.

_{i}*i*th rule towards the target is determined as (Jang

*et al.*1997):where, is the output of layer 3 and {

*p*,

_{i}*q*,

_{i}*r*} is the parameter set.

_{i}*et al.*1997):

To calibrate the premise parameters set {*a _{i}*,

*b*,

_{i}*c*} and consequent parameters set {

_{i}*p*,

_{i}*q*,

_{i}*r*} of the ANFIS, the conjunction of least squared and gradient descent methods is used as a hybrid calibration algorithm. (Aqil

_{i}*et al.*2007). Back propagation (Kurian

*et al.*2006) and hybrid learning algorithms (Ebtehaj & Bonakdari 2014) are two common algorithms widely used in various applications to train ANFIS and show somewhat accurate predictions. The Sugeno FIS was employed in this study among other fuzzy inference engines (Jang 1993).

### Support vector regression (SVR)

*x*is the input vector,

_{i}*d*is the actual value and

_{i}*N*is the total number of data patterns), the general SVR function is Wang

*et al.*(2013):where

*φ*(

*x*) indicates feature spaces, non-linearly mapped from input vector

_{i}*x*. Regression parameters of b and

*w*may be determined by assigning positive values for the slack parameters of

*ξ*and

*ξ** and by minimization of the objective function as:where is the weights vector norm and

*C*refers to the regularized constant specifying the exchange between the empirical error and the regularized term.

*ε*is called the tube size and is equivalent to the approximation accuracy placed within the training data points. The optimization problem can be changed to a dual quadratic optimization problem by defining Lagrange multipliers

*α*and

_{i}*α*. Vector

_{i}^{*}*w*in Equation (6) can be computed after solving the quadratic optimization problem as Wang

*et al.*(2013):

*n*stands for the number of data patterns,

*k (*

*,*

*)*is the kernel function performing the non-linear mapping into feature space and

*b*is the bias term. One commonly used kernel function is the Gaussian Radial Basis Function (RBF) as:where

*γ*is the kernel parameter. Figure 5 shows the architecture of SVM algorithms. For more details about SVR, readers are referred to Wang

*et al.*(2013) and Raghavendra & Deka (2014).

### Ensemble unit

In peak values modeling or in cases where different models have better results at various conditions and intervals, it is expected that the accuracy of a predicted time series will be improved by combining (ensemble) the outputs from several prediction models. In this regard, utilizing combined predictions is less risky than depending on an individual method. According to several types of research in various fields of engineering, it has been already confirmed that combining the results of different models could improve overall prediction accuracy (e.g. see Zhang & Berardi 2001; Kasiviswanathan *et al.* 2013; Nourani *et al.* 2020a).

In this study, two techniques were employed to combine the utilized models’ outcomes to enhance modeling results: (i) the linear ensemble technique, which includes SA (Equation (18)) and WA (Equation (19)); and (ii) the nonlinear ensemble technique (i.e. NA).

Unlike linear ensemble techniques, in the nonlinear ensemble technique a FFNN is trained which takes the outputs of single AI models as input and provides the ensemble output (Sharghi *et al.* 2019). In neural averaging (NA), the outputs acquired by individual models (i.e. FFNN, ANFIS, and SVR) are combined together to create a new model and train through the FFNN model to obtain an ensemble output (see Figure 6). It is worth mentioning that other AI models like ANFIS and SVR can be similarly employed for NA.

### Efficiency criteria

*n*,

*,*and are data number, observed data, averaged value of the observed data and calculated values, respectively. DC ranges between – ∞ and 1, with a perfect score of 1.

*m*is the number of input-output patterns and

*npar*is the number of parameters to be identified. While RMSE is expected to progressively improve as more parameters are added to the model, the BIC penalizes the model for having more parameters and therefore tends to result in more parsimonious models. The smaller the BIC value, the better the model performs.

## RESULTS AND DISCUSSION

### Results of the first scenario (single station modeling)

In scenario 1, all individual stations were modeled by their historical data using FFNN, ANFIS, SVR, and ensemble techniques. For this purpose, different lagged time series of streamflow and SSL were used to construct various input combinations to feed into the input layers of the models since streamflow and SSL demonstrate Markovian behavior, meaning that the parameters' values at the current time step could be related to their values at previous time steps. It is worth mentioning that since the soil type of each station is not a dynamic characteristic of the system in the single station scenario, the black-box models learn and then include this feature (i.e. soil type and land use) by self-learning. However, one of the major goals of this study is to observe the effect of different soil types and land use of different stations on the results, as the soil types and land uses differ from one station to the others in the multi-station scenario. To check the relations between each station's SSL, discharge and SSL of previous time steps, their MI values were computed and tabulated in Table 4.

Time series . | . | . | . |
---|---|---|---|

1.137 | 1.265 | 1.369 | |

1.098 | 1.225 | 1.328 | |

1.212 | 1.513 | 1.386 | |

1.201 | 1.508 | 1.385 | |

1.195 | 1.509 | 1.391 |

Time series . | . | . | . |
---|---|---|---|

1.137 | 1.265 | 1.369 | |

1.098 | 1.225 | 1.328 | |

1.212 | 1.513 | 1.386 | |

1.201 | 1.508 | 1.385 | |

1.195 | 1.509 | 1.391 |

As is clear from Table 4, the SSL time series of each station has a high relation with its previous SSL and discharge time series, but to determine which variables must be used together to lead to the best result, the trial and error procedure should be employed. Identifying and selecting the best input combination for the models is one of the crucial tasks of any modeling and as tabulated in Table 5, different combinations were developed and utilized to estimate the SSL value at the current time step at both daily and monthly time scales for the first scenario.

Modeling scale . | Combination No. . | Inputs . |
---|---|---|

Daily | 1 | , |

2 | , , , | |

3 | , , , | |

Monthly | 4 | , |

5 | ,, , | |

6 | , , , |

Modeling scale . | Combination No. . | Inputs . |
---|---|---|

Daily | 1 | , |

2 | , , , | |

3 | , , , | |

Monthly | 4 | , |

5 | ,, , | |

6 | , , , |

It is worth mentioning that, monthly scale (SSL value of the same month in the previous year) was also considered in the input layer to capture the seasonality of the process. To avert the overfitting problem in training AN FFNN, the network's proper architecture selection (i.e. number of the hidden layer's neurons and iteration epoch number) is an essential step in the modeling. Thus, the ranges of 10–100 and 1–10 were examined for epoch number and number of neurons of hidden layers, respectively. Accordingly, the FFNNs were trained by employing the tangent sigmoid as activation function of hidden and output layers using the gradient descent BP training algorithm (Haykin 1994) and the best structure specified via trial and error procedure for modeling the SSL of each individual station.

Due to its ability to handle the uncertainties and complex nonlinear behavior through fuzzy concepts, the ANFIS model was used in this study as another AI model. In training the ANFIS, the Sugeno fuzzy inference system was applied and Gaussian, Trapezoidal, and Triangular-shaped MFs exhibited better results than other MFs for modeling the hydroclimatic process (Nourani *et al.* 2020a). In addition, a constant MF was used in the output layers of the models. The number of training epochs was assessed along with the number of MFs in order to attain the finest ANFIS model. With this aim, the ranges of 5–100 and 2–4 were respectively examined to find the optimum number of epochs and MFs.

Afterward, SVR models were developed based on the RBF kernel which shows more reliable performance compared to other kernels (e.g. sigmoid and polynomial kernels) when considering smoothness assumption (Noori *et al.* 2011). By adjusting the RBF kernel parameters, optimal SVR was obtained for every input combination for each station. The attained results of the best structure for the AI models in daily and monthly scale are summarized in Table 6.

Modeling scale . | Station . | Inputs . | Model . | Structure^{a}
. | DC . | RMSE^{b}. | ||
---|---|---|---|---|---|---|---|---|

Calibration . | Verification . | Calibration . | Verification . | |||||

Daily | A | Comb. 1 | FFNN | (2-1-1, 20) | 0.86 | 0.79 | 0.04 | 0.04 |

ANFIS | Trapezoidal-2 | 0.90 | 0.82 | 0.03 | 0.03 | |||

SVR | 0.5, 0.1, 50 | 0.88 | 0.81 | 0.03 | 0.04 | |||

Comb. 2 | FFNN | (4-1-1, 40) | 0.88 | 0.82 | 0.03 | 0.04 | ||

ANFIS | Gaussian-2 | 0.91 | 0.80 | 0.02 | 0.03 | |||

SVR | 0.66, 0.1, 10 | 0.88 | 0.82 | 0.03 | 0.04 | |||

Comb. 3 | FFNN | (4-2-1, 10) | 0.89 | 0.84 | 0.03 | 0.03 | ||

ANFIS | Trapezoidal-2 | 0.91 | 0.85 | 0.03 | 0.03 | |||

SVR | 0.25, 0.1, 10 | 0.88 | 0.84 | 0.03 | 0.03 | |||

B | Comb. 1 | FFNN | (2-2-1, 30) | 0.90 | 0.70 | 0.03 | 0.07 | |

ANFIS | Trapezoidal-2 | 0.82 | 0.71 | 0.02 | 0.04 | |||

SVR | 0.5, 0.1, 8 | 0.87 | 0.73 | 0.03 | 0.05 | |||

Comb. 2 | FFNN | (4-2-1, 20) | 0.90 | 0.76 | 0.02 | 0.06 | ||

ANFIS | Double-Gaussian-2 | 0.93 | 0.62 | 0.02 | 0.08 | |||

SVR | 0.7, 0.1, 10 | 0.90 | 0.90 | 0.03 | 0.04 | |||

Comb. 3 | FFNN | (4-8-1, 40) | 0.92 | 0.79 | 0.02 | 0.06 | ||

ANFIS | Triangular-2 | 0.95 | 0.65 | 0.02 | 0.08 | |||

SVR | 0.5, 0.2, 15 | 0.88 | 0.82 | 0.03 | 0.05 | |||

C | Comb. 1 | FFNN | (2-2-1, 30) | 0.89 | 0.87 | 0.04 | 0.05 | |

ANFIS | Gaussian-2 | 0.90 | 0.89 | 0.04 | 0.04 | |||

SVR | 0.333, 0.1, 15 | 0.88 | 0.86 | 0.04 | 0.05 | |||

Comb. 2 | FFNN | (4-5-1, 40) | 0.94 | 0.87 | 0.03 | 0.05 | ||

ANFIS | Trapezoidal-2 | 0.93 | 0.90 | 0.03 | 0.04 | |||

SVR | 0.333, 0.1, 10 | 0.92 | 0.85 | 0.03 | 0.06 | |||

Comb. 3 | FFNN | (4-2-1, 20) | 0.92 | 0.91 | 0.04 | 0.04 | ||

ANFIS | Triangular-2 | 0.95 | 0.86 | 0.03 | 0.05 | |||

SVR | 0.25, 0.1, 14 | 0.90 | 0.87 | 0.04 | 0.05 | |||

Monthly | A | Comb. 4 | FFNN | (2-2-1, 30) | 0.74 | 0.61 | 0.08 | 0.08 |

ANFIS | Trapezoidal-2 | 0.67 | 0.63 | 0.09 | 0.07 | |||

SVR | 0.9, 0.1, 15 | 0.60 | 0.59 | 0.10 | 0.08 | |||

Comb. 5 | FFNN | (4-1-1, 40) | 0.35 | 0.30 | 0.14 | 0.12 | ||

ANFIS | Triangular-2 | 0.42 | 0.25 | 0.13 | 0.11 | |||

SVR | 0.25, 0.1, 10 | 0.34 | 0.25 | 0.14 | 0.12 | |||

Comb. 6 | FFNN | (4-5-1, 20) | 0.90 | 0.74 | 0.05 | 0.06 | ||

ANFIS | Triangular-2 | 0.85 | 0.72 | 0.06 | 0.07 | |||

SVR | 0.15, 0.1, 15 | 0.70 | 0.67 | 0.09 | 0.08 | |||

B | Comb. 4 | FFNN | (2-2-1, 50) | 0.77 | 0.52 | 0.06 | 0.14 | |

ANFIS | Triangular-2 | 0.70 | 0.52 | 0.07 | 0.14 | |||

SVR | 0.125, 0.1, 5 | 0.70 | 0.50 | 0.07 | 0.14 | |||

Comb. 5 | FFNN | (4-2-1, 70) | 0.44 | 0.35 | 0.09 | 0.16 | ||

ANFIS | Trapezoidal-2 | 0.50 | 0.20 | 0.09 | 0.18 | |||

SVR | 0.05, 0.1, 5 | 0.33 | 0.16 | 0.10 | 0.19 | |||

Comb. 6 | FFNN | (4-4-1, 30) | 0.89 | 0.61 | 0.04 | 0.12 | ||

ANFIS | Gaussian-2 | 0.83 | 0.60 | 0.05 | 0.11 | |||

SVR | 0.05, 0.1, 10 | 0.76 | 0.62 | 0.06 | 0.12 | |||

C | Comb. 4 | FFNN | (2-2-1, 40) | 0.81 | 0.60 | 0.08 | 0.13 | |

ANFIS | Triangular-2 | 0.75 | 0.50 | 0.09 | 0.13 | |||

SVR | 0.9, 0.1, 15 | 0.80 | 0.65 | 0.08 | 0.13 | |||

Comb. 5 | FFNN | (4-2-1, 40) | 0.51 | 0.24 | 0.14 | 0.20 | ||

ANFIS | Triangular-2 | 0.45 | 0.005 | 0.14 | 0.23 | |||

SVR | 0.125, 0.1, 10 | 0.40 | 0.18 | 0.16 | 0.21 | |||

Comb. 6 | FFNN | (4-3-1, 60) | 0.89 | 0.63 | 0.07 | 0.13 | ||

ANFIS | Triangular-2 | 0.87 | 0.46 | 0.07 | 0.17 | |||

SVR | 0.25, 0.1, 10 | 0.85 | 0.66 | 0.08 | 0.13 |

Modeling scale . | Station . | Inputs . | Model . | Structure^{a}
. | DC . | RMSE^{b}. | ||
---|---|---|---|---|---|---|---|---|

Calibration . | Verification . | Calibration . | Verification . | |||||

Daily | A | Comb. 1 | FFNN | (2-1-1, 20) | 0.86 | 0.79 | 0.04 | 0.04 |

ANFIS | Trapezoidal-2 | 0.90 | 0.82 | 0.03 | 0.03 | |||

SVR | 0.5, 0.1, 50 | 0.88 | 0.81 | 0.03 | 0.04 | |||

Comb. 2 | FFNN | (4-1-1, 40) | 0.88 | 0.82 | 0.03 | 0.04 | ||

ANFIS | Gaussian-2 | 0.91 | 0.80 | 0.02 | 0.03 | |||

SVR | 0.66, 0.1, 10 | 0.88 | 0.82 | 0.03 | 0.04 | |||

Comb. 3 | FFNN | (4-2-1, 10) | 0.89 | 0.84 | 0.03 | 0.03 | ||

ANFIS | Trapezoidal-2 | 0.91 | 0.85 | 0.03 | 0.03 | |||

SVR | 0.25, 0.1, 10 | 0.88 | 0.84 | 0.03 | 0.03 | |||

B | Comb. 1 | FFNN | (2-2-1, 30) | 0.90 | 0.70 | 0.03 | 0.07 | |

ANFIS | Trapezoidal-2 | 0.82 | 0.71 | 0.02 | 0.04 | |||

SVR | 0.5, 0.1, 8 | 0.87 | 0.73 | 0.03 | 0.05 | |||

Comb. 2 | FFNN | (4-2-1, 20) | 0.90 | 0.76 | 0.02 | 0.06 | ||

ANFIS | Double-Gaussian-2 | 0.93 | 0.62 | 0.02 | 0.08 | |||

SVR | 0.7, 0.1, 10 | 0.90 | 0.90 | 0.03 | 0.04 | |||

Comb. 3 | FFNN | (4-8-1, 40) | 0.92 | 0.79 | 0.02 | 0.06 | ||

ANFIS | Triangular-2 | 0.95 | 0.65 | 0.02 | 0.08 | |||

SVR | 0.5, 0.2, 15 | 0.88 | 0.82 | 0.03 | 0.05 | |||

C | Comb. 1 | FFNN | (2-2-1, 30) | 0.89 | 0.87 | 0.04 | 0.05 | |

ANFIS | Gaussian-2 | 0.90 | 0.89 | 0.04 | 0.04 | |||

SVR | 0.333, 0.1, 15 | 0.88 | 0.86 | 0.04 | 0.05 | |||

Comb. 2 | FFNN | (4-5-1, 40) | 0.94 | 0.87 | 0.03 | 0.05 | ||

ANFIS | Trapezoidal-2 | 0.93 | 0.90 | 0.03 | 0.04 | |||

SVR | 0.333, 0.1, 10 | 0.92 | 0.85 | 0.03 | 0.06 | |||

Comb. 3 | FFNN | (4-2-1, 20) | 0.92 | 0.91 | 0.04 | 0.04 | ||

ANFIS | Triangular-2 | 0.95 | 0.86 | 0.03 | 0.05 | |||

SVR | 0.25, 0.1, 14 | 0.90 | 0.87 | 0.04 | 0.05 | |||

Monthly | A | Comb. 4 | FFNN | (2-2-1, 30) | 0.74 | 0.61 | 0.08 | 0.08 |

ANFIS | Trapezoidal-2 | 0.67 | 0.63 | 0.09 | 0.07 | |||

SVR | 0.9, 0.1, 15 | 0.60 | 0.59 | 0.10 | 0.08 | |||

Comb. 5 | FFNN | (4-1-1, 40) | 0.35 | 0.30 | 0.14 | 0.12 | ||

ANFIS | Triangular-2 | 0.42 | 0.25 | 0.13 | 0.11 | |||

SVR | 0.25, 0.1, 10 | 0.34 | 0.25 | 0.14 | 0.12 | |||

Comb. 6 | FFNN | (4-5-1, 20) | 0.90 | 0.74 | 0.05 | 0.06 | ||

ANFIS | Triangular-2 | 0.85 | 0.72 | 0.06 | 0.07 | |||

SVR | 0.15, 0.1, 15 | 0.70 | 0.67 | 0.09 | 0.08 | |||

B | Comb. 4 | FFNN | (2-2-1, 50) | 0.77 | 0.52 | 0.06 | 0.14 | |

ANFIS | Triangular-2 | 0.70 | 0.52 | 0.07 | 0.14 | |||

SVR | 0.125, 0.1, 5 | 0.70 | 0.50 | 0.07 | 0.14 | |||

Comb. 5 | FFNN | (4-2-1, 70) | 0.44 | 0.35 | 0.09 | 0.16 | ||

ANFIS | Trapezoidal-2 | 0.50 | 0.20 | 0.09 | 0.18 | |||

SVR | 0.05, 0.1, 5 | 0.33 | 0.16 | 0.10 | 0.19 | |||

Comb. 6 | FFNN | (4-4-1, 30) | 0.89 | 0.61 | 0.04 | 0.12 | ||

ANFIS | Gaussian-2 | 0.83 | 0.60 | 0.05 | 0.11 | |||

SVR | 0.05, 0.1, 10 | 0.76 | 0.62 | 0.06 | 0.12 | |||

C | Comb. 4 | FFNN | (2-2-1, 40) | 0.81 | 0.60 | 0.08 | 0.13 | |

ANFIS | Triangular-2 | 0.75 | 0.50 | 0.09 | 0.13 | |||

SVR | 0.9, 0.1, 15 | 0.80 | 0.65 | 0.08 | 0.13 | |||

Comb. 5 | FFNN | (4-2-1, 40) | 0.51 | 0.24 | 0.14 | 0.20 | ||

ANFIS | Triangular-2 | 0.45 | 0.005 | 0.14 | 0.23 | |||

SVR | 0.125, 0.1, 10 | 0.40 | 0.18 | 0.16 | 0.21 | |||

Comb. 6 | FFNN | (4-3-1, 60) | 0.89 | 0.63 | 0.07 | 0.13 | ||

ANFIS | Triangular-2 | 0.87 | 0.46 | 0.07 | 0.17 | |||

SVR | 0.25, 0.1, 10 | 0.85 | 0.66 | 0.08 | 0.13 |

For FFNN, the format “a-b-c, d” presents the number of input layers-hidden layers-output layer neurons, and epochs. For ANFIS, the different membership functions are represented as “MF-a”, indicating the MF and number of membership functions. For SVR, “a, b, c” denote γ, ε and c, respectively.

^{b}RMSE has no dimension since all data are normalized.

It can be seen from Table 6 that Comb. 3, in all models of daily time scale, led to slightly better results than the other two combinations which could be related to the incorporation of streamflow value in the current time step () in the input layer. The information provided by Table 6 indicates the outperformance of Comb. 6 in all models of monthly time scale as it includes current month streamflow value ( and SSL of the same month in the previous year ( which enable the model to capture the seasonality characteristics of the sedimentation process as well as an autoregressive property. Moreover, the success of Comb. 6 demonstrates the strong correlation between , , and .

In single station modeling, as Table 6 shows, the overall outputs indicate the supremacy and robustness of daily scale modeling compared to monthly modeling owing to the fact that a large amount of data are involved in daily scale training. Various studies have been carried out for SL modeling using AI models and they demonstrate different performances. Ghani *et al.* (2011) predicted the SL of three Malaysian rivers employing FFNN and the findings indicated the superiority of FFNN compared to traditional methods. Azamathulla *et al.* (2010) used SVM for modeling SL and proved its superiority to conventional SL models. Buyukyildiz & Kumcu (2017) developed ANN, SVR, and ANFIS models to predict SSL employing different input combinations. The results specified the superiority of SVR over others. Choubin *et al.* (2018) used classification and regression tree (CART), MLP neural network, ANFIS and SVM models to model the SSL of Haraz watershed in northern Iran. The results of that study demonstrated the superiority of classification-based models (CART and SVM) compared to others.

Both ANFIS and ANN were employed by Kumar *et al.* (2019) to predict the current day runoff and SSL in the Godavari basin and the outcomes demonstrated the outperformance of ANFIS to ANN. Comparing the reviewed studies in the literature, it can be understood that the accuracy of the AI models may vary for distinct case studies and different time scales. This is due to the SSL data stochasticity of each considered catchment as well as the capacity of the constructed AI-based models to handle the non-stationarity and nonlinearity in the data sets (Salih *et al.* 2019).

Like the previous studies reported in the literature, AI models demonstrated different performances for different stations in this study. The obtained results reveals that FFNN showed slightly better results for most of the stations. ANFIS, despite its ability to handle the uncertainties, demonstrated robust performance for some stations but also showed poor results for some others which could be related to the large number of ANFIS rules and parameters involved in the modeling and the nature of each station's data. For instance, it could not properly model station B in daily and station C in monthly time scales yet in other stations it showed as reliable results as the other models. Moreover, SVR presented the lowest performance compared to other techniques, except for station C in monthly scale that could be related to the used RBF kernel which is not as sensitive as the kernel used in FFNN (i.e. sigmoid). So from the given explanation, it could be understood that it is not an easy task to choose one model which provides the best prediction results for all stations and in all temporal scales. Consequently, combining the outcomes of all three models seems to be practical and could solve the problem of best model selection as well as improving the prediction results.

For the next step of modeling via scenario 1, the three previously described ensemble techniques were employed to combine the outputs of the single AI models for the best input combination of each time scale. An ensemble technique has been employed successfully in numerous fields of hydrology as a model combination technique (Elkiran *et al.* 2018; Sharghi *et al.* 2018; Choubin *et al.* 2019; Nourani *et al.* 2020b). Only the verification dataset was used to obtain the variables of both WA and NA techniques and, as for the single FFNN model, tangent sigmoid was considered as the activation function for the hidden and output layers in the NA technique. Then the neural ensemble network was trained using a scaled conjugate gradient scheme of BP algorithm and the best structure and epoch number of the network were determined through the trial and error procedure. The obtained results of daily and monthly ensemble models are tabulated in Table 7. For better comparison, Figures 7–9 present the observed versus computed daily SSL time series (computed by the NA ensemble and single models) of the verification phase along with the verification step scatter plots of the NA ensemble technique for all stations at daily scale.

Modeling scale . | Station . | Ensemble model . | Model structure^{a}
. | DC . | RMSE^{b}. | ||
---|---|---|---|---|---|---|---|

Calibration . | Verification . | Calibration . | Verification . | ||||

Daily | A | SA | – | 0.92 | 0.88 | 0.03 | 0.03 |

WA | 0.332, 0.336, 0.332 | 0.92 | 0.88 | 0.03 | 0.03 | ||

NA | 3-5-1, 60 | 0.95 | 0.93 | 0.02 | 0.02 | ||

B | SA | – | 0.94 | 0.86 | 0.02 | 0.06 | |

WA | 0.350, 0.288, 0.362 | 0.94 | 0.87 | 0.02 | 0.05 | ||

NA | 3-6-1, 50 | 0.95 | 0.90 | 0.02 | 0.04 | ||

C | SA | – | 0.95 | 0.91 | 0.03 | 0.04 | |

WA | 0.345, 0.325, 0.330 | 0.94 | 0.91 | 0.03 | 0.04 | ||

NA | 3-5-1, 10 | 0.95 | 0.92 | 0.03 | 0.04 | ||

Monthly | A | SA | – | 0.87 | 0.79 | 0.06 | 0.06 |

WA | 0.351, 0.336, 0.313 | 0.87 | 0.80 | 0.06 | 0.06 | ||

NA | 3-6-1, 10 | 0.92 | 0.77 | 0.05 | 0.07 | ||

B | SA | – | 0.87 | 0.65 | 0.04 | 0.12 | |

WA | 0.333, 0.328, 0.339 | 0.87 | 0.65 | 0.04 | 0.12 | ||

NA | 3-7-1, 80 | 0.92 | 0.68 | 0.04 | 0.12 | ||

C | SA | – | 0.90 | 0.69 | 0.06 | 0.13 | |

WA | 0.36, 0.262, 0.378 | 0.90 | 0.69 | 0.06 | 0.13 | ||

NA | 3-2-1, 30 | 0.92 | 0.70 | 0.06 | 0.13 |

Modeling scale . | Station . | Ensemble model . | Model structure^{a}
. | DC . | RMSE^{b}. | ||
---|---|---|---|---|---|---|---|

Calibration . | Verification . | Calibration . | Verification . | ||||

Daily | A | SA | – | 0.92 | 0.88 | 0.03 | 0.03 |

WA | 0.332, 0.336, 0.332 | 0.92 | 0.88 | 0.03 | 0.03 | ||

NA | 3-5-1, 60 | 0.95 | 0.93 | 0.02 | 0.02 | ||

B | SA | – | 0.94 | 0.86 | 0.02 | 0.06 | |

WA | 0.350, 0.288, 0.362 | 0.94 | 0.87 | 0.02 | 0.05 | ||

NA | 3-6-1, 50 | 0.95 | 0.90 | 0.02 | 0.04 | ||

C | SA | – | 0.95 | 0.91 | 0.03 | 0.04 | |

WA | 0.345, 0.325, 0.330 | 0.94 | 0.91 | 0.03 | 0.04 | ||

NA | 3-5-1, 10 | 0.95 | 0.92 | 0.03 | 0.04 | ||

Monthly | A | SA | – | 0.87 | 0.79 | 0.06 | 0.06 |

WA | 0.351, 0.336, 0.313 | 0.87 | 0.80 | 0.06 | 0.06 | ||

NA | 3-6-1, 10 | 0.92 | 0.77 | 0.05 | 0.07 | ||

B | SA | – | 0.87 | 0.65 | 0.04 | 0.12 | |

WA | 0.333, 0.328, 0.339 | 0.87 | 0.65 | 0.04 | 0.12 | ||

NA | 3-7-1, 80 | 0.92 | 0.68 | 0.04 | 0.12 | ||

C | SA | – | 0.90 | 0.69 | 0.06 | 0.13 | |

WA | 0.36, 0.262, 0.378 | 0.90 | 0.69 | 0.06 | 0.13 | ||

NA | 3-2-1, 30 | 0.92 | 0.70 | 0.06 | 0.13 |

Comb. 3 and Comb. 6 are the selected input combinations of daily and monthly scales, respectively.

For WA, the format “a, b, c,” denotes the weights of the FFNN, ANFIS, SVR, models. For NA, the format “a-b-c, d” in represents the number of input layers-hidden layers-output layer neurons, and epochs.

^{b}The RMSE has no dimension since all data are normalized.

Modeling scale . | Station . | Inputs . | Model . | BIC . | |
---|---|---|---|---|---|

Calibration . | Verification . | ||||

Daily | A | Comb. 3 | FFNN | −19,065.48 | −6,305.35 |

ANFIS | −18,833.09 | −6,102.58 | |||

SVR | −18,962.19 | −6,215.23 | |||

NA | −21,171.07 | −6,948.10 | |||

B | Comb. 3 | FFNN | −20,973.12 | −4,769.30 | |

ANFIS | −21,119.43 | −4,371.66 | |||

SVR | −18,962.19 | −5,282.46 | |||

NA | −21,128.04 | −5,644.86 | |||

C | Comb. 3 | FFNN | −17,492.14 | −5,780.04 | |

ANFIS | −18,901.94 | −5,229.89 | |||

SVR | −17,388.86 | −5,282.46 | |||

NA | −18,953.59 | −5,682.41 | |||

Monthly | A | Comb. 3 | FFNN | −378.25 | −41.87 |

ANFIS | −340.24 | −28.54 | |||

SVR | −303.61 | −49.18 | |||

NA | −378.25 | −32.63 | |||

B | Comb. 3 | FFNN | −449.57 | −24.86 | |

ANFIS | −414.60 | −34.17 | |||

SVR | −376.59 | −24.86 | |||

NA | −392.45 | 20.18 | |||

C | Comb. 3 | FFNN | −380.00 | −44.62 | |

ANFIS | −312.49 | 24.70 | |||

SVR | −324.81 | −20.05 | |||

NA | −449.29 | −77.37 |

Modeling scale . | Station . | Inputs . | Model . | BIC . | |
---|---|---|---|---|---|

Calibration . | Verification . | ||||

Daily | A | Comb. 3 | FFNN | −19,065.48 | −6,305.35 |

ANFIS | −18,833.09 | −6,102.58 | |||

SVR | −18,962.19 | −6,215.23 | |||

NA | −21,171.07 | −6,948.10 | |||

B | Comb. 3 | FFNN | −20,973.12 | −4,769.30 | |

ANFIS | −21,119.43 | −4,371.66 | |||

SVR | −18,962.19 | −5,282.46 | |||

NA | −21,128.04 | −5,644.86 | |||

C | Comb. 3 | FFNN | −17,492.14 | −5,780.04 | |

ANFIS | −18,901.94 | −5,229.89 | |||

SVR | −17,388.86 | −5,282.46 | |||

NA | −18,953.59 | −5,682.41 | |||

Monthly | A | Comb. 3 | FFNN | −378.25 | −41.87 |

ANFIS | −340.24 | −28.54 | |||

SVR | −303.61 | −49.18 | |||

NA | −378.25 | −32.63 | |||

B | Comb. 3 | FFNN | −449.57 | −24.86 | |

ANFIS | −414.60 | −34.17 | |||

SVR | −376.59 | −24.86 | |||

NA | −392.45 | 20.18 | |||

C | Comb. 3 | FFNN | −380.00 | −44.62 | |

ANFIS | −312.49 | 24.70 | |||

SVR | −324.81 | −20.05 | |||

NA | −449.29 | −77.37 |

As is obvious from Table 7 and Figures 7–9, performance of the NA ensemble technique for stations A and C are slightly better than for station B at daily scale but at monthly scale all stations showed almost similar accuracy. According to the presented results in Table 7, the NA ensemble outperformed the other ensemble techniques in most of the stations. The overall performance of different ensemble techniques is almost the same but NA performed slightly better than the SA and WA ensemble techniques. Moreover, as shown by Figures 7–9, almost all models presented reliable outcomes but one model in one particular station captured the maximum amounts while in another it seized the minimum values well. In short, depending on the specific hydrologic features and spatial zone, the performance of a model may significantly or slightly fluctuate. Spread distribution of the dataset has a higher CV value and this means it is harder for the model to learn and predict. So when CV increases, the performance of AI models would decrease. Comparing the different stations’ performances and their CV values, it can be seen that in daily modeling all stations exhibited almost similar performance as they have the same CV values according to Table 2. Yet in monthly modeling, station B demonstrated slightly poor performance since its CV is greater than the CV values of stations A and C.

For better comparison of the models performance and to check the results for overfitting issue, the BIC values of single models for best combinations as well as NA model are tabulated in Table 8. According to Table 8, in both phases of calibration and verification, the obtained BIC values are close to each other. However, in most cases the NA model has a lower BIC value and leads to better results when compared to the single models.

### Results of the modeling via second scenario (multi-station modeling)

The main goal of modeling through the second scenario was to estimate the SSL of station C through a multi-station approach, meaning that the discharge and SSL of the other two upstream stations were used in the SSL modeling of station C (). The spatial locations of these three stations indicate a high correlation among them. To this end, different input combinations based on different physical interpretations were developed and thoroughly investigated. For instance, one combination only considered the data of station A as the determining station, while the other took station B, whilst another considered both of these stations. Table 9 presents the considered input combinations of the multi-station scenario for daily and monthly scales.

Modeling scale . | Combination No. . | Inputs . |
---|---|---|

Daily | 1 | , |

2 | , , | |

3 | , | |

4 | , , | |

5 | , , , | |

6 | , , , , | |

7 | , , | |

8 | , , | |

Monthly | 9 | , , |

10 | , , | |

11 | , | |

12 | , , | |

13 | , , , | |

14 | , , , , | |

15 | , , | |

16 | , , |

Modeling scale . | Combination No. . | Inputs . |
---|---|---|

Daily | 1 | , |

2 | , , | |

3 | , | |

4 | , , | |

5 | , , , | |

6 | , , , , | |

7 | , , | |

8 | , , | |

Monthly | 9 | , , |

10 | , , | |

11 | , | |

12 | , , | |

13 | , , , | |

14 | , , , , | |

15 | , , | |

16 | , , |

The optimum lag times of the combinations were obtained via trial and error process and as Table 9 indicates, the discharge and SSL of station A were considered with greater lags in comparison with station B which could be related to the remoteness of station A. In other words, it takes more time for the discharge and sediment load of station A to reach station C. Each combination was modeled by all methods, and then the dominant input combination was obtained through a modeling performance ranking. Results of FFNN, ANFIS, and SVR modeling for the best input combinations in both daily and monthly temporal scales are demonstrated in Table 10.

Modeling scale . | Inputs . | Model . | Structure^{a}
. | DC . | RMSE^{b}. | ||
---|---|---|---|---|---|---|---|

Calibration . | Verification . | Calibration . | Verification . | ||||

Daily | Comb. 8 | FFNN | (3-4-1, 50) | 0.82 | 0.70 | 0.06 | 0.08 |

ANFIS | Gaussian-2 | 0.81 | 0.70 | 0.06 | 0.08 | ||

SVR | 0.05, 0.1, 5 | 0.77 | 0.66 | 0.06 | 0.08 | ||

Monthly | Comb. 16 | FFNN | (3-2-1, 20) | 0.88 | 0.77 | 0.07 | 0.11 |

ANFIS | Triangular-2 | 0.86 | 0.74 | 0.08 | 0.12 | ||

SVR | 0.9, 0.1, 10 | 0.85 | 0.74 | 0.09 | 0.12 |

Modeling scale . | Inputs . | Model . | Structure^{a}
. | DC . | RMSE^{b}. | ||
---|---|---|---|---|---|---|---|

Calibration . | Verification . | Calibration . | Verification . | ||||

Daily | Comb. 8 | FFNN | (3-4-1, 50) | 0.82 | 0.70 | 0.06 | 0.08 |

ANFIS | Gaussian-2 | 0.81 | 0.70 | 0.06 | 0.08 | ||

SVR | 0.05, 0.1, 5 | 0.77 | 0.66 | 0.06 | 0.08 | ||

Monthly | Comb. 16 | FFNN | (3-2-1, 20) | 0.88 | 0.77 | 0.07 | 0.11 |

ANFIS | Triangular-2 | 0.86 | 0.74 | 0.08 | 0.12 | ||

SVR | 0.9, 0.1, 10 | 0.85 | 0.74 | 0.09 | 0.12 |

For FFNN, the format “a-b-c, d” presents the number of input layers-hidden layers-output layer neurons, and epochs. For ANFIS, the different membership functions are represented as “MF-a”, indicating the MF and number of membership functions. For SVR, “a, b, c” denote γ, ε and c, respectively.

^{b}RMSE has no dimension since all data are normalized.

The main reason for examining various combinations in the modeling was to provide a tool for the physical interpretation of the outcomes to evaluate the interaction degree between stations and to see which river transports the major SSL to station C. Throughout the modeling with different input combinations, it was found out that employing the sole data of stations A and B (input combinations 2 vs 4 and 10 vs 12) demonstrated almost similar performance. This indicates the same interconnection of stations A and B with station C, unlike the early speculations. The initial surmise was that station A may have less correlation with station C (compared to station B) because of the long distance between them. This anomalous outcome could be related to the poor condition of land cover which consequently adds to the sediment load through the soil erosion process. Furthermore, considering the verification phase DC values of modeling by different input combinations, Comb. 6 and Comb. 8 in daily modeling and Comb. 14 and Comb. 16 at monthly scale showed better performance among other input combinations. Although Comb. 6 and Comb. 8 outperformed the other two input combinations, the 8th and 14th input combinations were chosen as the best data fusion for daily and monthly modeling, respectively, since they have fewer parameters and no SSL data were involved in their formulation. It is worth noting that, unlike single station modeling in which daily modeling exhibited better results, in the multi-station scenario, the monthly outcomes showed more reliable results than the daily modeling as the stations in monthly modeling interact with each other in a long lag time.

Prediction results for Comb. 8 and Comb. 16 were fed into ensemble models, as in the ensemble methodology conducted in scenario 1. Table 11 shows the obtained results of the ensemble approach for the chosen input combinations and Table 12 shows the BIC values of single models and NA ensemble for better comparison of the models’ performances. Additionally, Figures 10 and 11 respectively demonstrate daily and monthly observed versus computed values of SSL for multi-station modeling of station C in the verification step, and scatter plots of all proposed models for station C are presented in Figure 12.

Modeling scale . | Inputs . | Ensemble Method . | Model Structure^{a}
. | DC . | RMSE^{b}. | ||
---|---|---|---|---|---|---|---|

Calibration . | Verification . | Calibration . | Verification . | ||||

Daily | Comb. 8 | SA | _ | 0.81 | 0.71 | 0.06 | 0.07 |

WA | 0.340, 0.340, 0.320 | 0.81 | 0.71 | 0.06 | 0.07 | ||

NA | (3-2-1, 80) | 0.83 | 0.73 | 0.06 | 0.07 | ||

Monthly | Comb. 16 | SA | _ | 0.86 | 0.76 | 0.09 | 0.12 |

WA | 0.342, 0.329, 0.329 | 0.86 | 0.76 | 0.09 | 0.12 | ||

NA | (3-3-1, 40) | 0.93 | 0.79 | 0.06 | 0.10 |

Modeling scale . | Inputs . | Ensemble Method . | Model Structure^{a}
. | DC . | RMSE^{b}. | ||
---|---|---|---|---|---|---|---|

Calibration . | Verification . | Calibration . | Verification . | ||||

Daily | Comb. 8 | SA | _ | 0.81 | 0.71 | 0.06 | 0.07 |

WA | 0.340, 0.340, 0.320 | 0.81 | 0.71 | 0.06 | 0.07 | ||

NA | (3-2-1, 80) | 0.83 | 0.73 | 0.06 | 0.07 | ||

Monthly | Comb. 16 | SA | _ | 0.86 | 0.76 | 0.09 | 0.12 |

WA | 0.342, 0.329, 0.329 | 0.86 | 0.76 | 0.09 | 0.12 | ||

NA | (3-3-1, 40) | 0.93 | 0.79 | 0.06 | 0.10 |

For WA, the format “a, b, c,” denotes the weights of the FFNN, ANFIS, SVR, models. For NA, the format “a-b-c, d” in represents the number of input layers-hidden layers-output layer neurons, and epochs.

Modeling scale . | Inputs . | Model . | BIC . | |
---|---|---|---|---|

Calibration . | Verification . | |||

Daily | Comb. 3 | FFNN | −15,205.80 | −4,454.27 |

ANFIS | −15,231.62 | −4,476.80 | ||

SVR | −15,214.41 | −4,461.78 | ||

NA | −15,291.87 | −4,773.20 | ||

Monthly | Comb. 3 | FFNN | −421.54 | −87.40 |

ANFIS | −330.00 | −28.95 | ||

SVR | −329.57 | −45.33 | ||

NA | −423.33 | −72.64 |

Modeling scale . | Inputs . | Model . | BIC . | |
---|---|---|---|---|

Calibration . | Verification . | |||

Daily | Comb. 3 | FFNN | −15,205.80 | −4,454.27 |

ANFIS | −15,231.62 | −4,476.80 | ||

SVR | −15,214.41 | −4,461.78 | ||

NA | −15,291.87 | −4,773.20 | ||

Monthly | Comb. 3 | FFNN | −421.54 | −87.40 |

ANFIS | −330.00 | −28.95 | ||

SVR | −329.57 | −45.33 | ||

NA | −423.33 | −72.64 |

According to Table 11, ensemble techniques enhanced the prediction results in both calibration and verification phases and the results indicated the superiority of the NA ensemble compared to the other two ensemble techniques. NA ensemble has been successfully employed as a multi-model combination technique in hydrological modeling (e.g. Elkiran *et al.* 2018; Sharghi *et al.* 2018; Nourani *et al.* 2021). According to Figures 10–12 almost all models presented acceptable estimations considering the fact that in this scenario none of the models employed the observed SSL values of station C as input data and also the selected combinations for the ensemble were only based on streamflow data. However, this multi-station scenario could be employed when the SSL data of station C is not available due to financial or technical issues. This multi-station modeling outlined the high correlation of the stations and indicated the important role of the physical aspects and geomorphology features of the study area (e.g. land use and soil type) to be taken into consideration in SSL and erosion modeling. For instance, as Table 1 states, the upstream basin of station A has less forest and agricultural lands compared to station B which causes the basin of station A to be subjected to severe soil erosion and to carry more sediment to station C. This was confirmed through the multi-station modeling scenario of this study as stations A and B demonstrated almost the same correlation with SSL value of station C, despite the great distance of station A.

To have a visualized comparison of the models’ performances in the verification step, Radar diagrams of station C at daily scale for both modeling scenarios are presented in Figure 13.

## CONCLUSIONS

In this research, due to the practical restrictions of physical and conceptual models, soft computing methods of FFNN, ANFIS, SVR, as well as ensemble techniques were employed for predicting SSL. To this end, discharge and SSL data of three gauging stations in two different hydrological regions of the United States were used via two distinct modeling scenarios. The first scenario, using single station modeling, was developed to predict the SSL in the current time step () for all stations, using each station's own data. In the next modeling scenario, using multi-station modeling, observed data of two upstream stations (i.e. A and B) were utilized to estimate the SSL of a downstream station C () without using the sediment record of station C but employing its streamflow data.

Several input combinations were examined for both scenarios and each input combination's performance was evaluated through modeling. Then the outcomes of the best input combination were fed into ensemble units using the SA, WA, and NA ensemble techniques to increase the modeling efficiency and to address the issues with choosing one single model which is the best fit for all stations in any circumstance. Since every single model has its own merits and demerits, opting for the best method is not an easy task and this was proved by the acquired outcomes of this study. The results indicated that each AI method may exhibit different accuracies depending on the spatial and temporal variations of the stations. In individual daily modeling of station B for example, ANFIS showed poor performance in the verification step compared to the FFNN and SVR models but in stations A and C it performed as well as other methods. All in all, the obtained results of this study indicated the outperformance of FFNN among other AI models.

For both modeling scenarios, employing ensemble techniques improved the prediction efficiency and the NA ensemble showed more reliable performance than SA and WA techniques since linear averaging always gives a result that is higher than the minimum value and less than the maximum value in the set. Ensemble models affected the monthly scale estimation more than the daily modeling. Ensemble techniques could increase the prediction accuracy of single station modeling up to 7% in training and up to 20% in the verification phase of the daily modeling and up to 10 and 20% in the training and verification phases of monthly modeling. Yet in multi-station modeling, it improved the prediction results up to 5% in the training and verification step of the daily time scale and up to 8% in both the training and verification steps of the monthly modeling. Furthermore, in this study, daily modeling exhibited more robust outcomes in the single station scenario while in the multi-station scenario monthly scale outperformed the daily modeling which highlights the significant function of SSL modeling in various time scales. The findings of this study could pave the way for reliable prediction of SSL by employing AI-based ensemble models and facilitate the evaluation of rivers’ fluvial and morphological processes. Moreover, the findings provide technical support for river hydraulics and environment as well as for river management.

Every model has its own limitations and disadvantages. In this case, the main disadvantage of the applied AI models of this study is their being black-box models in which the computing system is opaque and consequently does not allow the user to sequentially eliminate possible explanatory variables that do not contribute to model fitting (Kisi *et al.* 2012). Moreover, the AI models are case sensitive meaning that they demonstrate different performances for different case studies. Also they are highly dependent on the quality and quantity of the used dataset. Another limitation of this study is using only black-box models in developing the ensemble unit. So for future studies, it is also suggested to use physical-based models along with AI models to develop a comparative study and to employ the advantages of the models, as well as modeling other tributaries of the study area in order to strengthen the credibility of the proposed methodology. Also, it is suggested to employ other AI methods like genetic programming (GP) and to ensemble AI models with conceptual methods.

## CONFLICTS OF INTEREST

The authors report no conflicts of interest

## DATA AVAILABILITY STATEMENT

Data cannot be made publicly available; readers should contact the corresponding author for details.