## Abstract

The European Union Water Framework Directive obliges each country to monitor the groundwater level as it is an important source of drinking water, but also an important part of agriculture. A water budget is used for assessing the accuracy of the groundwater level determination. The computations of the water budget are based on evapotranspiration and the state of land surface hydrosphere. On the basis of the determined water budget, statistics and the prognosis for the next 12 months can be computed. In this paper, all the components of the water budget, such as precipitation, surface run-off and evapotranspiration, are studied for the three tested locations in Poland: Suwalki, Zegrzynski and Tarnow cells. The resultant water budget was also determined and presented graphically. On the basis of the water budget research, a prognosis was determined using AutoRegressive Integrated Moving Average (ARIMA) models with the parameters (2,0,2). A comparison between actual water budget data and a prediction prepared for 2015.08–2016.08 indicated that analysing a 12-month period provides a satisfactory prediction assessment.

## INTRODUCTION

A groundwater system tends to reach equilibrium, at least in the long-term perspective. The aim of this study was to evaluate a water budget derived from NOAH (National Centres for Environmental Prediction/Oregon State University/Air Force/Hydrologic Research Lab Model), which is the closest to *in-situ* measurements among all GLDAS (Global Land Data Assimilation System) sub-models (Pennemann *et al.* 2016), and analyse trends with ARIMA (AutoRegressive Integrated Moving Average) models. The final result is the research prepared in the Suwalki, Tarnow and Zegrzynski cells. The water budget is the quantity of water inflows and outflows, averaged throughout a period of time. Water stored in the global system is constant, but locally some variations can occur, depending on different environmental conditions. Changing climate, droughts and flooding put pressure on us to estimate the global water cycle budget for the purpose of getting to know how much water is on the move and where it is going (Rodell *et al.* 2015). For the purpose of water budget computation, the GLDAS model was used (Matrai *et al.* 2011).

GLDAS is prepared on the basis of merging different simulation and prediction models with observed data over a long-term period, presented in the form of a gridded meteorological dataset on a global scale (Chen *et al.* 2013; Ji *et al.* 2015). The purpose of GLDAS is to integrate observations from satellite and ground measurement systems and provide a land information system. The result is a set of estimated parameters, such as: water storage, energy storage and fluxes. Post-processing is carried out to limit any bias caused by the atmosphere and constraints that are a result of the state of the land surface (Koster *et al.* 2004; Rodell *et al.* 2004; Kumar *et al.* 2006).

The GLDAS model is globally distributed from January 1979, in the form of 28 different fields, which constitute the geophysical parameters. The model is published in two versions: version one consists of four sub-models (land surface models) – VIC (Variable Infiltration Capacity Model), NOAH, Mosaic Land Surface Model and CLM (Common Land Surface model), while version two contains only the NOAH model. It contains land surface states (e.g. soil moisture and temperature) and fluxes (e.g. evapotranspiration). The concept of the GLDAS is to distribute energy and water budget components (Rodell *et al.* 2004; Zaitchik *et al.* 2010).

An essential part of groundwater monitoring and assessment is water budget estimation. The water budget can be understood as the quantity of water flux in the tested area. This needs to be calculated for a specific period of time and in a time scale of dry-dry or wet-wet. The purpose of water budget determination is to indicate the quantity of water flow, which can be extremely important in water monitoring assessment, especially groundwater monitoring.

One of the methods of statistical analysis of the observed data is using ARIMA models. This method allows us to make a prognosis using previous data to perform a simulation of the phenomenon. The ARIMA model is the most effective method for a time series prognosis and distinguishes the following parameters: autoregressive parameter (p), differentiate order (d) and moving average (MA) parameter (q). This consists of the following steps of analysis: identification, estimation and diagnosis. During identification, the number and type of parameters should be identified. A minimum of fifty observations have to be taken into account. Afterwards, all of the parameters are estimated. During estimation, the least squares method is used and the prognosis model is chosen during diagnosis.

The purpose of this paper is to evaluate the water budget in Poland in chosen locations, in equal intervals: month-by-month and in one degree spatial resolution. The second part is to analyse trends using ARIMA models and assess prediction possibilities.

Nearly the entire territory of Poland is covered by the data. It was decided to choose the following areas for the purpose of the research, presented in Figure 1:

near Suwalki (54°06′04″N 22°55′57″E), in this paper called Suwalki cell (SUW), northern Poland;

near Tarnow (50°01′45″N 20°58′18″E), in this paper called Tarnow cell (TAR), southern Poland;

near Zegrzynski Bay (52°27′01″N 21°02′29″E), in this paper called Zegrzynski cell (ZEG), central Poland.

