Estimating a Reliable Water Budget at A Basin Scale: A Comparison between the Geostatistical and Traditional Methods (Foro River Basin, Central Italy)

.


Introduction
The reliable estimation of the water budget at a river basin scale is crucial for proper water-management practices and the sustainable multi-purpose exploitation of water resources, especially considering the evident and increasing impact of climate change on water availability [1][2][3][4][5][6].However, the quantification of available water resources connected to the direct rainfall recharge, as in river basins [7], is often affected by the problem of the spatial representativeness of the data (i.e., rainfall and temperature) collected in sparse weather stations across areas that are hundreds to thousands of square kilometers, and which are usually recorded as time series and then interpolated to obtain a spatial distribution [8].As a matter of fact, the distribution of weather data is often not optimal according to both a spatial and an altimetric point of view [9].
To overcome this limitation related to the spatial representativeness of discrete measurements from a sparse weather monitoring network, simple spatialization approaches, such as triangulation to obtain isolines or the Thiessen method to obtain representative values within areas of influence of a single station, have traditionally been used [10,11].However, even though they are quite robust, these approaches to the spatialization of rainfall and temperature data and then to the quantitative estimation of the water budget at a basin scale do not provide any reliable estimation of the uncertainty due to the inherent spatial variability of the meteorological phenomena, which cannot be physically measured by a sparse monitoring network.
In recent decades, different approaches have been used to tackle the spatialization issues, like data fusion with satellite-based information on precipitation [12,13], the use of machine learning techniques [14,15], and geostatistical methods, which have been developed and tested in several hydrogeological and environmental applications [16,17].In particular, geostatistical techniques are effective when it comes to studying spatially relevant natural phenomena, as they take advantage of the spatial correlation among measurements over a specific study area to provide reliable estimates of the variables of interest [18][19][20][21], additionally providing the quantification of the associated uncertainty [22].
In the present work, we aimed to test the feasibility of using a stationary geostatistical technique to spatially estimate rainfall and temperatures to provide a reliable estimation of the Foro river regime over an average year.To achieve this purpose, we compared the water budget estimation in terms of river discharge obtained from the spatialization of monthly rainfall and temperature with both the Ordinary Kriging and the Thiessen methods (see the flowchart in Figure 1 for a better comprehension of the theoretical approach).The original datasets were obtained by calculating the monthly average rainfall and temperature from a 33-year-long time series (i.e., from 1986 to 2019), collected by the Hydrographic Service of Abruzzo Region, at each weather station of its monitoring network.
Since the objective of the study is to compare, in the Foro valley test area, the results obtained by applying the two techniques using only and exclusively the real data of the available monitoring network, neither virtual weather station was used in the calculations using the traditional method to cover unmonitored areas (generally corresponding to altitudes greater than 1000 m a.s.l.(above sea level)), nor were geostatistical correlations with altitude performed in the geostatistical estimation process.
As a comparison term, the flow rates of the Foro river monitored in the same statistical time interval in two monitoring stations located approximately in the middle and at the end of the river course were used, which are linked to an 87 km 2 and a 232 km 2 wide sub-basin, respectively.

Study Area
The study area is the Foro alluvial river basin located in the Periadriatic area of the Abruzzo Region (Figure 2), which is mainly characterized by foredeep deposits and alluvial deposits of the main rivers.Figure 3 shows the geological-hydrogeological framework in detail: in the southwestern part, at the highest altitudes, the basin is characterized by calcareous-marly deposits, while in the other portions, it is characterized by Plio-Pleistocene clays with sandygravelly levels, and Quaternary alluvial and continental deposits are present.The Foro basin is about 236 km 2 wide, and from the altimetric point of view, it ranges between the sea level and about 2000 m a.s.l.The most permeable complexes can be found in the SW portion and along the Foro riverbed, while the less permeable ones are observed in the middle part of the basin.
Along the Foro river, the anthropogenic factors are relevant, and a significant amount of water resources is exploited for drinking purposes and for fields' irrigation.Furthermore, a hydroelectric plant can be found within the Foro basin, which takes from the river and then returns it downstream.
This kind of plants can be found in almost all the catchments of the rivers of central Italy with Adriatic drainage [23][24][25] and in the Apennine intra-mountain basins [26].

