O Boulariah, PA Mikhailov, A Longobardi, AN Elizariev, SG Aksenov

1 Department of Civil Engineering, University of Salerno, Via Giovanni Paolo II, 123, 84084, Italy.

2 Ufa State Aviation Technical University, 450000, Russian Federation.

Abstract: Stochastic modelling of hydrological time series with insufficient length and data gaps is a serious challenge since these problems significantly affect the reliability of statistical models predicting and forecasting skills. In this paper, we proposed a method for searching the seasonal autoregressive integrated moving average (SARIMA) model parameters to predict the behavior of groundwater time series affected by the issues mentioned. Based on the analysis of statistical indices, 8 stations among 44 available within the Campania region (Italy) have been selected as the highest quality measurements. Different SARIMA models, with different autoregressive, moving average and differentiation orders had been used.By reviewing the criteria used to determine the consistency and goodness-of-fit of the model, it is revealed that the model with specific combination of parameters, SARIMA (0,1,3) (0,1,2) 12, has a high R2 value,larger than 92%, for each of the 8 selected stations. The same model has also good performances for what concern the forecasting skills, with an average NSE of about 96%. Therefore, this study has the potential to provide a new horizon for the simulation and reconstruction of groundwater time series within the investigated area.

Keywords: Groundwater level forecast; Stochastic modelling; Southern Italy; Seasonality; Homogeneity

Introduction

Groundwater evaluation has a profound impact on the management of water supplies, as groundwater is a major source of drinking water as well as irrigation water (Semiromi and Koch, 2019). The historical time series of groundwater level are the necessary tools for the study of groundwater drought features, to establish a relationship with the meteorological input and to potentially quantify the impact of climate conditions on the groundwater level by Longobardi and Van (2018).Depending on the temporal extension, lengths and frequency of data gaps of the time series,availability of neighbouring data at the regional scale and the characteristics of the study area,different techniques can be used to analyze groundwater levels. These techniques vary from numerical to mathematical, conceptual, and physical models, relying on the information used to explain the dynamic behaviour of aquifers and associated hydrological and hydrogeological variables. The Autoregressive Integrated Moving Average (ARIMA) family of models have been extensively used for the analysis of a wide range of water related issues (Longobardi and Villani, 2006;Faruk, 2010; Boulariah et al. 2019) including the groundwater level modelling (Ahn and Salas,1997; Mirzavand and Ghazavi, 2014; Oikonomou et al. 2018). The complexity and heterogeneity of aquifer behaviour and parameters along with the poor input data requirements, makes the flexible statistical modelling a very attractive tool for the purpose (Suryanarayani et al. 2014; Oikonomou et al. 2018). Considering the detrending and the seasonality, the Box-Jenkins model is the most common approach to forecasting time series(Young, 1999; Adamowski, 2008). It has high accuracy in short-term forecasts (Takafuji et al.2019) and allows forecasting in the conditions of data shortage (Boulariah et al. 2019). In Kashan plain, Isfahan province of Iran, Mirzavand et al.(2015) used time series models and proved that the autoregressive moving average (ARMA) model can predict groundwater level change with high accuracy. In order to simulate the groundwater level in a coastal aquifer, Yang et al. (2017)explored three time series analysis approaches, i.e.,Holt-Winters (HW), integrated time series (ITS)and seasonal autoregressive integrated moving average (SARIMA), showing that all three models can reliably predict the water table, but the SARIMA model showed highest reliability compared to the other two. Takafuji et al. (2019) compared geostatistical simulation (SGS) and the ARIMA model to forecast the depth of groundwater table in the Bauru Aquifer System (BAS), the results showed that the ARIMA model has better performance than the geostatistical model with higher accuracy and precision. In Italy, Nunno and Granata (2020)used the NARX neural network to predict groundwater level in Apulia region. They reported that very good results were obtained for the region of Apulia and suggested to use the NARX network in other areas characterized by Mediterranean climate and karst phenomena for groundwater prediction.

In this direction, the current paper aims to address the problem of predicting and forecasting groundwater table fluctuations in a specific area of the Mediterranean region, the Campania region, in southern Italy, characterized by temporally discontinuous observed time series, by using a statistical approach based on SARIMA tools.