## METHODS

The water balance tends to reach values equal to zero or near zero (i.e. NORMAL water balance). However, in reality, the trend is described as BELOW NORMAL (inflow of water is lower than discharge) or ABOVE NORMAL (inflow of water is larger than discharge) (Smith & Ummerow 2012). For the purpose of water cycle observation, the water budget is computed in the form (Rodell *et al.* 2015):

All elements, i.e. precipitation, surface run-off and evapotranspiration, are natural processes, taken directly from the GLDAS repository. It has been observed that precipitation has the largest influence on the final water budget assessment (Birylo *et al.* 2016). This categorises all water input into a hydrological cycle. When analysing precipitation, an annual trend is observed in Poland. Although the annual trend is calculated in the water budget, for proper assessment, precipitation has to be analysed month-by-month (as the mean value of the annual water budget is not accurate enough). Evapotranspiration is the effect of vaporisation from the land surface and plants. As such, it depends on the type of soil and vegetation. When considering surface run-off, both the slope and neighbouring large-scale basins need to be taken into consideration, as this is essential for the run-off calculation within the tested area (Michalska & Wachta 2015; Bianchi *et al.* 2016).

In this paper, for the purpose of time series analysis, the ARIMA method was used. This method may be used when the time series can be characterised by trend and seasonality. A time series is a statistical way of representing an observation set from equal intervals of time. The aim of a water budget time series analysis is to prepare a prognosis for a future water balance, which can help in making strategic decisions. The first part of the analysis is to separate some elements in the time series, e.g. equal average phenomenon level, trend, long-term cycles, short-term fluctuations, interventions and random noise. The aim of research using ARIMA models is to prepare a prognosis of the further phenomenon. Moreover, it becomes possible to answer questions about the reason for the mechanism causing changes in the water cycle. The prognosis model is dependent upon the time series elements. If a time series consists of systematic elements (equal level and random noise), only the fluctuation coefficient is computed. If the coefficient is lower than 0.05, MA or the so-called naive method is used for smoothing. If only a trend is used for time series identification, for processing, the Holt model or Winters model is used. Moreover, a time series with a trend can be smoothed by linearization. If seasonal fluctuations are observed, Fourier harmonic analysis and the autocorrelation function need to be used. If all systematic elements are distinguished, they should be selected so that decomposition needs to be prepared (Sokołowski 2010; Chen *et al.* 2013; Yan & Ma 2016).

## RESULTS AND DISCUSSION

For the locations mentioned, each parameter needed for water budget computation is plotted in Figures 2 and 3. Basic statistics are summarized in Tables 1–4. It can be stated that for the studied areas, total snowfall and surface run-off had the least influence on the final result. The values of rainfall and evapotranspiration are ten times larger than the surface run-off. On the other hand, in the tested areas, snowfall is five times smaller than evapotranspiration and rainfall. However, as total precipitation is the sum of snowfall and rainfall, it can be concluded that precipitation is a crucial part of water budget assessment. Analysing all plots, we can state that the rainfall and evapotranspiration in the tested locations are almost the same. However, the climate and natural conditions of these places are rather different, with a different influence of the climatic zones and thus different vegetation period ranges. It is very surprising that surface run-off has the largest value in SUW. As expected, total snowfall has the highest value in SUW. Rainfall oscillates from approximately 0^{*}10^{−5} kg/m^{2}/s up to approximately 8^{*}10^{−5} kg/m^{2}/s. The lowest values are observed in June, while most of the rainfall occurs in March. This trend is similar for all tested epochs. From 1979 to 1989, the rate of snowfall (in winter months only) was up to 4^{*}10^{−5} kg/m^{2}/s, and from 1989, we can notice snowfall in a range of 0–2^{*}10^{−6} kg/m^{2}/s, and this is the same in all three locations. Total evapotranspiration is almost the same for the three locations, and the values are apparently periodic, having strong annual variations. Evapotranspiration is the most similar parameter for the tested locations. Minimum values are reached when vegetables acquire more water than they give back by evapotranspiration. When taking into consideration the values of surface run-off, one can see no consistency between the studied locations. All peaks are observed at the beginning of the year. In the years 2000–2009, surface run-off was from 0 to 2.5^{*}10^{−6} kg/m^{2}/s. The highest peaks were observed in 1985–1997 and these values reached 2.5^{*}10^{−5} kg/m^{2}/s.

Cell . | Mean [kg/m^{2}/s]
. | St. Dev. [kg/m^{2}/s]
. | Minimum [kg/m^{2}/s]
. | Maximum [kg/m^{2}/s]
. |
---|---|---|---|---|