Meteorological Network
The rainfall and temperature datasets used have been collected by the Hydrographic Service of Abruzzo Region database for a 33-year period, from 1986 to 2019.Table 1 summarizes the selected monitoring stations and their main features.
In the traditional approach, data from just 10 weather stations inside or immediately outside the Foro basin were used (Figure 3), while for the geostatistical analyses, data from the whole Abruzzo and Molise regions were used to estimate the spatial distribution of both rainfall and temperature (Figure 4).Eight out of ten of the stations considered in the traditional approach had both a rainfall gauge and a temperature sensor, while two of them (i.e., ASL and CLD) had just the rainfallmeasurement system.For the latter, monthly and annual temperatures were estimated using linear regression between temperature and elevation.As can be observed in the Altitude column in Table 1 and in Figure 3, weather stations are not homogeneously distributed and cannot be found at altitudes over 1280 m a.s.l., nor between 700 and 1300 m a.s.l.
Two hydrometers, Foro a Molino Galasso and Foro a Ponte di Vacri, were considered inside the Foro basin in order to define two sub-basins and compare measured discharge with calculated ones.

Geostatistical Method
To spatialize both the rainfall and the temperatures measured at the weather station of the regional monitoring network, we used the stationary technique called Ordinary Kriging (OK) [19,20,27].OK estimates the target variable (z *  ) at each location of the selected spatial domain ( ) through an unbiased and optimal estimator called the Best Linear Unbiased Estimator (BLUE), which is defined by the following equation: In this equation, λ represents the weights assigned to the variable measurements (z  ) within a certain distance called a neighborhood ( ).
The unbiased OK estimator imposes the following condition to ensure that the estimated values are the most optimal and unbiased (i.e., E z *  z  0): This condition represents a constraint within the Kriging equation system (Equation ( 3)), which consists of a set of N 1 linear equations: In the OK equation system, μ is a Lagrangian multiplier, while γ  ,  and γ  ,  are the variograms related to pairs of measurements and to pairs of points that include the unsampled location ( ).
The variograms are described by a function that incorporates the spatial dependency of a given random variable of interest and describes the relation between semi-variance (γ  ) and distance in terms of a separation vector or lag ().Variograms are defined by the equation (Equation ( 4)) defined below: where z  and z   are a pair of distinct measurements separated by a lag  at a specific location within the spatial domain ( ), and N  is the number of pairs separated by the lag.
To solve the OK linear equation system (Equation ( 3)), the experimental variogram (obtained from actual measurements) is fitted using a variogram model.
In addition to the predicted value at each target location on the gridded domain, OK allows quantifying the uncertainty associated with the estimate in terms of Kriging variance (σ  ): It is important to highlight that the Kriging variance, as defined in Equation ( 5), and the corresponding standard deviation can be used as a local measure of error only when the variable of interest has a Gaussian statistical distribution, as the prediction may be non-linear and not optimal to overcome this limitation.All the monthly rainfall and temperature data were transformed into standardized variables (i.e., mean equal to 0, and variance equal to 1) through the Gaussian Anamorphosis [22].This function converts a Gaussian variable (Z Φ Y ) into a non-Gaussian one by fitting a polynomial expansion, as defined below: In this equation, H Y are the Hermite polynomials, while Ψ are the Hermite coefficients.
Once the Gaussian Anamorphosis function is defined, it is possible to use its inverted version to transform a non-Gaussian variable into a standardized one (Equation ( 7)) as follows: In this study, all the raw monthly rainfall and temperature data obtained from the 33-year-long time series collected over the whole regional monitoring network were previously transformed into standardized Gaussian variables and then used to fit the variogram model and eventually interpolated across the entire Abruzzo region.Finally, the predictions were back-transformed to obtain the monthly rainfall and temperature distributions within the selected domain through the Gaussian Anamorphosis function.Back transformation was applied to 95% confidence interval limits (Lower Limit-LL, and Upper Limit-UL) maps as well, obtained using the following relation, to provide a quantification of the uncertainty associated with the rainfall and temperature estimates: where σ is the OK standard deviation, whereas n is the optimal number of measurement locations in the neighborhood.The rainfall and temperature values were estimated via OK through Equation (1) on a grid (i.e., support) as large as the entire Abruzzo region, with a cell size of 100 × 100 m, and then cut with a polygon corresponding to the Foro basin to compare the estimates at a catchment scale.