1 Methodology

1.1 Time series homogeneity tests

Different statistical tests were used for homogeneity and trend analysis on time series data.Pettit’s test and the student’s test have been used in this study.

1.1.1 Pettitt's Test

In some cases, where the exact time of displacement is unknown, the use of the Pettitt’ s test allows us to determine the moment of change in the time series. So, the Pettit’s test is based on the following rules:

Test statisticsKand the corresponding confidence level ρ for the sample (n) can be described as:

When a certain value is more than the confidence levelρ

The null hypothesis will be rejected. For the change point, the approximate probability of significance p is determined as follows:

Thus, in the place where the change point is located, the series is divided into two parts (Jaiswal et al. 2015).

1.1.2 Student's Test

The t-test refers to a parametric test in which the difference between the theoretical value and the average value is indicated by an alternative hypothesis, while the null hypothesis indicates that these values are equal to each other. This test has more statistical power than the nonparametric tests.t criterion is calculated by the formula (Cleophas et al. 2016):

Where: µ is the specified value;x¯ is the mean of the samples;sis the standard deviation of the samples;nis the size of the samples.

These methods of testing the hypothesis on data homogeneity are selected on the basis of accuracy as well as the simplicity of utilization of results.Data analysis was performed by using the statistical software R studio.

1.1.3 Time Series prediction and forecasting models

The fundamental steps in forecasting and time series modeling, as described by (Valipour et al.2012), are to plot the time series and determine its basic characteristics. Regarding various features such as trend and seasonality, it is important to look for potential signs that the time series has changed over the investigation period. Next, trend or seasonality should be eliminated by differentiating or by fitting model. The main goal of this process is to produce a set of stationary residuals (Valipour et al. 2013). As a final step, a forecasting model for the residuals should be developed. Occasionally, several potential models can be found, and in this case, an additional analysis must be performed to determine the best model for deployment.

Sometimes, potential models can be eliminated based on how well they fit over historical data. It is possible that a poorly adapted model may produce good forecasts. Montgomery et al. (2015) has shown that it is necessary to validate the performance of a model from the previous step. This will likely involve some type of split sampling or crossvalidation procedure. The objective of this step is to select a suitable forecast model.

In addition, it is interesting to know the differences between the original data and the forecasted values by the original scale model. The methodology for determining the optimal seasonal ARIMA model consists of six consecutive steps as described above. Fig. 1 demonstrates the step-bystep process to forecast time series from data analysis to the estimation of the corresponding model.

Fig. 1 Flowchart of applied methodology for forecasting time series

The use of ARMA models has proven itself in the forecasting and analysis of hydrological time series, which is confirmed by Cleophas and Zwinderman (2016), Balasmeh et al. (2019) and Mombeni et al. (2013). As mentioned above, the ARMA model is a combination of the autoregressive (AR) model and the moving average(MA) model (Kumar and Rathnam, 2019). AR part indicates that the variable of interest regresses according to its previous values. MA part shows that a linear combination of a limited number of residuals is a regression error. In other words, the deviation of a dependent variable from its average value is a linear combination of current and past values of the random perturbation vector. Thus, the ARMA model allows one to make a forecast that depends on both the current and past values of variable as well as the magnitude of normal disturbances of the current and past values. The stationarity of data usage is an essential criterion for the use of such models.

The chosen model is based on manually selected parameters that depends on the values of ACF and PACF graphs, which minimizes the Akaike and Swartz criteria, and is also based on essential coefficients. The model will not be accepted if the resulting model does not meet certain mathematical conditions (Valipour et al. 2013a). The forecasted values have a predictive interval that increases as the forecast horizon increases. The forecast intervals are based on the theory that there is no autocorrelation of the residues and their normal distribution. If one of the conditions is not satisfied, the forecast may turn out to be incorrect.Therefore, A graph and residual histogram are built before the forecast, and together with the Lung-Box test, a visual analysis that allows us to draw conclusions whether the model fits or not.

The Lung-Box test allows us to check the hypothesis that the autocorrelation coefficient differs from zero. The null hypothesis of the test indicates that the remnants of the model are white noise. An alternative hypothesis suggests that the data is not random. The formula for calculating Qstatistics is as follows :lag length;rkisthe sample autocorrelation coefficientforlagk;

