Water Flow and Sediment Flux Forecast in the Chókwè Irrigation Scheme, Mozambique

Lateiro Salvador De Sousa^{1}^{*},
Raphael Muli Wambua^{2},
James Messo Raude^{3},
Benedict Mwavu Mutua^{4}

Show more

1. Introduction

Water flow and suspended sediment forecasting is an important problem in water resources management [1] [2] [3]. In modelling of the river streams and unlined irrigation canals the evacuation of suspended sediment is quite significant. To increase the active lifetime of irrigation canals, the dead volume capacity should be minimized. A considerable amount of sediment carried by river flow can be transported through the irrigation canal from upstream to downstream. Controlling the sediment flushing, which is, limiting the sediment concentration of the sluiced waters, is the main mitigation measure that can be adopted to decrease the impact of flushing on the downstream environment. However, on the one hand, controlling the sediment concentration while flushing a reservoir is operationally difficult and, on the other hand, thresholds for sediment concentration balancing technical, economic, and environmental aspects is controversial. Forecasting of the suspended sediment can provide information for management of dam sluiceway, as well as, for the river streams and irrigation canals. Suspended sediment is a key variable in a river because it is related to contaminant transport, water quality, reservoir sedimentation, silting, soil erosion and loss, and has clear ecological and recreational impacts. The forecasting studies in the literature are restricted to one-step ahead forecasting [4]. Within the fluvial environment, sediment and its associated contaminants can be in transport in suspended form or in temporary channel bed storage. Contaminants stored in channel beds are, however, considered to be important for benthic habitat quality. Concentrations of pollutants in suspended and/or temporarily stored sediment have been shown to be controlled by proximity to sediment sources [5] [6]. Additionally, there exist many papers in the literature that compare the performance of models for different engineering problems [4] [7] [8] [9].

Water flow and sediments on the canal stream can behave differently according to the changes observed in space and time. Within space alterations, it can be seen in two levels, across distance from one point to the other, as well as with depth at surface, middle and bed points [10]. On the other hand, within time, trends may change with time (days, months, season, years or even decades). Different studies have been carried out to find out the spatial and temporal variability of water flow and sediments. In a study by [11] it was found that the intensity of transport varies with channel location, and that coarse bed load particles appeared to be differentially transported in some canal’s areas than on others [12]. Moreover, research work found that the presence of vegetation in canals increased the grain resistance that transports sediment [13] [14] [15]. The authors concluded that the grain resistance is a function of vegetation density, and therefore the dense the vegetation, the less the erosion [16]. This study examines the employment of two methods, Mann-Kendall test and ARIMA model for forecasting of water flow and sediment flux in the Chókwè Irrigation Scheme, in Mozambique. The autoregressive integrated moving average (ARIMA) model is considered for one-step ahead forecasting of sediment series.

1.1. Hydraulic Functions

Sediment concentration in canals varies with time and location. Some canals have very little or no sediments while others suffer from high concentrations throughout the year or in certain months [17]. Two aspects are considered when dealing with water flow in the irrigation canals. Firstly, the operational aspects where the water flow becomes non-uniform and unsteady due to changing water requirements and gates operation to fulfill the water demand and to keep water level as it is required for the field’s needs [18] [19] [20]. Secondly, the sediment transport aspects, where the changes in water flow in time and space are faster than changes in morphology of sediment. Therefore, the interrelation between water flow and sediment transport can be illustrated as one-dimensional event without changing the cross-sectional shape [21]. Additionally, the trends-dynamics of water and sediment fluxes can also be described as presented in the following equations:

The following continuity equation for water movement was put forth by [21]:

$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=0$ (1)

where: *A* = area (m^{2}); *Q* = discharge (m^{3}); *t* = time (s); *x* = distance (m).

Additionally, a dynamic equation for water movement was presented by the same author:

$\frac{\partial h}{\partial x}+\frac{{V}^{2}}{{C}^{2}R}+\frac{\partial z}{\partial x}+\frac{V}{g}\times \frac{\partial V}{\partial x}+\frac{1}{g}\times \frac{\partial V}{\partial t}=0$ (2)

where: *h* = water depth (m); *V* = mean velocity (m^{1/2}/s); *C* = Chezy’s coefficient (m/s^{2}); *R* = hydraulic radius (m); *z* = bottom level above datum (m); *g* = gravity acceleration (m^{2}/s); *t* = time (s); and x = distance (m).

Equations (3) and (4) describe the conservation of mass and momentum. They are also known as Saint-Venant equations for continuity and dynamic unsteady flows and water-flow related factors. For sediment related factors, the equations are as follows:

1.2. Sediment Models

In this case, the friction factor predictor is given as a function of

$C=f\left({d}_{50},V,h,{S}_{o}\right)$ (3)

and the continuity equation for sediment transport is

$\left(1-p\right)\times B\times \frac{\partial z}{\partial t}+\frac{\partial {Q}_{s}}{\partial x}=0$ (4)