Water Budget Estimation
The water budget is traditionally defined as where P is the total rainfall related to a certain area; ET is real evapotranspiration; and O is outflow, defined as sum of runoff (R) and infiltration (I), which is correlated with the potential infiltration coefficient (IR).In this work, rainfall and temperature were analyzed through two approaches, the traditional method based on the Thiessen polygons and the geostatistical spatial estimation.
For the first method, the ten thermo-pluviometric stations' positions were considered to draw the Thiessen polygons (Figure 3); the corresponding rainfall and temperature data were cumulated to a monthly and annual resolution and then averaged to obtain datasets representative of the whole 33-year-long time series.
In the second approach, the Ordinary Kriging was applied to rainfall and temperature data from gauging stations all over Abruzzo and Molise regions (Figure 4) for the same 33-year period, obtaining the estimated spatial distribution of monthly and annual rainfall and temperature.For each variable, three maps were carried out, one for the estimated values and two for the relative errors identified as Upper (UL) and Lower (LL) Limits.
In both approaches, the real evapotranspiration (ET ) was calculated using the Turc, Turc modified [28], and Thornthwaite and Mather [29] methods; mean real evapotranspiration values related to a statistically significant period (i.e., over at least 30 years) were provided from both methods and can be assumed as representative of the local meteoclimatic condition.
The Turc method provides yearly ET values through the following relation in Equation (10): where L is the evaporative potential of the atmosphere 300 25T 0.05T and T is the mean yearly temperature of air (°C).The Turc modified is also based on Equation (10), but it considers L as defined by 300 25T 0.05T , with T ∑ , and P and T are the rainfall and air temperature values related to the ith month, respectively.This method quantifies evapotranspiration without considering seasonal variation in the total amount of water returned to the atmosphere, either to affect air temperature (i.e., evaporation) or for plant life and growth (i.e., transpiration).
The Thornthwaite and Mather method [29] offers a more accurate estimation of the evapotranspiration by calculating potential evapotranspiration in relation to the ith month ET through an exponential equation (Equation ( 11)): .
Monthly ET values were compared with the residual water content within the shallower portion of the soil, where plant roots influence the water budget, to estimate the monthly evapotranspiration values (ET ).In this way, the yearly ET value was estimated while considering the seasonal variability and the actual availability of water in the topsoil.
After calculating the amount of water returning to the atmosphere, the outflow was calculated according to Equation (9).In order to quantify runoff and infiltration, Potential Infiltration Coefficients (IR) derived from the most complete geological map of the study area [30] were considered.IR values were assigned to every hydrogeological complex [24] according to the predominant lithotype.
In both methods, Equations ( 12) and ( 13) have been used to calculate infiltration (I) and runoff (R), respectively: In order to compare the calculated water budget with the discharge measured using the two hydrometers in Figure 3, two sub-basins were considered, and the relative runoff was estimated, as a reference for the conversion of runoff into river discharge equivalent.

Traditional Water Budget Method
Based on the principles described in Section 2.4, "Water Budget Estimation", Table 2 summarizes the results of the traditional approach to estimate the water budget: monthly and yearly runoff were calculated for both the sub-basins and then converted into m 3 /s to be compared with measured discharge.In general, the traditional water budget shows that about 65% of the inflows return to the atmosphere, while 35% is available for surface runoff and infiltration.
Table 2. Monthly rainfall, evapotranspiration, and runoff values for the two sub-basins (Figure 3) carried out with the traditional method.

Geostatistical Water Budget Estimation
The geostatistical analyses were carried out for the twelve average monthly datasets available following the method explained in Section 2.3, "Geostatistical Method".In   Examples of monthly outflow, infiltration, and runoff maps can be observed in Figure 6, related to Molino a Galasso sub-basin.The water budget estimation has been carried out using the theoretical principles explained in Section 2.4 and applied to the rainfall and temperature geostatistical maps obtained by applying the Ordinary Kriging.
As can be seen, the geostatistical method allows us to obtain a more accurate estimation of the water budget terms all over the study area, with a resolution corresponding to the chosen cell size.As a result of the application of the water budget, the outflow (upper part in Figure 7) appears substantially connected to the distribution of precipitation and temperature, which, in turn, are essentially conditioned by the elevation and the orography.Instead, both runoff and infiltration maps reflect the lithologies of the area, which are directly correlated to the IR coefficient used in Equation ( 12) for the calculation.Accordingly, the two distributions appear complementary to each other.A more intense infiltration and consequently lower runoff can be seen in the SW area of the sub-basin, where the IR is higher (i.e., about 80-90%), whereas infiltration is less intense in the central portion of the sub-basin.
In order to compare runoff values obtained with the geostatistical method and the results in Table 2, zonal statistics have been applied to each map.This GIS tool allowed us to obtain statistical parameters, such as minimum, maximum, and mean values, as well as the sum of each pixel value for every raster map.This approach was applied to both subbasin results.