TAR | 1.631533e-05 | 1.575977e-05 | −6.99e-06 | 4.979e-05 |

ZEG | 1.554403e-05 | 1.528896e-05 | −1.087e-05 | 5.07e-05 |

SUW | 1.674843e-05 | 1.689836e-05 | −1.022e-05 | 5.405e-05 |

Cell . | Mean [kg/m^{2}/s]
. | St. Dev. [kg/m^{2}/s]
. | Minimum [kg/m^{2}/s]
. | Maximum [kg/m^{2}/s]
. |
---|---|---|---|---|

TAR | 1.631533e-05 | 1.575977e-05 | −6.99e-06 | 4.979e-05 |

ZEG | 1.554403e-05 | 1.528896e-05 | −1.087e-05 | 5.07e-05 |

SUW | 1.674843e-05 | 1.689836e-05 | −1.022e-05 | 5.405e-05 |

Cell . | Mean [kg/m^{2}/s]
. | St. Dev. [kg/m^{2}/s]
. | Minimum [kg/m^{2}/s]
. | Maximum [kg/m^{2}/s]
. |
---|---|---|---|---|

TAR | 5.365487e-07 | 1.472229e-06 | 0 | 1.024e-05 |

ZEG | 5.062832e-07 | 1.31413e-06 | 0 | 1.032e-05 |

SUW | 8.122345e-07 | 2.12159e-06 | 0 | 2.449e-05 |

Cell . | Mean [kg/m^{2}/s]
. | St. Dev. [kg/m^{2}/s]
. | Minimum [kg/m^{2}/s]
. | Maximum [kg/m^{2}/s]
. |
---|---|---|---|---|

TAR | 5.365487e-07 | 1.472229e-06 | 0 | 1.024e-05 |

ZEG | 5.062832e-07 | 1.31413e-06 | 0 | 1.032e-05 |

SUW | 8.122345e-07 | 2.12159e-06 | 0 | 2.449e-05 |

Cell . | Mean [kg/m^{2}/s]
. | St. Dev. [kg/m^{2}/s]
. | Minimum [kg/m^{2}/s]
. | Maximum [kg/m^{2}/s]
. |
---|---|---|---|---|

TAR | 1.938341e-06 | 4.710764e-06 | 0 | 4.829e-05 |

ZEG | 1.635796e-06 | 3.417389e-06 | 0 | 3.15e-05 |

SUW | 2.912832e-06 | 5.443264e-06 | 0 | 4.194e-05 |

Cell . | Mean [kg/m^{2}/s]
. | St. Dev. [kg/m^{2}/s]
. | Minimum [kg/m^{2}/s]
. | Maximum [kg/m^{2}/s]
. |
---|---|---|---|---|

TAR | 1.938341e-06 | 4.710764e-06 | 0 | 4.829e-05 |

ZEG | 1.635796e-06 | 3.417389e-06 | 0 | 3.15e-05 |

SUW | 2.912832e-06 | 5.443264e-06 | 0 | 4.194e-05 |

Cell . | Mean [kg/m^{2}/s]
. | St. Dev. [kg/m^{2}/s]
. | Minimum [kg/m^{2}/s]
. | Maximum [kg/m^{2}/s]
. |
---|---|---|---|---|

TAR | 2.003794e-05 | 1.327048e-05 | 0 | 8.224e-05 |

ZEG | 1.718148e-05 | 1.13828e-05 | 0 | 9.579e-05 |

SUW | 2.0775e-05 | 1.269313e-05 | 0 | 7.033e-05 |

Cell . | Mean [kg/m^{2}/s]
. | St. Dev. [kg/m^{2}/s]
. | Minimum [kg/m^{2}/s]
. | Maximum [kg/m^{2}/s]
. |
---|---|---|---|---|

TAR | 2.003794e-05 | 1.327048e-05 | 0 | 8.224e-05 |

ZEG | 1.718148e-05 | 1.13828e-05 | 0 | 9.579e-05 |

SUW | 2.0775e-05 | 1.269313e-05 | 0 | 7.033e-05 |

Using the above-mentioned parameters and Equation (1), the water budget was computed (Figure 4 and Table 5) for the epochs 1979.01–2016.07. Additionally, a linear trend was computed for all locations. The highest amplitudes were observed in SUW: −4^{*}10^{−5}–6^{*}10^{−5} kg/m^{2}/s. It needs to be stated that water budget changes are very similar in each tested location, but the smallest peaks are observed in TAR. In the years 1979–1995, maxima were noted in SUW, while minima were noted in ZEG. Since 1995, the situation has become very changeable between the three locations, and none of them can be distinguished, as the peaks vary every season. The trend has a downward characteristic (Figure 4). It should be noted that for SUW and TAR, the trend lines are parallel, and the changes are rather large (from 1.25^{*}10^{−5} kg/m^{2}/s in 1979 to 0 kg/m^{2}/s in 2015). In ZEG, the respective changes are smaller (from 5.0^{*}10^{−6} kg/m^{2}/s in 1979 to 0 kg/m^{2}/s in 2015).