where: *t* = time (s); *Q _{s}* = sediment discharge (m

Finally, the sediment transport equation is presented as a function of

${Q}_{s}=f\left({d}_{50},V,h,{S}_{o}\right)$ (5)

These equations are related to each other since sediment transport and water flow are interrelated. For example, the roughness coefficient is influenced by water flow, while sediment transport is affected by the water flow [22]. The unsteady flow condition, on the other hand, in irrigation canal is assumed to be quasi-steady; hence ∂*A*/∂*t* and ∂*V*/∂*t* can be neglected. The continuity and dynamic equations become

$\frac{\partial Q}{\partial x}=0$ (6)

$\frac{\text{d}h}{\text{d}x}=\frac{{S}_{o}-{S}_{f}}{1-F{r}^{2}}$ with $Fr=\frac{V}{\sqrt{g\times h}}$ (7)

For the uniform flow, there is no change in water depth. Hence,

$\frac{\text{d}h}{\text{d}x}=0\to {S}_{o}={S}_{f}$ (8)

making the uniform flow new equation

$v=\frac{1}{n}\times {R}^{1/3}\times {S}_{f}^{2/3}$ or $v=C\sqrt{R\times {S}_{f}}$ (9)

where: *A* = area (m^{2}); *Q* = discharge (m^{3}); *h* = water depth (m); *V* = mean velocity (m^{1/2}/s); *C* = Chezy’s coefficient (m/s^{2}); *R* = hydraulic radius (m); *z* = bottom level above datum (m); *g* = gravity acceleration (m^{2}/s); *t* = time (s); *Q _{s}* = sediment discharge (m

1.3. Mann-Kendall Models

Trend detection in hydrologic and water quality time series has received considerable attention in the recent past. In a number of studies on water quality data in lakes and streams [23] [24] [25] and streamflow data [26] [27] [28] a number of parametric and non-parametric tests have been applied for trend detection. Both parametric and non-parametric tests are commonly used [29] [30]. Parametric trend tests are more powerful than non-parametric ones, but they require data to be independent and normally distributed [31]. On the other hand, non-parametric trend tests require only that the data be independent and can tolerate outliers in the data [32].

One of the widely used non-parametric tests for detecting trends in the time series is the Mann Kendall test [33] [34]. The Mann-Kendall trend test is derived from a rank correlation test for two groups of observations proposed by [34]. In the Mann-Kendall trend test, the correlation between the rank order of the observed values and their order in time is considered. The null hypothesis for the Mann-Kendall test is that the data are independent and randomly ordered, that is, there is no trend or serial correlation structure among the observations. However, in many real situations the observed data are autocorrelated.

The autocorrelation in observed data will result in misinterpretation of trend test results. [35] states that: “*Positive* *serial* *correlation* *among* *the* *observations* *would* *increase* *the* *chance* *of* *significant* *answer*, *even* *in* *the* *absence* *of* *a* *trend*.” A closely related problem that has been studied is the case where seasonality exists in the data [36]. By dividing the observations into separate classes according to seasons and then performing the Mann-Kendall trend test on the sum of the statistics from each season, the effect of seasonality can be eliminated. This modification is called the seasonal Kendall test [36] [37] [23]. Although the seasonal test eliminates the effect of dependence between seasons, it does not account for the correlation in the series within seasons [37]. The same problem exists when yearly data are analysed since they are often significantly autocorrelated.

The rank correlation test [34] for two sets of observations
$X={x}_{\text{l}},{x}_{\text{2}},\cdots ,{x}_{n}$ and
$Y={y}_{l},{y}_{2},\cdots ,{y}_{n}$ is formulated as follows. The statistic *S* is calculated as in Equation (10):

$S={\displaystyle {\sum}_{i<j}{a}_{ij}\times {b}_{ij}}$ (10)

where:

${a}_{ij}=\mathrm{sgn}\left({x}_{j}-{x}_{i}\right)=\{\begin{array}{l}1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}{x}_{i}<{x}_{j}\\ 0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{x}_{i}={x}_{j}\\ -1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}{x}_{i}>{x}_{j}\end{array}$ (11)

and *b _{ij}* is similarly defined for the observations in

$E\left(S\right)=0$ (12)

$var\left(S\right)=n\left(n-1\right)\left(2n+5\right)/18$ (13)

If the values in *Y* are replaced with the time order of the time series *X*, *i.e.* 1, 2, ∙∙∙,* n*, the test can be used as a trend test [33]. In this case, the statistic S reduces to that given in Equation (14):

$S={\displaystyle {\sum}_{i<j}{a}_{ij}}={\displaystyle {\sum}_{i<j}\mathrm{sgn}\left({x}_{j}-{x}_{i}\right)}$ (14)

with the same mean and variance as in Equations (12) and (13). [34] gives a proof of the asymptotic normality of the statistic *S*. The significance of trends is tested by comparing the standardized test statistic *Z* = *S*/[*var*(*S*)]^{0.5} with the standard normal variate at the desired significance level. The derivation of the mean and variance of *S* is discussed in detail by [34]. If *X* is normally distributed with mean *g* and variance *σ*^{2}, then (*x _{j}* −

1.4. ARIMA Models

In theory, Auto-Regressive Integrated Moving Average (ARIMA) models, are the most general class of models for forecasting a time-series. An ARIMA model can be stationarized by transformations such as differencing and logging. In fact, the easiest way to think of ARIMA models is as fine-tuned versions of random-walk and random-trend models: the fine-tuning consists of adding lags of the differenced series and/or lags of the forecast errors to the prediction equation, as needed to remove any last traces of autocorrelation from the forecast errors. Lags of the differenced series appearing in the forecasting equation are called “autoregressive” terms; lags of the forecast errors are called “moving average” terms; and a time-series which needs to be differenced to be made stationary is said to be an “integrated” version of a stationary series [38]. Random-walk and random-trend models, autoregressive models, and exponential smoothing models (that is, exponential weighted moving averages) are all special cases of ARIMA models. A non-seasonal ARIMA model is classified as an “ARIMA (*p,* *d,* *q*)” model, where: *p* is the number of autoregressive terms, *d* is the number of non-seasonal differences and *q* is the number of lagged forecast errors in the prediction equation [39]. To identify the appropriate ARIMA model for a time-series, it begins by identifying the order(s) of differencing needed to stationarize the series and to remove the gross features of seasonality, perhaps in conjunction with a variance-stabilizing transformation such as logging or deflating. If it stops at this point and predicts that the differenced series is constant, it will have merely fitted a random-walk or random-trend model.

However, the best random-walk or random-trend model may still have autocorrelated errors, suggesting that additional factors of some kind are needed in the prediction equation. ARIMA model forecasting includes three basic steps: model identification, parameter estimation and forecasting. ARIMA model parameter selection is based on the autocorrelation -function linear relation between observation pairs; and partial-autocorrelation-function conditional correlation with intervening observations removed. According to [40], the general procedure for ARIMA model selection and calibration includes: 1) stationarity conditions checking; 2) autocorrelation function checking to choose the *p* value; 3) partial autocorrelation function checking to choose the Q value 4) identifying the ARIMA (*p*, *q*) model; 5) estimation; 6) residual diagnostics. In this study, several trials were made to choose the optimal ARIMA model parameters. The model parameters that meet the statistical residual diagnostic checking were chosen in the ARIMA forecasting model [38].