Comparison between Traditional and Geostatistical Methods
In Figure 8, water budget results from both methods are compared to each other and with direct discharge measures for each sub-basin.
For the Ponte di Vacri station, there is a good correspondence between both the geostatistical and traditional methods.On the other hand, an underestimation is highlighted for the Molino Galasso station, where the traditional method gives discharge values consistently lower than those of OK.In any case, in both sub-basins estimated water budgets have the same trend, with a maximum during the wet season and a zero value during the dry one.This last consideration derives from the fact that Thornthwaite's method, applied to groundwater-dependent areas, such as the one under study, does not consider the modulating effect of the soil and aquifers, which often affects the presence of outflow even during the dry season (a local occurrence of this phenomenon is reported in [31]).In any case, the fact that OK also allowed the estimation of 95% confidence interval limits makes the results obtained with the geostatistical analysis more similar to the ones obtained through the traditional approach: at Ponte di Vacri, values from the traditional approach are almost always within the 95% confidence interval, whereas at Molino Galasso, they are close to the Lower Limit.
The comparison between calculated and measured discharge shows some differences; a shift probably caused by a delay in the natural system (i.e., infiltration and groundwater flow) can be observed.During the dry season, a basal flux is evident in river discharge data, but this is not present in estimated values.This evidence may be connected to the presence of arenaceous and alluvial deposits, which usually host local-to-regional aquifers and provide a constant water supply to the river even during the dry season [32].Furthermore, because of calculation assumptions, only direct contributions (i.e., precipitation and temperatures) were considered in the estimation of the monthly river discharge through both the geostatistical and traditional methods.
Moreover, in the Molino Galasso sub-basin, hydrometer data are influenced by measure availability because only four years were recorded.Despite this short monitoring period, measured discharge from January to April shows good correspondence with the discharge calculated with the geostatistical estimation data, and the shifting observed in Ponte di Vacri is less pronounced.This difference is probably because Ponte di Vacri subbasin (Figure 3) is closer to the calcareous complex in the southwest part of the basin, while the Molino Galasso one is located close to the Adriatic Sea, where the infiltration delay is less evident.2).
The yearly runoff estimations for the two sub-basin and for the whole Foro river basin are shown in Table 3.For each method, the yearly water budget was calculated using the Turc, Turc modified, and Thornthwaite and Mather approaches.
The results obtained from both the traditional and geostatistical methods were compared to each other and with annual measured discharge, such as for monthly data.In this case, hydrometer measures rose by 0.5 m 3 /s, corresponding to the amount of water drawn annually for drinking purposes [33].In the Ponte di Vacri sub-basin, both methods underestimate the measured discharge.This evidence can be explained by considering that, by subtracting discharge estimated with the geostatistical method (about 0.3 m 3 /s) from the measured one (1.1 m 3 /s), a 0.8 m 3 /s of surplus is obtained.Comparing this result with the Molino Galasso one, the difference between measured and estimated is still 0.8 m 3 /s.This comparison thereby supports the hypothesis that 0.8 m 3 /s is an external contribution to the estimated discharge, considering that literature data suggest, for exactly the carbonate aquifer in the southwestern side of the study area, an infiltration rate of 0.029 m 3 /s/km 2 [31,34,35].This infiltration rate, related to the 0.8 m 3 /s external contribution, corresponds to a 27 km 2 of additional area adjacent to the considered catchment.The situation described is very common in the carbonate aquifers of the Apennines [36][37][38].
Also in this case, the values estimated through a traditional approach either fall within the 95% confidence interval or are very close to one of the two interval limits.This evidence suggests that the geostatistical approach provides reliable estimates of the water budget, as it quantifies the uncertainty related to the fact that the measurement of both rainfall and temperature is discrete and the monitoring network is too sparse to be able to effectively describe the spatial variability in meteorologic phenomena at a basin scale.Nevertheless, the data availability is one of the most critical factors for the application of geostatistical techniques, as these methods need an appropriate number of measures.Also, the traditional approach could benefit from a higher number of point data.The recent developments obtained with the use of weather RaDAR data are encouraging [9,39,40] and may represent a valuable additional source of information to be integrated into the water budget estimation, especially through an advanced geostatistical approach (e.g, Multi-Collocated Co-Kriging or Kriging with External Drift).The use of weather Ra-DAR data would allow estimating in a more reliable way the spatial distribution of rainfall on a finer grid mesh and with a lower associated uncertainty.