Cell . | Mean [kg/m^{2}/s]
. | St. Dev. [kg/m^{2}/s]
. | Minimum [kg/m^{2}/s]
. | Maximum [kg/m^{2}/s]
. |
---|---|---|---|---|

TAR | 5.124403e-06 | 1.491108e-05 | −3.967e-05 | 5.147e-05 |

ZEG | 2.766969e-06 | 1.499913e-05 | −3.686e-05 | 5.611e-05 |

SUW | 6.127168e-06 | 1.751155e-05 | −4.258e-05 | 5.824e-05 |

Cell . | Mean [kg/m^{2}/s]
. | St. Dev. [kg/m^{2}/s]
. | Minimum [kg/m^{2}/s]
. | Maximum [kg/m^{2}/s]
. |
---|---|---|---|---|

TAR | 5.124403e-06 | 1.491108e-05 | −3.967e-05 | 5.147e-05 |

ZEG | 2.766969e-06 | 1.499913e-05 | −3.686e-05 | 5.611e-05 |

SUW | 6.127168e-06 | 1.751155e-05 | −4.258e-05 | 5.824e-05 |

Firstly, decomposition was prepared, taking into account the annual signal (presented in Figures 5–7). The seasonal component is found by smoothing the seasonal sub-series (e.g. the series of all January values, all February values etc.). If the data is periodic, the smoothing is effectively replaced by taking the mean, which was the case in the water budget data. The seasonal values are removed and the remainder is smoothed to find the trend. The mean level is removed from the seasonal component and added to the trend component. This process is iterated a few times. The remainder component is the residuals from the seasonal plus trend fit. The grey bars to the right of each panel on Figures 5–7 show the relative scales of the components. Each grey bar represents the same height but the bars vary in size because the plots are on different scales. The large grey bar in the trend panel shows that the variation of the trend component is small compared to the variation in the data on the top panel, whose bar is much shorter. That is the case for all three tested cells.

It can be concluded that decomposition of the signal was very similar in each tested cell (values between −1^{*}10^{−5} – 1^{*}10^{−5} kg/m^{2}/s). When compared to the raw signal, the seasonal element is surprisingly small (values of a raw signal −4^{*}10^{−5} – 4^{*}10^{−5} kg/m^{2}/s); while remainders reached high values, especially in SUW (values between −4^{*}10^{−5} – 4^{*}10^{−5} kg/m^{2}/s in SUW, −2^{*}10^{−5} – 6^{*}10^{−5} kg/m^{2}/s in ZEG, −2^{*}10^{−5} – 2^{*}10^{−5} kg/m^{2}/s in TAR).

It can be seen from Figures 5–7 that the dominant frequency of the data is the annual frequency. It is possible to use periodogram analysis to perform spectral analysis of the data and find if the other frequencies are present in the series. The periodogram values may be computed by means of the Fast Fourier Transform (Shumway & Stoffer 2011). Such an analysis was carried out for all three cells and it could be concluded that only the annual frequency is significant in the data. It was also proven by autocorrelation analysis. The analysis allows us to conclude that original, full signals are strongly correlated (up to |0.5|) for lag equal to 12. It was observed that the signal values were much higher than the significance level. When annual frequency was turned out, the autocorrelations are much smaller (approximately 0.1) and they are under the level of significance. It can be concluded that the annual frequency is the most important.

ARIMA models, together with exponential smoothing and structural models, are one of the most often used methods of forecasting future behaviour of time series (Valipour *et al.* 2013). It is well-known that in an autoregressive model, the variable of interest is predicted using a linear combination of past values of the variable (AR part of the model). It is called an AR(p) model if p past values of the time series were used for the combination. A MA part uses past errors in a regression-like model. We refer to this as an MA(q) model when errors from q past epochs are used. Additionally, the differencing process may be used (d times), especially for data containing a significant linear trend. The full ARIMA(p,d,q) model is obtained through combining all the above models. To find the values of p, d and q, the R function *auto.arima* (Hyndman & Khandakar 2008) was used. This function computes the best ARIMA model according to indicated parameters. It performs a search over possible models with small values of p, d and q.