Generally, for the Box-Jenkins methodology [39], the auto-regressive moving average (ARMA (*p*, *q*)), or auto-regressive integrated moving average (ARIMA (*p*, *d*, *q*)) models are often applied for time series forecasting. However, the application of ARMA model requires the time series to be stationary; that is to say, the algorithm of ARMA assumes that the process remains in equilibrium about a constant mean level. If the series are nonstationary or have obvious trend variability, the ARIMA model based on difference process can be used [41]. In this work, the augmented Dickey-Fuller (ADF) test [42] is used to test the stationarity in the original annual runoff time series and the decomposed annual runoff time series [43].

As it is known, the ARMA model consists of three main steps: model identification, parameter estimation and application. Among these three steps, the identification step is important, and includes two stages: 1) if it is necessary, appropriate differencing of the series is performed to achieve stationary and normality, 2) the order of the AR and MA parts of ARMA model is identified. [39] employed the autocorrelation function (ACF) and the partial autocorrelation function (PACF) of the sample data as the basic tools to identify the order of the ARIMA model. If sample data is an AR (*p*) model, the PACF cut-off is at lag *p*, On the other hand, if sample data is a MA (*q*) process, the ACF has a cut-off at lag *q*. However, the PACF and ACF method is not useful when dealing with mixed ARMA processes. Simple inspection of the graphs of the ACF and the PACF would not, in general, give clear values of *p* and *q* for mixed models [44]. Some other identification methods have been presented based on the information-theoretic approaches, such as the Akaike’s Information Criterion (AIC) [45], the Bayesian Information Criterion (BIC) [46], the final prediction error criterion (FPE) [47] and others. In this work, the best-fitted model is selected according to AIC. Once an appropriate model is chosen, the parameters of the model must be estimated. This can be accomplished using a nonlinear optimization procedure. The application and analysis are presented as follows:

The MA (1) model is in the form off:

${x}_{i}={e}_{t}+\theta {e}_{t-1}$ (15)

where *e _{t}* is
$N\left(0,{\sigma}_{e}^{2}\right)$. In this case the autocorrelation function

$\rho \left(i\right)=\{\begin{array}{l}1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}i=0\\ \frac{\theta}{1+{\theta}^{2}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}i=1\\ 0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}i>1\end{array}$ (16)

The AR (1) model is of the form of:

${x}_{t}=\phi {x}_{t-1}+{e}_{t}$ (17)

In this case, the autocorrelation function is given by:

$\rho \left(i\right)={\phi}^{\left|i\right|}$ (18)

2. Materials and Methods

2.1. Study Area

This study was undertaken at the Chókwè Irrigation Scheme (CIS), which is in the Limpopo River Basin (LRB), Chókwè District, Gaza Province in Mozambique. The scheme is located at the Lower Limpopo River Sub-Basin (LLRSB), covers an area of approximately 84,981 km^{2}, and lies from latitudes 24˚04'32'' South to 25˚01'35'' South and longitudes 32˚40'11'' East to 33˚37'14'' East. The CIS is largely dry, with rainfall averaging between 500 and 600 mm/year (with a mean of 530 mm/year). Rainfall events are concentrated between October and March. The population density is 88.78 inhabitants/km^{2}, according to the 2017 Census results [48]. The Limpopo River originates from Central southern Africa and flows generally eastwards to the Indian Ocean, traversing a terrain encompassing an altitude of 1600 m in South Africa (Drakensberg Mountains) relative to the sea level in Mozambique [49]. Its length and drainage area are estimated to be 1750 km long and 430,000 km^{2}, respectively, while the mean annual discharge at its mouth in Mozambique is 170 m^{3}/s [50]. Figure 1 presents the map for the study site.

2.2. Data Acquisition

The CIS is the main irrigation scheme in Mozambique and sources water from Limpopo River at approximately 45 m^{3}/s. Water is diverted to unlined canals benefiting more than 12,000 farmers tilling approximately 33,000 hectares for food production [51] [52]. CIS is used to deviate, store, manage and distribute water to the local producers, which is possible using two hydraulic structures: Massingir dam and Macarretane weir, both located at upstream. Agriculture is the main activity in the region and constitutes the backbone of the district, producing rice, maize, and vegetables. Nearly 90% of the irrigation scheme is irrigated by gravity. Gravity flow system is the main form of water application through furrow and flood methods. The main crops grown in the region are rice (for wet season), vegetables (dry season) and maize (in both seasons).

Nine sampling points were established across the main canal of CIS for field data collection on the bedload sediment. Three stations where established at each section of the canal, namely at upstream (Montante section), midstream (Sul section) and downstream (Rio section). At the Montante section, sampling points were Macarretane Weir Intake, Railways-Node and FIPAG bridges. At

Figure 1. Map of Chókwè district in Gaza province, Mozambique.

the Sul section were Lionde, Massawasse and Conhane bridges. At the Rio section Nico, Muianga and Marrambajane bridges were the sampling stations. The study covered an estimated distance of around 60 km. The study did not cover some parts of the CIS radius, including Canal Esquerdo, Canal of Nwachicoluane and the most downstream of Rio section, below Marrambajane bridge towards the Chiguidela region.

2.3. Sediment Concentration

The determination of temporal water flow distribution trends at CIS was performed using historical data on the water flow from 2004 to 2019 available at HICEP. From the water flow data, correlation was established to generate a relation function between both data (water flow and sediment). These outcomes were used to generate the temporal sediment amounts at CIS. To estimate the suspended-sediment concentration by the interpolation method, USGS proposed the following equation:

${Q}_{s}={Q}_{w}\times {C}_{s}\times k$ (19)

where: *Q _{s}* = suspended-sediment discharge (tonnes per day);

For determination of water flow and sediment the study considered different distances from the main intake toward the lowest point in the canals. Mann-Kendall test was also considered for determination of temporal trends. Two hypotheses where tested for MK: null hypothesis (H_{0}), there was no trend in the series, and the alternative hypothesis (H_{1}), there was a trend in the series. When the computed p-value showed to be greater than the significance level alpha = 0.05, then the null hypothesis (H_{0}) was not rejected. And when the computed p-value was lower than the significance level alpha = 0.05, the null hypothesis H_{0}, would be rejected, and the alternative hypothesis (H_{1}) accepted.