Conclusions
At first sight, the traditional and geostatistical analyses of input data for the water budget could not be more different.The traditional method is based only on point observations of rainfall and temperature.On the other hand, the geostatistical method is built on spatial variability models (i.e., variograms) and allows taking advantage of the spatial correlation among observations to provide reliable estimates and uncertainty quantification.
The results compared in the graphs in Figure 7 show a similar trend: a most intense discharge during the wet season, with a maximum in December, and a slow decrease from January to May until a zero value is reached during the dry season.The best correspondence between the two methods can be observed in the Ponte di Vacri sub-basin, while in the Molino Galasso one, the discharge calculated with the traditional method is lower than that of the geostatistical ones during the wet season.However, the possibility of also calculating the 95% confidence interval limits with the Ordinary Kriging makes the results obtained with the two considered approaches more similar to each other.In fact, the values estimated through a traditional approach either fall within the 95% confidence interval or are very close to one of the two interval limits, suggesting that this geostatistical technique provides reliable estimates of the water budget.
The comparison between the discharge values calculated with the two methods and the measured one appears more pronounced in the dry season.This is mainly due to the presence of local-to-regional arenaceous and/or alluvial aquifers, which provide a constant water supply to the whole hydro (geo)logical system.These discrepancies between monthly measured discharge and estimated values can also be explained using the water budget calculation method, which does not take into account additional inflows, such as the water subtracted or added by human activities or the contributions from other aquifers.Moreover, the geomorphologic features are not considered in this work, but they can indeed influence infiltration and runoff.Infiltration was assumed as a net loss for the river basin system, but it is a dynamic resource over a 30-year statistical period, especially in groundwater-dependent systems and in the presence of river-aquifer hydraulic connections.
The comparison between annual runoff obtained through estimation and measurement pointed out the presence of an external contribution of 0.8 m 3 /s, which may be related to the local carbonate aquifer in the southwestern side of the study area.In the calculated water budget, this additional inflow is likely related to a volume of water previously lost as infiltration and then returned with a delay.In addition, even though the obtained results are encouraging, it is important to point out that the two methods had to overcome some issues, such as inhomogeneous databases through time, the impact of human activities along the Foro river in terms of water utilization and partial return, and the non-overlapping between the hydrographic and hydrogeological catchments.
In conclusion, the application of the Ordinary Kriging technique to rainfall and temperature measurements proved to provide reliable estimates of the water budget at a basin scale, very similar to the ones that can be obtained using the traditional approach.However, the geostatistical method is additionally able to quantify the uncertainty related to discrete measurements of both rainfall and temperature and to a sparse monitoring network.For both approaches, data availability is one of the key factors, and the integration of other and more continuous sources of data, such as weather data, would be beneficial to estimate the water budget in a reliable way.
Pragmatical aspects of the research can be summarized as follows: (1) different methods can quantify the single water budget terms; (2) uncertainty can be determined; (3) the detailed knowledge of the catchment framework, such as the hydrogeological setup and the anthropization degree (how much water is exploited and released), is crucial.The lack of this information does not allow the comparison between different methods.

Figure 1 .
Figure 1.Flow chart showing the research methodology (for symbols see Section 2.3 and 2.4).
coefficient for the latitude; T is the air temperature related to the ith month (in °C), and a 0.49239 1792 • 10 I 771 • 10 I 675 • 10 I is the exponent of Equation(11), which is based on the yearly heat index I ∑ .
Figure 5, as an example, fitted variogram models related to the Gaussian-transformed rainfall and temperature data for January are shown.

Figure 5 .
Figure 5. Variogram (variance vs. distance) examples.To the (left): rain data; to the (right): temperature data-both for January.The numbers on the variogram curves indicate the number of pairs.In Figure6, rainfall and temperature interpolations for January are shown; from left to right, Upper Limit, estimated values, and Lower Limit maps can be observed.Rainfall interpolations highlight the rainfall distribution typical of this climatic area: more intense

Figure 6 .
Figure 6.Geostatistical interpolations for 1986-2019 period in Foro a Molino Galasso sub-basin for January.(A) Rainfall, from left to right: Upper Limit, estimation, and Lower Limit.(B) Temperature, from left to right: Upper Limit, Estimation, and Lower Limit.

Figure 7 .
Figure 7. From top to bottom, example of monthly outflow, infiltration, and runoff in the Foro a Molino Galasso sub-basin.For each term (from left to right), the Upper Limit, estimation, and Lower Limit maps are shown.

Figure 8 .
Figure 8.Comparison between discharge values calculated with traditional and geostatistical methods and measured using hydrometers.(a) Ponte di Vacri sub-basin; (b) Molino Galasso sub-basin (see also Table2).

Table 1 .
Gauging station selected and their annual mean features.The asterisk indicates stations, where temperature has been obtained from regression lines.P: pluviometer; T: thermometer; H: hydrometer (see Figure3for the location).

Table 3 .
Yearly runoff values for the two sub-basins and for the whole basin calculated with traditional and geostatistical methods.The Turc, Turc modified, and Thornthwaite and Mather methods were applied for evapotranspiration (Thorn, Thornthwaite; LL, Lower Limit; Estim, Estimation; UL, Upper Limit).