Where: T is the sample size; m is the maximum

The critical value is given as follows:

The Lung-Box test is based on the Box-Pierce statistics. Therefore, it has the same asymptotic distribution and gives comparatively similar results for relatively large numbers. The criterion does not lose its consistency even if the process does not have a normal distribution. When the data under study is non-stationary, the box-Jenkins models are used. To reduce the data to a stationary form, a method of taking consecutive differences is used.If the series becomes stationary after these transformations, we can use the ARIMA model to forecast. In this case, the initial step of differentiation can be applied one or several times to eliminate the unsteadiness (Valipour et al.2013b). There are four consecutive stages to determine the Box- Jenkins model:

(1) Model identification-the process of selecting a model that is best suited to the process in question.

(2) Model assessment-the use of regression methods to obtain the parameters included in the model.

(3) Model test-check the adequacy of the model using tests for normality and autocorrelation of residues.

(4) Model application-use of the model focus forecasting. (The design of the technique is presented in Fig. 2).

Fig. 2 ARIMA model prediction method (Chatfield et al.1973)

In the case where the data contain a seasonal component, the SARIMA model is used for forecasting. The difference between SARIMA and ARIMA model is to add three parameters (P, D, Q)to predict data seasonality. The equation of the ARIMA/SARIMA model in general looks like this:

Where: AR(p) is the operator that converges to provide the operator (MA) with stationary conditions, which is a polynomial of the orderqin B, such that each value is the sum of the previous white noise values ofq+1.

MA (q) is the operator which converges to ensure the reversibility

1.1.4 Goodness of fit

In order to analyze the performance of these models, various efficiency indices, such as Nash-Sutcliffe Efficiency (NSE) (Nash and Sutcliffe,1970), R2determination coefficient, root mean square error (RMSE) (Singh et al. 2005), Mean Absolute Error (MAE) (Willmott et al. 2005),Index of Agreement (Willmott, 1981), Mean Squared Derivative Error (Willmott et al. 2005)and the percent bias (%) (Singh et al. 2005). In addition, Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) were used to evaluate the maximum likelihood of models.Moriasi et al. (2007) suggested a rating method for this purpose to classify models into four types:very good, good, satisfactory and unsatisfactory.

2 Study area and datasets

The Campania region in southern Italy is situated between 40.0°-41.5° N and 13.5°-16.0° E,extending from the Apennine Mountains to the Mediterranean Sea with steadily decreasing elevations from the interior to the coastline,covering approximately 14 000 km2(Fig. 3, right panel). The climate of study area is usually seasonal, with some apparent variations depending on the region. With the greater amount of precipitation observed during the winter months,the seasonality is very marked. The mean annual rainfall ranges from 600 mm to 2 400 mm in the study region (Fig. 3, left panel), while the average annual temperature is about 17 °C.

Fig. 3 Digital elevation model with location of water wells (right panel) and mean annual precipitation (left panel) in the study area

Available data consist of monthly groundwater level time series for over 44 sites located across the Campania region and the Lazio region, the data were reviewed and analyzed in this paper from 1926 to 1994. The water wells database was provided by the national hydrographic and mareograph service SIMN. After a data quality control performed using different statistical tests,only 8 out of these stations were effectively used in the area under investigation.

3 Results and discussion

To determine the parameters of the model, the auto ARIMA function of the R software was used, the function allows us to automatically select the model with the lowest AIC, BIC indices, which significantly reduces the time for the selection of parameters and model order. When evaluating the observed values, a restriction in the enumeration of parameters was used. Models with a step of differentiation (d＞1, D＞1) were excluded. This restriction is based on studies conducted by Hamilton (1994) and Mohammadi et al. (2005),which recommended that the differentiation order for time series should not be more than 1,otherwise the results will be far from reality. The following is an example of estimating the observed values for data from a measuring station in Casamicciola. We applied the SARIMA model(1,1,2) (0,1,1)12, with the lowest AIC and BIC(AIC = -1 226.74 and BIC = -1 207.5). A comparison of the predicted and observed values is presented in Fig. 4.