Prior to implementing any time series analysis, in ARIMA test, the data was evaluated for any dominating trend signals. A bars/column, auto-correlogram function (ACF) and partial auto-correlogram function (PACF) technique was employed to determine the linear trend, if any, and the seasonality of the data. If present, the linear trend was removed by a simple linear regression technique. Since seasonality can easily be identified in the domain, it was not necessary to remove the seasonal signal prior to further analysis. Therefore, the time domain and frequency domain analyses were performed using the linear detrended data to quantify the persistence of total water discharge in the scheme.

3. Results and Discussion

3.1. Mann-Kendall Trends for Water Discharge

Figures 2(a)-(r) presents Mann-Kendall (MK) test at 95% of confidence interval, whereby the Sen’s slope and intercept are also given for all the nine sampling

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o) (p) (q) (r)

Figure 2. Mann-Kendall test and Sen’s slope analysis for water discharge at CIS all stations.

stations in the dry (DS) and wet (WS) seasons.

The sampling stations that were tested performed well by presenting an indication of existence of trends. However, from the results only Lionde WS, Massawasse DS, Muianga WS and Marrambajane WS, showed trends that were statistically significant. Among the stations with significant trends, only Lionde and Massawasse presented increasing trend, whilst others did not. These findings are consistent with of [53], where decreases in streamflow were observed in Urmia Lake, in Iran. Another study by [30], found a decreasing trend of the annual streamflow yet at the 5% significance level at the upper Senegal River Basin. Each sampling station showed difference in terms of the highest water flow amount passing through it. For example, it can be noted that the only stations with water flow above 2 m^{3}/s were the Intake and the Lionde’s. This is because these two points are, respectively, the general inlet to the whole scheme, and a diversion point to other channels, therefore, dealing with significant amount of water. The remaining stations did not present considerable water flow. Sampling stations as Node, Massawasse, Conhane and Marrambajane had a maximum water flow of 2 m^{3}/s, whilst Nico and Muianga, just below 1 m^{3}/s. Finally, FIPAG station was well below 0.1 m^{3}/s.

Each MK plot presents a Sen’s slope, which captures the magnitude of the changes within the trends along the studied year. So, Sen’s slope offers the insight of the magnitude of the trend with time. In this case, for the significant trends, they are either positive or negative. From Table 1, at Lionde, Sen’s slope was approximately 0.001 m^{3}/s (1 l/s). While in Massawasse, Sen’s slope was around 0.00005 m^{3}/s, which inclines towards zero (or 50 ml/s). The Muianga and Marrambajane, on the other hand, presented negligible value for their Sen’s slope. This means that in these stations, changes in magnitude were small and negligible in face of the total amount of water in the system. From these results, one can see that the increase in water discharge is of higher magnitude than the

Table 1. Mann-kendall test for water discharge trends for both seasons from 2004-2005 to 2018-2019 period.

*The hydrologic year in Mozambique starts in September/October at the second semester of the year and ends up in August/September of following year. Therefore, where 2004 is stated, this means that the hydrologic year is 2004-2005, and 2018-2019. **Test interpretation: Null hypothesis (H_{0}): There is no trend in the series; Alternative hypothesis (H_{a}): There is a trend in the series. When the computed p-value is *greater* than the significance level alpha = 0.05, one cannot reject the null hypothesis H_{0}. And when the computed p-value is *lower* than the significance level alpha = 0.05, one should reject the null hypothesis H_{0}, and accept the alternative hypothesis H_{a}.

decrease in water discharge. Also, it shows the upstream of the scheme having increasing discharge while the downstream presents decreasing discharges.

3.2. Mann-Kendall Trends for Sediment Flux

The maximum sediment flux generated at each station was, respectively, 1000 ton/day and 70 ton/day, for the Intake and Offtake stations. MK test for sediment discharge trends were not significant at 95% significance level, except for the Offtake in WS. However, though with no significant trend, the Sen’s slope tends to increase at the Intake and decrease at the Offtake. The magnitude of these variations for Intake were 0.027 and 0.041 ton/day (27 and 41 kg/day), during the WS and DS, respectively. However, for Offtake, the magnitudes were −0.009 and −0.007 ton/day (−9 and −7 kg/day), for WS and DS, respectively. Marrambajane station was the only station that presented significant trends for both water discharge and sediment. These results allow us to infer that in other stations such as Lionde in the WS, Massawasse in the DS, Muianga in the WS, significant trends could be observed as well. Figure 3 presents Mann-Kendall Test and Sen’s Slope analysis of sediment discharge for Intake and Offtake. Additional information about it is provided in Table 2 which presents the MK analysis for sediment discharge.

(a) (b) (c) (d)

Figure 3. Mann-Kendall test and Sen’s slope analysis of sediment discharge for intake and offtake.

Table 2. Mann-kendall test for sediment discharge trend at intake and Offtake for both seasons for the 2004-2005 to 2018-2019 period.

*The hydrologic year in Mozambique starts in September/October at the second semester of the year and ends up in August/September of following year. Therefore, where 2004 is stated, this means that the hydrologic year is 2004-2005, and 2018-2019. **Test interpretation: Null hypothesis (H_{0}): There is no trend in the series; Alternative hypothesis (H_{a}): There is a trend in the series. When the computed p-value is *greater* than the significance level alpha = 0.05, one cannot reject the null hypothesis H_{0}. And when the computed p-value is *lower* than the significance level alpha = 0.05, one should reject the null hypothesis H_{0}, and accept the alternative hypothesis H_{a}.

3.3. ARIMA for Water Discharge