The ARIMA calibration period was admitted as January 1979 to July 2015 (440 epochs) and it was used to produce predicted data for the period of August 2015 - August 2016 (12 epochs). Computed ARIMA model parameters are given in Table 6. For the data analysed, the best statistics were obtained for p = 2, d = 0, q = 2, for all three cells. It is equivalent to fitting ARIMA(2,0,2) models to the data. Two coefficients of the AR part of the model (ar1, ar2) and two of the MA part (ma1, ma2) are presented in Table 6. It is worth mentioning that logarithmic transformation of the data was applied before ARIMA model fitting. To avoid logs of negative values, the data had also been shifted by a constant value (5.0^{*}10^{−6} was used).

Cell . | ARIMA parameters . | Coefficient . | |||||
---|---|---|---|---|---|---|---|

ar1 . | ar2 . | ma1 . | ma2 . | intercept . | σ^{2}
. | ||

TAR | (2,0,2) | 1.7315 | −1.0000 | −1.726 | 0.9937 | −9.8457 | 0.06157 |

ZEG | (2,0,2) | 1.7321 | −0.9999 | −1.7376 | 1 | −9.8952 | 0.05153 |

SUW | (2,0,2) | 1.7319 | −1.0000 | −1.7278 | 0.9933 | −9.8488 | 0.06287 |

Cell . | ARIMA parameters . | Coefficient . | |||||
---|---|---|---|---|---|---|---|

ar1 . | ar2 . | ma1 . | ma2 . | intercept . | σ^{2}
. | ||

TAR | (2,0,2) | 1.7315 | −1.0000 | −1.726 | 0.9937 | −9.8457 | 0.06157 |

ZEG | (2,0,2) | 1.7321 | −0.9999 | −1.7376 | 1 | −9.8952 | 0.05153 |

SUW | (2,0,2) | 1.7319 | −1.0000 | −1.7278 | 0.9933 | −9.8488 | 0.06287 |

*t*), previous (

*t*− 1) and one before previous (

*t*− 2).

On the basis of the above results, a prognosis was determined for the three cells. Having models given by equations (Equation (3)) it is possible to forecast future values of the time series. Forecasting was performed using the *forecast* function (Hyndman & Khandakar 2008). The results are presented in Figure 8, with predictions marked in grey. As previously mentioned, the data were logarithmically transformed before the forecasting process (see the ranges in Figure 8).

For the purpose of accuracy assessment, a comparison between actual water budget data and prediction was prepared for 2015.08–2016.08 (Figure 9). It should be emphasised that the epochs taken for the comparison had not been used for the ARIMA model determination. It can be stated that analysing a 12-month period provides a satisfactory prediction assessment. The shape of the prediction chart is less detailed than the actual data, but it covers the water budget curve.

The actual and predicted values are highly correlated, the correlation coefficients are equal to 0.78 (TAR), 0.64 (ZEG) and 0.81 (SUW). Mean differences between the compared values are −2.7^{*}10^{−6}, 1.3^{*}10^{−6} and 5.1^{*}10^{−7} kg/m^{2}/s for the three cells. It can be concluded that using computed ARIMA models of the data allows us to predict the future curves of the water budgets in Tarnow, Zegrzynski and Suwalki.

## CONCLUSIONS

The EU Water Framework Directive obliges each European country to monitor and manage groundwater storage (Kanakoudis & Tsitsifli 2010). For the purpose of assessment of groundwater, the water budget was computed, and then further states were predicted with ARIMA models. The tests were performed in three locations: SUW, TAR and ZEG. For the purpose of water budget computation, the following parameters are required: precipitation (which is the sum of total rainfall and total snowfall), evapotranspiration and surface run-off. All of the parameters were analysed in the mentioned locations for the epochs: 1979.01–2016.07. Afterwards, the water budget was computed. The greatest influence on the final water budget values are snowfall and rainfall (precipitation). The water budget values are in the range of −4^{*}10^{−5} kg/m^{2}/s to +4^{*}10^{−5} kg/m^{2}/s in the tested area.

A comparison between the real water budget data and the prediction for 12 months for the period of 2015.08–2016.08 proved the assumption that the prediction using ARIMA models can give correct results – charts of predicted and actual data show a good agreement, at least in general behaviour. Quantitative evidence of this is high cross-correlation values between predicted and actual data. It should be recalled here, that the actual data used for comparisons were not used for the ARIMA models calibration. As in seasonal decomposition, there are higher frequencies buried in the actual signals, thus ARIMA models of higher orders could also be tested in future.

## ACKNOWLEDGEMENTS

This work was supported by the National Science Centre, agreement number: UMO-2015/17/B/ST10/03927 from March 2016.