Fig. 4 Predicted and observed values of the ground water level

Fig. 5 Autocorrelation function graph for the residues of the SARIMA model (1,1,2) (0,1,1)12

As it can be seen from the picture (Fig. 4), the data have an increasing trend, which is confirmed by the tests for stationarity. According to the Pettitt test, the series is heterogeneous and a shift is observed in August 1946. The model describes the data well enough, which is confirmed by the efficiency indices between the estimated and observed values. A high correlation value between graphs indicates that the prediction is highly accurate. The adequacy of the model was verified using the residual analysis shown in Fig. 5.

Residuals look like white noise and correspond to a normal distribution and the values of the autocorrelation lags are within the confidence interval, which indicates the absence of autocorrelation in the residues. The histogram of residuals shows symmetric distribution relative to zero,which indicates the absence of systematic errors.Those results suggest that this model is well suited for evaluating the data in question and can be used for prediction and forecasting. In a similar way,models were found for predicting all 8 stations, the results are presented in Table 1.

The analysis of the resulting models allowed us to determine the three most frequent forecasting models, which are ARIMA (0,1,3) (0,1,2)12,ARIMA (1,1,2) (0,1,1)12, ARIMA (1,1,1) (0,1,1)12.The three recurring models were applied for each of the stations. The model parameters were calibrated using various statistical indices of model quality, such as (NSE) (Nash and Sutcliffe, 1970),correlation coefficient R2, root mean square error(RMSE) (Singh et al. 2005) and percentage bias(%), etc. The goodness-of-fit assessment of the models was carried out on the basis of 11 statistical criteria presented below (Table 2).

Table 1 Forecasting models Results based on AIC and BIC criteria

The analysis of Table 2 shows high values of Pearson correlation for each model and the average value is 0.94. Such a high correlation shows that each of the models describes the data well enough to ensure high accuracy of prediction.

However, the SARIMA model (0,1,3) (0,1,2)12is the best for predicting the groundwater level for each of the 8 stations in the Campania region, since it has the highest correlation between the observed and forecasted values and the lowest errors than two other models. Using this model to forecastgroundwater level behavior for 24 months,the obtained results are shown in Fig. 6, where blue and gray colors highlight the 80% and 95%confidence intervals respectively.

Table 2 The comparison of different SARIMA models based on accuracy measures

Fig. 6 Forecasting using SARIMA(0,1,3)(0,1,2)12; (a) Acera Capomazzo, (b) Casamicciola, (c) Cassano di Sessa Aurunca, (d) Forio(Calitto), (e) Forio(Pontone), (f) Forio(Umberto I), (g) Nocelleto di Carinola, (h) Parete(tre ponti) scrivi cosa sono le fasce colorate

As shown in Fig. 6, this model is very suitable for predicting time series in each region; data is well predicted despite the presence of trend and seasonality. For clarity, and according to classification by Moriasi et al. (2007), Table 3 shows the values of some statistical criteria and their average values when using the SARIMA model (0,1,3) (0,1,2)12.

The average value of the indices for the selectedmodel shows that this model describes the data accurately and allows us to predict and forecast the selected time series with high accuracy and precision.

Table 3 Statistical index values for the selected SARIMA model (0,1,3) (0,1,2)12

4 Conclusions

Ground water level forecasting is an integral tool in the management, implementation, and operation of drinking water supply systems. The use of the Box-Jenkins technique is a reasonably accurate approach for groundwater time series modeling.The results show that SARIMA models predict monthly groundwater levels with high accuracy.The seasonal autoregressive integrated moving average models are sensitive to extreme values and data periodicity. The procedure used to carry out the work simplifies the search for the desired combination of prediction parameters. In fact, the most suitable parameters for predicting the groundwater level behavior of the Campania region is the (0,1,3) (0,1,2)12, which is confirmed by the analysis of statistical indices. The SARIMA models obtained have furthermore a sufficiently high forecasting accuracy although there is basically no way to improve them, since they select a significant part of the information from the data. The relevance of this model for predicting the groundwater level in each of the 8 selected stations proves the possibility of using the model to simulate the situation in the entire region using the data of remaining stations.