A complementary analysis of Auto Regressive Integrated Moving Average (ARIMA (*p*, *d*, *q*)) was performed, using the model parameters of *p* = 1/*d* = 0/*q* = 0/*P* = 0/*D* = 0/*Q* = 0/*s* = 0 for model training, at a confidence interval of 95%. Therefore, the performed analysis was ARIMA (1, 0, 0) also called first-order autoregressive model for pre-modelling. Where *p* is the number of autoregressive terms, *d* is the number of non-seasonal differences needed for stationarity, and *q* is the number of lagged forecast errors in the prediction equation. The test considered preliminary estimation based on Yule-Walker, with the optimization likelihood (Convergence = 0.00001/Iterations = 500). The number of validation and prediction was 5. Once the model training was performed, then the actual modelling was considered, the results are as presented in Figure 4, where WS stands for wet season and DS for season.

In this work, the Augmented Dickey-Fuller (ADF) test was used to test the stationarity in the original annual water discharge time series and the decomposed annual water discharge series. It was observed that all the series were non-stationary. Therefore, these series were eligible for ARIMA modelling to be applied over them, except for series represented by Node and Muianga, in both seasons, and Lionde in the WS. This therefore indicate that, except for these five stations, all other dataset had no unit root effect in the sample dataset, and it needed no differentiation in the original annual water discharge time series. Hence, requiring no decomposition of annual water discharge time series from nine of the sampling stations. They were suitable for ARMA model. After the first order difference of the original data of Node, Muianga and Lionde, the ADF test showed that the transformed data was stationary. Therefore, a process of model identification in the second stage was performed. Once the stationarity test was performed, then the following step was to identify the model for each time series, as shown on Table 3. The adopted structure and parameters for the ARIMA models for water discharge are presented in Table 4.

The ARMA model consisted of model identification, parameter estimation and application. Among these steps, the identification step is important, and includes two stages: 1) if it is necessary, appropriate differencing of the series is performed to achieve stationarity and normality, 2) the order of the AR and MA parts of ARMA model is identified. [54] employed the autocorrelation function (ACF) and the partial autocorrelation function (PACF) of the sample data as the basic tools to identify the order of the ARIMA model. When the sample data is an AR (*p*) model, the PACF cut-off is at lag *p*. On the other hand, if sample data is a MA (*q*) process, the ACF has a cut-off at lag *q*. However, the PACF and ACF method was not useful when dealing with mixed ARMA processes. Simple inspection of the graphs of the ACF and the PACF would not, in general, give clear values of *p* and *q* for mixed models [54]. These can be seen from Figure 5-13, for the nine sampling stations. Some other identification methods have been presented based on the information-theoretic approaches, such as the Akaike’s

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o) (p) (q) (r)

Figure 4. ARIMA test analysis for each sampling stations.

Figure 5. ACF, PACF, ACFR and PACFR for water flow at Intake in the DS and WS.

Table 3. Dickey-fuller test for stationarity assessment of water discharge dataset relevant for ARIMA model.

WS = Wet season; DS = Dry season.

Information Criterion (AIC) [55] [56] [57] the Bayesian Information Criterion (BIC) [54], the Final Prediction Error criterion (FPE) [47], among others. In this work, the best fitted model was selected according to AIC, whereby, the smaller the coefficient, the better the model. Once an appropriate model was chosen, the parameters of the model were estimated by reading the value of the best fitted

Table 4. Adopted structure and parameters of ARIMA model for water discharge at CIS based on AIC.

Test performed at confidence intervals of 95%; Optimization at Likelihood (Convergence = 0.00001/Iterations = 500).

*p*- and *q*-values, in a trial and error procedure. Therefore, models for Node (DS and WS), Massawasse (DS), Conhane (WS) and Marrambajane (WS) presented a good fitted model, with AIC of 16.582, 53.315, 12.206, 29.828 and 13.145, respectively. Additionally, Massawasse (WS) and Lionde (DS) with 101.235 and 408.291, respectively also offered a good fit. Other models were relatively acceptable since the negative values were found.

The results of the auto correlation analysis for the period 2004-2019 are presented in Figures 6-14, for the annual water flow, in the Intake, Node, FIPAG, Lionde, Massawasse, Conhane, Nico, Muianga and Marrambajane sampling stations, respectively, in the dry and wet seasons. From these figures, we can conclude that statistically, some series were auto correlated at 5% with significance level at certain lags. For example, the autocorrelation is significant at 5% at lag 1, 2 and 3 levels for the Intake (in both seasons), FIPAG (DS), Lionde and Nico

Figure 6. ACF, PACF, ACFR and PACFR for water flow at Node in the DS and WS.

Figure 7. ACF, PACF, ACFR and PACFR for water flow at FIPAG in the DS and WS.

Figure 8. ACF, PACF, ACFR and PACFR for water flow at Lionde in the DS and WS.

Figure 9. ACF, PACF, ACFR and PACFR for water flow at Massawasse in the DS and WS.

Figure 10. ACF, PACF, ACFR and PACFR for water flow at Conhanein the DS and WS.

Figure 11. ACF, PACF, ACFR and PACFR for water flow at Nico in the DS and WS.

Figure 12. ACF, PACF, ACFR and PACFR for water flow at Muianga in the DS and WS.

Figure 13. ACF, PACF, ACFR and PACFR for water flow at Marrambajane in the DS and WS.

Figure 14. ARIMA test analysis for sediment discharge at Intake and Offtake stations.

(WS). Also, significance was found at lag 1 and 2 levels for FIPAG, Massawasse and Muianga (WS), and for Lionde (DS). Additionally, the autocorrelation was significant at 5% at lag 1 level for the Node (WS), Massawasse, Nico and Muianga (DS), Conhane (both seasons).

3.4. ARIMA for Sediment Flux

After the ARIMA test for the sediment discharges, it was found that at the Intake, for DS and WS, sediments followed well the ARIMA model gave good results for the sediments, and indicated good fit between the observed and the predicted ARIMA model data. The adopted structure and parameters of ARIMA model for sediment discharge at CIS based on AIC has a good fit for AR (*p* = 1), whereby, at the Intake the ARIMA *p*-value was 0.822 and 0.932, for WS and DS, respectively. Whilst for the Offtake, the ARIMA *p*-value was 0.877 and 0.893, respectively for WS and DS. These parameters were used to assess the ARIMA model structures and *p,* *d,* and *q* variables calculation, by trial and error procedures, as shown in Figure 14 and Table 5.

The Auto-correlogram (ACF) and Partial Auto-correlogram (PACF) and Auto-correlogram Residuals (ACFR) and Partial Auto-correlogram Residuals (PACFR), all of them before the identification models, are presented. Once the models were identified, the next step produced new models and, therefore, new ACF, PACF, ACFR and PACFR, respectively for each sampling station during both seasons. The graphs for ACF, PACF, ACFR and PACFR for the ARIMA model for sediment discharge, at confidence interval of 95%, are presented in Figure 15 and Figure 16. The auto correlation analysis for the period 2004-2019 are presented in Figure 15 and Figure 16 for the annual sediment flux, in the Intake and Offtake sampling stations, respectively, in the dry and wet seasons. From the figures, can be seen that the autocorrelation is significant at 5% at lag 1, 2 and 3 levels for the Intake (WS), and significance at lag 1 and 2 levels for the Intake (DS). Moreover, the autocorrelation was significant at 5% at lag 1, 2, 5

Figure 15. ACF, PACF, ACFR and PACFR for sediment flux at Intake in the DS and WS.

Figure 16. ACF, PACF, ACFR and PACFR for sediment flux at the Offtake in the DS and WS.

Table 5. ARIMA structure and parameters of model for sediment discharge based on AIC.

Test performed at confidence intervals of 95%; Optimization at Likelihood (Convergence = 0.00001/Iterations = 500).

and 6 levels for the Offtake station (WS). A similar result was observed in terms of significant at 5% at lag 1, 2, 3, 4, 5, 6, 7 and 8 levels in the Offtake (DS).

4. Conclusion

Water flow and sediment forecast have been dealt with in this work for the CIS from the upstream to downstream between 2004 and 2019 period. Positive trend forecast was found in DS and WS. Mann-Kendall tests allowed identifying main trends. All stations showed higher magnitudes and peaks of water flow and sediment flux during wet season over dry season, for every year. Water discharge was found to decrease as it flows downstream of the canal. ARIMA model appeared with good fit between observed and predicted values. These findings will be relevant for irrigation scheme management where water flow and sediment flux constitutes challenge.

Acknowledgements

The authors are grateful to the Hidráulica de Chókwè, Empresa Pública (HICEP) for allowing access to the Chókwè Irrigation Scheme and to the relevant historical data set. We also acknowledge Q-Point, B.V., for encouragement and support in undertaking this research work.

Funding

This research was funded by Instituto Superior Politécnico de Gaza (ISPG) and the Netherlands Initiative for Capacity Development in Higher Education (NICHE) through the NICHE/MOZ/150 project.

References

[1] Zemzami, M. (2016) Flow Forecasts in a Non-Perennial River of an Arid Basin Using Neural Networks. Journal of Applied Water Engineering and Research, 4, 92-101.

https://doi.org/10.1080/23249676.2015.1072851

[2] Frémion, F., Bordas, F., Mourier, B., Lenain, J.-F., Kestens, T. and Courtin-Nomade, A. (2016) Influence of Dams on Sediment Continuity: A Study Case of a Natural Metallic Contamination. Science of the Total Environment, 547, 282-294.

https://doi.org/10.1016/j.scitotenv.2016.01.023

[3] Xia, X.H., Dong, J.W., Wang, M.H., Xie, H., Xia, N., Li, H.S., et al. (2016) Effect of Water-Sediment Regulation of the Xiaolangdi Reservoir on the Concentrations, Characteristics, and Fluxes of Suspended Sediment and Organic Carbon in the Yellow River. Science of the Total Environment, 571, 487-497.

https://doi.org/10.1016/j.scitotenv.2016.07.015

[4] Pektas, A.O. and Cigizoglu, H.K. (2017) Long-Range Forecasting of Suspended Sediment. Hydrological Sciences Journal, 62, 2415-2425.

https://doi.org/10.1080/02626667.2017.1383607

[5] Zhang, X.T., et al. (2015) Bioavailability of Pyrene Associated with Suspended Sediment of Different Grain Sizes to Daphnia Magna as Investigated by Passive Dosing Devices. Environmental Science & Technology, 49, 10127-10135.

https://doi.org/10.1021/acs.est.5b02045

[6] Yadav M.ISH, S.M. and Samtani M.ISH, B.K. (2009) Estimation of Bedload by Kalinske’s Equation Using Tapi River Data. ISH Journal of Hydraulic Engineering, 15, 97-104.

https://doi.org/10.1080/09715010.2009.10514963

[7] Myrhaug, D. and Holmedal, L.E. (2001) Bedload Sediment Transport Rate by Nonlinear Random Waves. Coastal Engineering Journal, 43, 133-142.

https://doi.org/10.1142/S0578563401000311

[8] Sasal, M., Kashyap, S., Rennie, C.D. and Nistor, I. (2009) Artificial Neural Network for Bedload Estimation in Alluvial Rivers. Journal of Hydraulic Research, 47, 223-232.

https://doi.org/10.3826/jhr.2009.3183

[9] Wu, F.-C. and Shen, H.W. (1999) First-Order Estimation of Stochastic Parameters of a Sediment Transport Model. Journal of Hydraulic Research, 37, 213-227.

https://doi.org/10.1080/00221689909498307

[10] Rovira, A., Alcaraz, C. and Ibánez, C. (2012) Spatial and Temporal Dynamics of Suspended Load At-a-Cross-Section: The Lowermost Ebro River (Catalonia, Spain). Water Research, 46, 3671-3681.

https://doi.org/10.1016/j.watres.2012.04.014

[11] Clayton, J.A. and Pitlick, J. (2007) Spatial and Temporal Variations in Bed Load Transport Intensity in a Gravel Bed River Bend. Water Resources Research, 43, Article ID: W02426.

https://doi.org/10.1029/2006WR005253

[12] Ugbaje, S.U., Odeh, I.O.A., Bishop, T.F.A. and Li, J.J. (2017) Assessing the Spatio-Temporal Variability of Vegetation Productivity in Africa: Quantifying the Relative Roles of Climate Variability and Human Activities. International Journal of Digital Earth, 10, 879-900.

https://doi.org/10.1080/17538947.2016.1265017

[13] Afzalimehr, H., Riazi, P., Jahadi, M. and Singh, V.P. (2019) Effect of Vegetation Patches on Flow Structures and the Estimation of Friction Factor. ISH Journal of Hydraulic Engineering, 25, 1-11.

https://doi.org/10.1080/09715010.2019.1660920

[14] Kulkarni, A.A. and Nagarajan, R. (2018) Drone Survey Facilitated Weeds Assessment and Impact on Hydraulic Efficiency of Canals. ISH Journal of Hydraulic Engineering, 24, 1-6.

https://doi.org/10.1080/09715010.2018.1520653

[15] Ngoupayou, J.R.N., Dzana, J.G., Kpoumie, A., Ghogomu, R.T., Takounjou, A.F., Braun, J.J. and Ekodeck, G.E. (2016) Present-Day Sediment Dynamics of the Sanaga Catchment (Cameroon): From the Total Suspended Sediment (TSS) to Erosion Balance. Hydrological Sciences Journal, 61, 1080-1093.

https://doi.org/10.1080/02626667.2014.968572

[16] Wu, X., Wen, J., Xiao, Q., You, D., Liu, Q. and Lin, X. (2018) Forward a Spatio-Temporal Trend Surface for Long-Term Ground-Measured Albedo up Scaling over Heterogeneous Land Surface. International Journal of Digital Earth, 11, 470-484.

https://doi.org/10.1080/17538947.2017.1334097

[17] Osman, I.S.E. (2015) Impact of Improved Operation and Maintenance on Cohesive Sediment Transport in Gezira Scheme, Sudan.

http://edepot.wur.nl/353547

[18] Goyal, M.K., Sharma, A. and Katsifarakis, K.L. (2017) Prediction of Flow Rate of Karstic Springs Using Support Vector Machines. Hydrological Sciences Journal, 62, 2175-2186.

https://doi.org/10.1080/02626667.2017.1371847

[19] Gu, L., Dai, B., Zhu, D.Z., Hua, Z., Liu, X., van Duin, B. and Mahmood, K. (2017) Sediment Modelling and Design Optimization for Stormwater Ponds. Canadian Water Resources Journal, 42, 70-87.

https://doi.org/10.1080/07011784.2016.1210542

[20] Gyimah, R.A.A., Gyamfi, C., Anornu, G.K., Karikari, A.Y. and Tsyawo, F.W. (2020) Multivariate Statistical Analysis of Water Quality of the Densu River, Ghana. International Journal of River Basin Management, 18, 1-11.

https://doi.org/10.1080/15715124.2020.1803337

[21] Sutama, N.H. (2010) Mathematical Modelling of Sediment Transport and Its Improvement in Bekasi Irrigation System, West Java, Indonesia. Master of Science, IHE Delft Institute for Water Education, Delft.

[22] Mendez, V.N.J. (1998) Sediment Transport in Irrigation Canals. (Doctor of Philosophy). Wageningen Agricultural University, Balkema, Rotterdam, Netherlands.

[23] Zetterqvist, L. (1991) Statistical Estimation and Interpretation of Trends in Water Quality Time Series. Water Resource Research, 27, 1637-1648.

https://doi.org/10.1029/91WR00478

[24] Bouchard, A. and Haemmerli, J. (1992) Trend Detection in Water Quality Time Series of LRTAP-Quebec Network Lakes. Water Air, and Soil Pollution, 62, 89-110.

https://doi.org/10.1007/BF00478455

[25] Yu, J.J., Qin, X.S., Larsen, O. and Chua, L.H.C. (2014) Comparison between Response Surface Models and Artificial Neural Networks in Hydrologic Forecasting. Journal of Hydrologic Engineering, 19, 473-481.

https://doi.org/10.1007/BF00478455

[26] Mitosek, H.T. (1992) Occurrence of Climate Variability and Change within the Hydrologic Time Series: A Statistical Approach. CP-92-05, International Institute for Applied Systems Analysis, Luxemburg.

[27] Chiew, F.H.S. and McMahon, T.A. (1993) Detection of Trend or Change in Annual Flow of Australian Rivers. International Journal of Climate, 13, 643-653.

https://doi.org/10.1002/joc.3370130605

[28] Burn, D.H. (1994) Hydrologic Effects of Climatic Change in West Central Canada. Journal of Hydrology, 160, 53-70.

https://doi.org/10.1016/0022-1694(94)90033-7

[29] Baria, P.M. and Yadav, S.M. (2019) Investigating Extreme Rainfall Non-stationarity of Upper Tapi River Basin, India. ISH Journal of Hydraulic Engineering, 25, 1-9.

https://doi.org/10.1080/09715010.2019.1627679

[30] Diop, L., Yaseen, Z.M., Bodian, A., Djaman, K. and Brown, L. (2018) Trend Analysis of Streamflow with Different Time Scales: A Case Study of the Upper Senegal River. ISH Journal of Hydraulic Engineering, 24, 105-114.

https://doi.org/10.1080/09715010.2017.1333045

[31] Qazi, N.Q. and Rai, S.P. (2018) Spatio-Temporal Dynamics of Sediment Transport in Lesser Himalayan Catchments, India. Hydrological Sciences Journal, 63, 50-62.

https://doi.org/10.1080/02626667.2017.1410280

[32] Yu, G.L. and Lim, S.-Y. (2003) Modified Manning Formula for Flow in Alluvial Channels with Sand-Beds. Journal of Hydraulic Research, 41, 597-608.

https://doi.org/10.1080/00221680309506892

[33] Mann, H.B. (1945) Nonparametric Tests against Trend. Econometrica, 13, 245-259.

https://doi.org/10.2307/1907187

[34] Kendall, M.G. (1955) Rank Correlation Methods. Griffin, London.

[35] Cox, D.R. and Stuart, A. (1955) Some Quick Sign Tests for Trend in Location and Dispersion. Biometrika, 42, 80-95.

https://doi.org/10.1093/biomet/42.1-2.80

[36] Hirsch, R.M., Slack, J.R. and Smith, R.A. (1982) Techniques of Trend Analysis for Monthly Water Quality Data. Water Resources Research, 18, 107-121.

[37] Hirsch, R.M. and Slack, J.R. (1984) A Non-Parametric Trend Test for Seasonal Data with Serial Dependence. Water Resources Research, 20, 727-732.

https://doi.org/10.1029/WR020i006p00727

[38] Nabipour, N., Dehghani, M., Mosavi, A. and Shamshirband, S. (2020) Short-Term Hydrological Drought Forecasting Based on Different Nature-Inspired Optimization Algorithms Hybridized With Artificial Neural Networks. IEEE Access, 8, 15210-15222.

https://doi.org/10.1109/ACCESS.2020.2964584

[39] Box, G.E.P. and Jenkins, G.M. (1976) Time Series Analysis Forecasting and Control. Holden-Day, San Francisco.

[40] Kahraman, R., Riella, M., Tabor, G.R., Ebrahimi, M., Djordjevic, S. and Kripakaran, P. (2020) Prediction of Flow Around a Sharp-Nosed Bridge Pier: Influence of the Froude Number and Free-Surface Variation on the Flow Field. Journal of Hydraulic Research, 58, 582-593.

https://doi.org/10.1080/00221686.2019.1631223

[41] Box, G.E.P. and Pierce, D.A. (1970) Distribution of Residual Autocorrelation in Autoregressive-Integrated Moving Average Time Series Models. Journal of the American Statistical Association, 65, 1509-1526.

https://doi.org/10.1080/01621459.1970.10481180

[42] Schmalz, B., Zhang, Q., Kuemmerlen, M., Cai, Q., Jahnig, S.C. and Fohrer, N. (2015) Modelling Spatial Distribution of Surface Runoff and Sediment Yield in a Chinese River Basin without Continuous Sediment Monitoring. Hydrological Sciences Journal, 60, 801-824.

https://doi.org/10.1080/02626667.2014.967245

[43] Shibata, R. (1976) Selection of The Order of an Autoregressive Model by Akaike’s Information Criterion. Biometrika, 63, 117-126.

https://doi.org/10.1093/biomet/63.1.117

[44] Khan, M., Sharma, A. and Goyal, M.K. (2019) Assessment of Future Water Provisioning and Sediment Load under Climate and LULC Change Scenarios in a Peninsular River Basin, India. Hydrological Sciences Journal, 64, 405-419.

https://doi.org/10.1080/02626667.2019.1584401

[45] Sati, S.P., Sharma, S., Sundriyal, Y.P., Rawat, D. and Riyal, M. (2020) Geo-Environmental Consequences of Obstructing the Bhagirathi River, Uttarakhand Himalaya, India. Geomatics, Natural Hazards and Risk, 11, 887-905.

https://doi.org/10.1080/19475705.2020.1756464

[46] Montanher, O.C., De Morais Novo, E.M.L. and De Souza Filho, E.E. (2018) Temporal Trend of the Suspended Sediment Transport of the Amazon River (1984-2016). Hydrological Sciences Journal, 63, 1901-1912.

https://doi.org/10.1080/02626667.2018.1546387

[47] Bozdogan, H. (1987) Model Selection and Akaike’s Information Criterion (AIC): The General Theory and Its Analytical Extensions. Psychometrika, 52, 345-370.

https://doi.org/10.1007/BF02294361

[48] Mozambique: Administrative Division (Provinces and Districts)—Population Statistics, Charts and Map. (2020). MOZAMBIQUE: Administrative Division.

https://citypopulation.de/en/mozambique/admin/

[49] World Meteorological Organization (2012) Limpopo River Basin: A Proposal to Improve the Flood Forecasting and Early Warning System.

http://www.wmo.int/pages/prog/hwrp/chy/chy14/documents/ms/Limpopo_Report.pdf

[50] Magombeyi, M.S., Taigbenu, A.E. and Barron, J. (2013) Rural Poverty and Food Insecurity Mapping at District Level for Improved Agricultural Water Management in the Limpopo River Basin. Paper Presented at the CGIAR Challenge Program on Water and Food, Colombo, 54 p.

[51] Hidráulica de Chókwè, Empresa Pública (HICEP) (2012) Perimetro Irrigado do Regadio de Chókwè. Ministry of Agriculture (MINAG) of Mozambique, Mozambique.

[52] Munguambe, E.O., de Sousa, L.S. and Maluana, C. (2013) Comparison of Water Irrigation Efficiency on Distributors 8 (D8) and 10 (D10) on the Chókwè Irrigation System, Gaza Province-Mozambique. Instituto Superior Politécnico de Gaza, Chókwè.

[53] Fathian, F., Dehghan, Z., Bazrkar, M.H. and Eslamian, S. (2016) Trends in Hydrological and Climatic Variables Affected by Four Variations of the Mann-Kendall Approach in Urmia Lake Basin, Iran. Hydrological Sciences Journal, 61, 892-904.

https://doi.org/10.1080/02626667.2014.932911

[54] Yurekli, K. and Ozturk, F. (2003) Stochastic Modeling of Annual Maximum and Minimum Streamflow of Kelkit Stream. Water International, 28, 433-441.

https://doi.org/10.1080/02508060308691721

[55] Aydin, D. and Memmedli, M. (2012) Optimum Smoothing Parameter Selection for Penalized Least Squares in Form of Linear Mixed Effect Models. Optimization, 61, 459-476.

https://doi.org/10.1080/02331934.2011.574698

[56] Fujikoshi, Y., Yanagihara, H. and Wakaki, H. (2005) Bias Corrections of Some Criteria for Selecting Multivariate Linear Models in a General Nonnormal Case. American Journal of Mathematical and Management Sciences, 25, 221-258.

https://doi.org/10.1080/01966324.2005.10737651

[57] Bozdogan, H. (2000) Akaike’s Information Criterion and Recent Developments in Information Complexity. Journal of Mathematical Psychology, 44, 62-91.

https://doi.org/10.1006/jmps.1999.1277