## Abstract

Decision-making is a significant tool in water resources management applications. This work addresses the global management decision dilemma for the sustainability of the groundwater resources of a watershed: should stakeholders use groundwater for irrigation and human consumption or should they construct infrastructure, for example water reservoirs, for irrigation purposes? The former constitutes an easy but non-sustainable solution, while the latter protects the groundwater body from overpumping, avoids the associated overpumping penalties, and utilizes both surface and groundwater watershed resources. The main question arising in the second case relates to the amount of surface water that can be used taking into consideration water scarcity and potentially dry hydrological years. Therefore, this proposed decision-making framework will provide the best management solution for the water needs of an area based on the balanced use of surface and groundwater resources, considering the ecosystem sustainability and the surface and groundwater sustainability. In addition, this work can help decision-makers to examine and compare various scenarios using different approaches before making a decision regarding the cost and the capacity of a hydrologic/hydraulic project, and the varied economic charges that water table limit violations can cause inside an audit interval.

## Introduction

Groundwater sustainability relies on the optimal management of water resources. To address the decision dilemma, the required volume of water must be identified first and the available sources second. As mentioned above, the use of groundwater is the easy option, but the combined use of surface and groundwater ensures the sustainability of groundwater resources (Pulido-Velázquez *et al*., 2006; Wu *et al*., 2016). The main problem with surface water, in the absence of lakes, is the collection of water from rivers, streams, runoff and rainfall for which the development of water reservoirs is suggested. A second way to ensure groundwater sustainability is to impose penalties for overexploitation (Zekri, 2009; Madani & Dinar, 2012; Rinaudo *et al*., 2012; Mechlem, 2016; Molle *et al*., 2018).

A decision-making framework for sustainable water resources management is developed in this work. The framework addresses the optimal collection of the required water volume by taking into consideration the needs for irrigation (including industrial use) and human consumption, the availability of surface water in the watershed, the statistical characteristics of surface water availability with time as well as the state and availability of groundwater. The least-cost approach considers the action with the lowest economic and environmental cost, leading to the most advantageous decision. For the alternative of constructing a reservoir to collect surface water runoff, which reduces the impact of overexploitation of groundwater bodies, different reservoir volumes and the associated cost will be examined. In addition, overpumping penalties will be applied as a safety measure to protect groundwater from overexploitation. In the case where the use of groundwater pumping is more economical than the reservoir construction but results in substantial drop-down of groundwater levels, an alternative approach based on using both the collection of surface water and reduced groundwater pumping will be examined to ensure sustainability.

A combination of two approaches is proposed to ensure the sustainability of groundwater. The first is based on a risk analysis in terms of financial and environmental criteria and the second on a hydrological risk analysis. In the first approach, Bayesian decision analysis that will consider the prior behavior of the aquifer in terms of safe groundwater levels and the prior behavior of surface water will be applied for the decision-making process. In the second approach, statistical hydrology tools will be applied to evaluate the groundwater volume that must be pumped in the case where the reservoir fails to cover the water needs from the surface water.

Recently, a published work (Varouchakis *et al*., 2016) related these two topics to the sustainable groundwater management of a watershed. Considering the space–time aquifer behavior and the historical dry years of a watershed in Crete, Greece, Varouchakis (2016) developed a global decision-making tool under the principles of Bayesian risk analysis that aids the decision-making in groundwater management problems and ensures groundwater sustainability. The tool considered the construction of a reservoir and a set of penalty fines for overpumping violation.

The decision dilemma in this previous work considered only a scaled parabolic variation of the overpumping penalties vs the construction of a water reservoir and at a second stage the probability of dry years that will not provide sufficient water availability to the reservoir. The groundwater level variations of the aquifer as well as the rainfall variations at the basin and analysis of dry and wet years were studied in the work of Varouchakis (2016). In addition, in Varouchakis *et al*. (2016) the full mathematical framework of the Bayesian risk analysis was presented and applied. However, this study extends this previous work by considering all the potential parameters that affect the decision-making under uncertainty, providing an innovative approach. Variables such as groundwater cost, the lost value of groundwater (LGV) as a source, the reservoir's usable volume and the annual operational cost, and overall, the hydrological variables (rainfall and runoff) that would affect the functionality of a reservoir were considered. In addition, it considers the combined use of surface and groundwater to cover the water needs in terms of sustainability.

There is not a specific guide to set overpumping penalties. In some cases, in Mediterranean countries, a ban on groundwater use applies as an alternative to penalties that result in financial loss for stakeholders. On the other hand, the set of penalties is case specific and depends on the availability of the aquifer, the importance of groundwater shortage, and the frequency of overpumping in specific areas (Lerche & Paleologos, 2001; Madani & Dinar, 2012; Molle *et al*., 2018). In some cases, in Portugal and Oman, a very high fixed penalty is imposed, while in areas of Spain scalable penalties apply depending on the frequency of overpumping (Zekri, 2009; Rinaudo *et al*., 2012). In different areas of the USA such as Nebraska, Kansas, and Texas, low daily overpumping penalties are imposed which change under a parabolic structure as overpumping increases (Brozović & Young, 2014). Therefore, the penalty fees should vary and depend on the availability/vulnerability of the local groundwater resources, the frequency of overpumping, and the financial situation of the country or area. In this case study, a recent decision of the local water management authorities was applied setting fines for overpumping at individual wells of up to 300 euros (Local Department of Water Resources Management, 2018).

Regarding the groundwater cost, it is also case specific and depends on the local cost of allocation of rural water (please see the case study). On the other hand, for the LGV there is a common practice that involves the assessment of the LGV based on the value difference of water between rural and potable. However, usually the values are also subject to country or district rules and conditions (Paleologos, 2008; Fenichel *et al*., 2016). Overall, the calculation of the aforementioned costs involves subjectivity. However, this work proposes the general framework of a decision-making procedure that can be applied for different case studies considering the appropriate national and/or local financial terms and conditions.

Such an integrated approach for a decision-making framework applied in the water resources sector consists of new feedback. The proposed work falls within the goals of UNEP (United Nations Environment Programme), FAO (Food and Agriculture Organization), and other international organizations because it relates to groundwater sustainability which is one of their priorities.

## Nature, scope, and objectives

Groundwater is invaluable for the hydrological cycle and for ecosystem viability. The sustainability of groundwater is based on optimal management. To prevent the overuse of groundwater resources, local authorities around the world have developed exploitation schemes that apply penalties to the end users in the case of overpumping. In most of the cases, the penalties vary depending on the level and the frequency of overexploitation. Therefore, a decision-making dilemma for the stakeholders would be the infrastructure development at a watershed instead of overpumping penalties to ensure the sustainable use of groundwater.

Usually water management infrastructure involves reservoir construction. However, in the decision-making framework the historical hydrological characteristics of the watershed should be considered as well as reservoir construction cost. The presence of drought years, the statistical behavior of rainfall, runoff, and evapotranspiration should be studied in the design of a reservoir to substitute for groundwater use. All the aforementioned factors would be appropriately considered in the decision-making framework in terms of loss functions.

The target is for the decision-making framework to provide the least-cost decision between infrastructure development or the use of groundwater only, which may lead to aquifer overexploitation, to cover the required water needs. Then, if the second prevails, the hydrological characteristics of the basin will be used to develop a balanced use of surface and groundwater in terms of a reservoir that will cover only a part of the water needs and the rest from groundwater. In the case that a reservoir is available, then a hydrological risk assessment would be applied to ensure the necessary use of groundwater in terms of the statistical hydrological characteristics in case the reservoir fails to cover the needs.

This framework gives the potential to authorities to formulate penalties depending on the different applications and evaluate the volume of a reservoir that can cover, or partly cover, their water needs to ensure the sustainability of groundwater. If they want to be strict with the end-users, they may select high penalties for groundwater overexploitation making the reservoir construction the basic alternative. If not, they may select lower varying penalties but select a smaller volume reservoir that will partly cover the needs, but the cost will be attractive and advantageous compared to the penalties. That solution would require the use of groundwater to cover the needs, but in lower volumes it will avoid overexploitation and ensure sustainability.

The expected results should indicate that the hydrological probability uncertainty is the driving issue that determines the least-cost decision and that, depending on how the unknown probability is handled, the methodology may lead to a different least-cost decision. Thus, in contrast to practices that assess the effect of each proposed action separately considering only current knowledge of the examined issue, this framework aids decision-making by considering prior information and the sampling distribution of future successful audits.

## Decision-making related investigations

Most of the recent decision-making research on groundwater sustainability regards multi-criterion decision-making approaches or multi-objective optimization to integrate different objectives into the planning, management, and decision-making process. A variety of criteria in terms of economic, social, and environmental dimensions are applied for the analysis. Different management scenarios are proposed that include reductions in irrigated areas, optimization of pumping, improved irrigation efficiencies, increased system loss for groundwater irrigation, and changes in cropping patterns (Pathak & Hiratsuka, 2011; Geng & Wardlaw, 2013; Kumar *et al*., 2014; Rothman & Mays, 2014; Yeh, 2015). However, the difference in these studies compared to the one proposed is that they suggest a solution based on the required water demands and the status of groundwater resources, optimizing their use, but without considering the prior and posterior hydrological behavior of the watershed in terms of surface waters and groundwater.

Furthermore, hydrological probability risk assessment has been applied to determine under uncertainty an optimal management decision, in subsurface flow and transport (Tartakovsky, 2007), in surface waters and especially flood events (Salas & Obeysekera, 2013; Tartakovsky, 2013; Efstratiadis *et al*., 2015). However, this approach considers only the prior information to estimate statistical characteristics and the uncertainty of the variable of interest using parametric or non-parametric methods, i.e., probability runoff and rainfall to cover the water needs using a reservoir, but without considering the future uncertainty.

On the other hand, this work provides a decision-making framework under uncertainty (using probability density functions) that provides a context on how groundwater should be managed in each watershed based on the hydrological characteristics and the water needs. Therefore, a decision can be made on the groundwater volume that is required in excess of the available surface water considering its sustainable use. In advance, it includes the design of a reservoir to handle surface water according to the hydrological potential of the watershed in contrast to the use of groundwater that has an economic impact (i.e., overpumping penalties, pumping costs, lost value of groundwater) and environmental impact. Therefore, herein, Bayesian decision theory in association with the statistical hydrology risk assessment can provide the least-cost decision for groundwater resources management.

An initial application of Bayesian decision theory in hydrology was for the assessment of the costs of overdesign of a flood level in the face of flood frequency uncertainty (Davis *et al*., 1972). Since then, it has been used in many applications. For example, it has been used to determine optimal groundwater sampling frequencies (Grosser & Goodman, 1985) and in decision analyses to engineer design projects, groundwater flow and transport, and monitoring networks in which the hydrogeological environment plays an important role (Freeze *et al*., 1990). It has been used to address the problem of permitting waste sites under conditions of imperfect information (Marin *et al*., 1989; Medina *et al*., 1989) and for the engineering design of a groundwater interception well used to capture a contaminant plume (Wijedasa & Kemblowski, 1993). Moreover, it has been used for selecting the best experimental design for groundwater modeling and management design under parameter uncertainty (McPhee & Yeh, 2006) and for investigating the value of collecting hydraulic conductivity data for optimal groundwater resources management (Feyen & Gorelick, 2005).

In most of the early decision analysis studies, it was assumed that decisions would be made by a rational, financially driven decision-maker, who might be risk averse, but who would otherwise make decisions that maximized his or her economic position. However, decisions are strongly influenced by the profile of the decision-maker. Thus, water resources management experts need to be aware of the complexity of the decision process, the close relationship that exists between the technical input and the risk term in a decision analysis, and the widely differing views toward the methodology and value of risk calculations (Freeze, 2015).

This work will consider the probability risk that the rainfall and runoff in a watershed are not sufficient to cover the designed reservoir capacity and the related uncertainty; then, the water volume that is necessary to cover the watershed needs will come from the groundwater but without exceeding a sustainable aquifer level threshold. Such an integrated decision-making approach has not been met in hydrological applications and it consists of a useful framework for the sustainable management of groundwater resources.

## Statistical analysis

During the presence of drought years, the statistical behavior of rainfall, runoff, and evapotranspiration is very important in designing a reservoir to substitute the groundwater use. In this regard, the hydrological features of the studied watershed should be analyzed. In the absence of surface water resources, the collection of water from rivers, streams, runoff, and rainfall is needed for the development of water reservoirs. For this reason, a characteristic statistical analysis of annual rainfall data for a pilot basin is presented. Data covering a period of 70 hydrological years between 1945/1946 and 2015/2016 have been collected for 12 catchments in a studied basin area, Messara Valley, in Crete, Greece (Figure 1).

StatsDirect (V2.7.2, StatsDirect, Ltd, Altrincham, Cheshire, UK) software package operating on a PC platform (Casper Excalibur, Intel® Core^{TM} i7-7700HQ CPU, 2.81 GHz, 16 GB of RAM, 64-bit) was used for the statistical analysis of the annual rainfall data. In all calculations, spreadsheets of Microsoft® Excel® 2016 (Microsoft Inc., Redmond, WA) running under the Windows 10 system was used as an open database connectivity data source. An alpha (*α*) level of 0.05 (or 95% confidence) was used to determine the statistical significance in the analysis. Various descriptive statistics, such as minimum and maximum values, medians, means, standard deviations, variance coefficients, standard errors of means, skewness, and Kurtosis values for each subset were analyzed in order to characterize the selected basin groups and identify the extreme values. The annual rainfall values used as quantitative hydrological measures of each catchment in the basin area were statistically evaluated using appropriate statistical tests.

Prior to applying parametric (unpaired or two-sample and matched-pair *t* tests) or non-parametric tests (the Mann–Whitney *U* (or the Wilcoxon rank-sum) test or the Kruskal–Wallis test with the Conover–Inman method), the Shapiro–Wilk *W* and Levene's tests were consecutively conducted as preconditions to ensure that the annual rainfall subsets had a normal or non-normal distribution, and variances (or standard deviations) of the paired groups were homogeneous or unequal. In the case of the collected datasets not being normally distributed, a non-parametric test was applied. The Mann–Whitney *U* test (named after Henry Berthold Mann and Donald Ramson Whitney, a non-parametric test used to compare two population means that come from the same population and also used to test whether two population means are equal or not) or the Kruskal–Wallis test (named after William Kruskal and W. Allen Wallis) were used to compare the the considered subsets.

*U*test involves the calculation of the following test statistics (

*U*

_{x},

*U*

_{y}, and

*U*) for datasets

*x*

_{i}= {

*x*

_{1},

*x*

_{2},…,

*x*

_{n}} and

*y*

_{i}= {

*y*

_{1},

*y*

_{2},…,

*y*

_{n}}: where

*n*

_{x}is the size of dataset

*x*;

*n*

_{y}is the size of dataset

*y*; and

*R*

_{x}is the adjusted rank sum for dataset

*x*, and

*R*

_{y}is the adjusted rank sum of dataset

*y*.

If the observed value of *U* ≤ *U*_{critical} then the test is significant (at the *α* level), the null hypothesis (*H*_{0}: there is a no difference between the ranks of the two samples) is rejected in favor of the alternative hypothesis (*H*_{1}: there is a difference between the ranks of the two samples).

*n*

_{x}> 20 and

*n*

_{y}> 20), the

*U*statistic is approximately normally distributed:

*N*(

*μ*,

_{U}*σ*), where

_{U}*μ*and

_{U}*σ*are the mean and standard deviation of

_{U}*U*. In that case, the standardized value

*z*statistic (i.e., the

*z*

_{critical}values for a two-tailed test, for a significance level of

*α*= 0.01 are

*z*

_{critical}= ±2.58) can be calculated and interpreted (for the

*z*

_{critical}value >0) as follows:

Test results were assessed with various descriptive statistics such as two-tailed *p* values to reflect the statistical significance between paired groups. Based on the other descriptive statistics (i.e., minimum, lower quartile (Q_{1}), median (Q_{2}), upper quartile (Q_{3}), maximum) of independent samples, box-and-whisker plots were also plotted within the framework of StatsDirect software to assess the annual rainfall statistics in a pictorial manner.

## Bayesian decision-making methodology and discussion

Groundwater sustainability depends on the availability of surface waters and due to ecosystem viability only a part of them can be used. Thus, initially, a hydrological design should be performed considering the historical data of the watershed to determine the hydrological balance. The users of the proposed framework should determine the water needs of their area and then from the historical hydrological characteristics the potential water volume that can be collected in a reservoir to cover the needs. Any difference will be covered by groundwater. In addition, from the hydrogeological characteristics a groundwater level threshold can be set to establish a sustainable aquifer level budget.

Thus, first the hydrological statistics will be assessed to identify if the required water volume of the selected reservoir can be fulfilled.

The decision-making process involves two stages: state estimation equations that express the proposed actions and decision-making. For state estimation, first, all the state parameters are defined. However, in the Bayesian approach, a state parameter is an unknown quantity and is considered a random variable that must be determined. The procedure of estimating each parameter involves previous knowledge on the examined issue and the use of the subjective prior distribution that expresses the prior information for each state parameter. Next, the Bayesian risk function is obtained to estimate the least-cost decision or the decision with the minimum expected risk. The latter also applies in terms of a cost-benefit analysis procedure and denotes the preferable action. The Bayesian decision-making process follows these four steps.

**Step 1:** Set up the decision-making problem by introducing the possible actions set and the parametric space. Establish the expected loss function for each decision.

In this work, two actions are considered: (i) Action *Α*(0): use only groundwater and (ii) Action *Α*(1): reservoir construction with a volume *V*_{Ω} to cover water demands.

The use of groundwater only as a major source can easily lead to overexploitation. Overpumping violation policy is proposed to be based on a scaled linear function with scaling coefficient (*K*) that varies with the frequency (*n*) of violations because of the importance of the problem. *Y* is a random variable that expresses the total number of overpumping violations during an auditing period, more specifically the *Y* variable indicates when the water table of the aquifer is below that of a threshold considered to be the limit that distinguishes whether there is violation or not.

*GC*denotes the groundwater cost (pumping and volume) and

*LGV*the lost value of groundwater as a sustainable source as soon as it is removed from the aquifer (Paleologos, 2008).

*A(1)*the following applies: where

*C*is the construction cost of the reservoir of a volume

*V*

_{Ω}(more than one reservoir can be considered if necessary); and AC is the annual operational cost for the examined auditing period. If one considers the case that there is a risk (probability)

*θ*

_{1}for the reservoir not to provide sufficient water during drought years or the watershed may not have the appropriate surface water resources an additional cost

*M*is applied, denoting that an additional water supply (i.e., water transport, groundwater) should be considered. Thus,

*θ*

_{1}expresses the probability that the reservoir does not cover the necessary water demands.

**Step 2:** Provide the state of the goal function. If, at this step, the parameters are considered known, then the decision process is called a cost-benefit analysis, and Step 4 is directly applied. If not, then both Steps 3 and 4 are applied.

*A(0)*the goal function is expressed as follows (Varouchakis

*et al*., 2016): where , and

*N*the total audits.

Positive ENPLV leads to decision *A(0*), but negative to decision *A(1)*.

**Step 3:** Develop the subjective prior distributions for each parameter quantifying the previous information.

*Y*and

*Y*are the state parameters considered unknown, then Bayesian analysis is applied in terms of the Bayesian risk function that considers prior information in terms of probability density functions to determine

_{1}*Y*and

*Y*: where

_{1}*π(θ)*denotes the conjugative prior distribution in each case that depends on the fitted probability density function to the data. The probability that over-pumping or a drought year would occur, or the necessary surface water would not be available is denoted as a ‘success’ and as a ‘failure’ correspondingly.

Parameters *r* and *t* are extracted from the use of historical data of the aquifer's water level. More especially, we set a threshold *Δh* as the maximum water level drop due to pumping. If water level is reduced more than *Δh*, a violation event occurs. Thus, the parameter *t* denotes the total number of water level historical data, whereas the parameter *r* denotes the fraction of water level data that exceeds *Δh*.

If *R* is positive, then the decision *A(1)* is more risky than the decision *A(0)*, and thus, we need to redesign the reservoir. On the other hand, for negative values of *R*, we estimate the volume of groundwater needed to cover the water demands (see flowchart in Figure 2). However, the appropriate volume, *G*, must not exceed the groundwater threshold (*GW _{threshold}*). If

*Δh*is greater than

*GW*, then either a water supply demand rebalance is required or an additional water volume

_{threshold}*WT*needs to be supplied occasionally in the reservoir as an extra water source to cover the needs. Then, the decision-making process is re-examined to obtain the least-cost approach.

The Bayesian risk analysis methodology involves the application of a sampling model distribution and a prior distribution to estimate the least-cost decision. In Bayesian methods, when the sampling model such as in this work is the binomial distribution, then its conjugate prior is a beta distribution, while the posterior is also a beta form distribution (Berger, 1985; Lee & Sabavala, 1987; Lerche & Paleologos, 2001; Ferson, 2005). In a previous work (Varouchakis *et al*., 2016), the beta distribution was tested on the groundwater level variations of the case study aquifer and indeed fitted to the beta distribution. Application of the posterior distribution is not considered by the Bayesian risk analysis methodology (Lerche & Paleologos, 2001; Ferson, 2005). However, an extension of this work, after presenting a framework to define the least-cost approach for the management of groundwater resources, could apply the posterior beta form distribution to provide the posterior condition of the water bodies. In another case study applying Bayesian risk analysis, the sampling model distribution may differ (e.g., exponential, Poisson, etc.) as well as the relative conjugate (e.g., gamma), based on the aim of the work, for instance, threshold violations as a function of time.

In the decision-making process described in Figure 2, a pumping optimization is required in the area of study to retain the sustainability of the aquifer before the risk function of overpumping is applied, while the transfer of extra water (e.g., waste water treatment discharge) is considered in the risk function of the reservoir operation as the source of water that will cover the usable volume in case it cannot be covered by surface water availability. In addition, a potential transfer of water is applied to cover the needs in case the excess groundwater required affects the set threshold of overpumping. By following the procedure described above, the user manages to cover water demands in an optimum manner, by using the appropriate statistical and prior information.

A water source such as the waste water treatment discharge may already exist in a study area, but it is not always accessible to exploit because either the plant might be away from the cultivation areas and the reservoir or there is no available network to transfer the water to the area of interest. In fact, both problems exist in the area of this case study. Therefore, it is more useful to include this amount of water as the additional source that will cover the needs and to consider the cost to transfer it in the reservoir. Reservoirs in rural areas are usually constructed close to a river in order to collect the river flow. The variation of annual river flow has been considered in this work to design the usable volume of the reservoir. In this work, there is a significant distance between the waste water treatment plant and the reservoir. Thus, the cost of transferring the water to the reservoir was considered in *R(A(1)*.

## Reservoir's technical aspects and cost

*V*

_{Ω}. A significant factor to be estimated is the exploitation rate

*α*. The exploitation rate is given by the following equation: where

*V*is the total water demands for a period

_{D}*t*and

*V*is the total amount of runoff for the same period of time.

_{runoff}*V*

_{Ω}is the usable volume (m

^{3}), is the average annual amount of runoff water (m

^{3}),

*D*is the ratio (m

^{3}), CV is the variation coefficient of annual runoffs,

*Z*

_{p}is the normalized random variable following the normal distribution

*N(0,1)*for various probabilities

*p*, and

*p*is the probability that reservoirs can cover the water demands.

To estimate the surface reservoir usable volume, historical data for runoff are needed. In Table 1, data of 41 annual river discharges from the Messara watershed in Crete, Greece are presented as a working example. It is noted that the hydrogeological and hydrological settings of the basin area were fully described in a previous study (Varouchakis, 2016). The value of total annual water demands is estimated as *V _{Demand}* = 5 Mm

^{3}. The value of exploitation rate (

*α*) was estimated for a period of

*t*

*=*41 years. The statistical and design parameters regarding the discharge data are as follows:

*V*

_{min}= 0,

*V*

_{max}= 65.90 Mm³,

*V*

_{average}= 13.97 Mm³,

*α*

*=*36.69%, , and CV = 1.19.

Year | V (Mm³) | Year | V (Mm³) |
---|---|---|---|

1967–1968 | 32.000 | 1987–1988 | 20.442 |

1968–1969 | 13.070 | 1988–1989 | 11.741 |

1969–1970 | 11.840 | 1989–1990 | 4.511 |

1970–1971 | 11.840 | 1990–1991 | 5.211 |

1971–1972 | 9.199 | 1991–1992 | 2.618 |

1972–1973 | 7.646 | 1992–1993 | 0.000 |

1973–1974 | 4.315 | 1993–1994 | 7.067 |

1974–1975 | 4.797 | 1994–1995 | 4.926 |

1975–1976 | 23.612 | 1995–1996 | 11.995 |

1976–1977 | 6.016 | 1996–1997 | 2.387 |

1977–1978 | 60.622 | 1997–1998 | 0.880 |

1978–1979 | 21.005 | 1998–1999 | 1.332 |

1979–1980 | 26.330 | 1999–2000 | 0.025 |

1980–1981 | 65.901 | 2000–2001 | 0.995 |

1981–1982 | 46.457 | 2001–2002 | 1.520 |

1982–1983 | 19.165 | 2002–2003 | 12.025 |

1983–1984 | 29.619 | 2003–2004 | 3.520 |

1984–1985 | 48.742 | 2004–2005 | 0.140 |

1985–1996 | 13.930 | 2005–2006 | 0.120 |

1986–1997 | 10.951 | 2006–2007 | 0.270 |

Total | 558.780 |

Year | V (Mm³) | Year | V (Mm³) |
---|---|---|---|

1967–1968 | 32.000 | 1987–1988 | 20.442 |

1968–1969 | 13.070 | 1988–1989 | 11.741 |

1969–1970 | 11.840 | 1989–1990 | 4.511 |

1970–1971 | 11.840 | 1990–1991 | 5.211 |

1971–1972 | 9.199 | 1991–1992 | 2.618 |

1972–1973 | 7.646 | 1992–1993 | 0.000 |

1973–1974 | 4.315 | 1993–1994 | 7.067 |

1974–1975 | 4.797 | 1994–1995 | 4.926 |

1975–1976 | 23.612 | 1995–1996 | 11.995 |

1976–1977 | 6.016 | 1996–1997 | 2.387 |

1977–1978 | 60.622 | 1997–1998 | 0.880 |

1978–1979 | 21.005 | 1998–1999 | 1.332 |

1979–1980 | 26.330 | 1999–2000 | 0.025 |

1980–1981 | 65.901 | 2000–2001 | 0.995 |

1981–1982 | 46.457 | 2001–2002 | 1.520 |

1982–1983 | 19.165 | 2002–2003 | 12.025 |

1983–1984 | 29.619 | 2003–2004 | 3.520 |

1984–1985 | 48.742 | 2004–2005 | 0.140 |

1985–1996 | 13.930 | 2005–2006 | 0.120 |

1986–1997 | 10.951 | 2006–2007 | 0.270 |

Total | 558.780 |

*D*and

*Z*;

_{p}To cover water demands, the values of *Z _{p}* are given as 2.326, 1.645, 1.036, 0.842, 0.524, and 0.39, respectively, for the probabilities (

*p*) of 99, 95, 85, 80, 70, and 65%. In Figure 3, the estimated values of

*V*

_{Ω}are presented for various probabilities to cover water demands as a function of demands in irrigation water.

*V*

_{Ω}and their respective construction costs

*C*, which is representative of this case study as the required reservoir is for rural use and its capacity falls in the range of reservoirs examined. The relation resulting from the regression is given in the equation below: where

*C*is the reservoir cost in M€; and

*V*

_{Ω}is the usable volume in Mm

^{3}.

## Statistical assessment of annual rainfall data

A detailed set of descriptive statistics that are associated with the annual rainfall data for 12 catchments in the basin area are summarized in Table 2. The statistical results indicated that the annual rainfall values of all the collected samples ranged between 218.50 and 2,162 mm, with an average value of 715.80 ± 283.2 mm. The mean annual rainfall in the Messara watershed was found to be in line with those reported in other studies (Tsanis & Apostolaki, 2009; Varouchakis, 2016). In addition to the overall descriptive statistics, it is noted that a comparison of the annual rainfall values is an important task for the interpretation of the precipitation trends of different catchments in the region. Therefore, an extensive parametric/non-parametric-based statistical analysis was also conducted to assess the annual rainfall variations of the catchments. This analysis provides a better evaluation of the possible differences in precipitation profile of the basin system.

Basin^{a} | Mean (mm) | Standard deviation (mm) | Variance coefficient | Standard error of mean | Skewness | Kurtosis | Maximum (mm) | Median (mm) | Minimum (mm) |
---|---|---|---|---|---|---|---|---|---|

AV | 924.060 | 194.190 | 0.210 | 28.630 | 0.190 | 2.100 | 1,301.000 | 927.550 | 551.800 |

AK | 533.150 | 165.770 | 0.310 | 24.440 | 1.560 | 8.120 | 1,230.800 | 505.000 | 251.100 |

AS | 597.070 | 164.500 | 0.280 | 23.990 | 0.900 | 3.770 | 1,105.500 | 562.300 | 341.600 |

VA | 540.110 | 192.370 | 0.360 | 31.620 | 1.290 | 4.760 | 1,107.500 | 520.000 | 278.600 |

VO | 1,135.000 | 334.200 | 0.290 | 47.740 | 1.000 | 4.290 | 2,162.000 | 1,082.000 | 588.100 |

GE | 869.820 | 226.790 | 0.260 | 33.440 | 1.040 | 5.150 | 1,696.900 | 858.250 | 465.000 |

ZA | 831.290 | 217.520 | 0.260 | 29.070 | 0.970 | 4.640 | 1,602.500 | 832.150 | 416.500 |

KA | 693.860 | 220.470 | 0.320 | 35.300 | 1.790 | 7.680 | 1,442.000 | 653.500 | 341.800 |

LA | 573.950 | 169.340 | 0.300 | 25.530 | 1.560 | 6.560 | 1,235.100 | 536.050 | 344.800 |

MO | 691.290 | 192.570 | 0.280 | 29.030 | 0.830 | 4.310 | 1,339.900 | 676.400 | 363.000 |

PO | 504.200 | 145.310 | 0.290 | 18.450 | 1.280 | 7.170 | 1,135.400 | 486.100 | 218.500 |

MI | 444.010 | 131.700 | 0.300 | 49.780 | 0.490 | 2.020 | 660.600 | 418.500 | 305.800 |

Basin^{a} | Mean (mm) | Standard deviation (mm) | Variance coefficient | Standard error of mean | Skewness | Kurtosis | Maximum (mm) | Median (mm) | Minimum (mm) |
---|---|---|---|---|---|---|---|---|---|

AV | 924.060 | 194.190 | 0.210 | 28.630 | 0.190 | 2.100 | 1,301.000 | 927.550 | 551.800 |

AK | 533.150 | 165.770 | 0.310 | 24.440 | 1.560 | 8.120 | 1,230.800 | 505.000 | 251.100 |

AS | 597.070 | 164.500 | 0.280 | 23.990 | 0.900 | 3.770 | 1,105.500 | 562.300 | 341.600 |

VA | 540.110 | 192.370 | 0.360 | 31.620 | 1.290 | 4.760 | 1,107.500 | 520.000 | 278.600 |

VO | 1,135.000 | 334.200 | 0.290 | 47.740 | 1.000 | 4.290 | 2,162.000 | 1,082.000 | 588.100 |

GE | 869.820 | 226.790 | 0.260 | 33.440 | 1.040 | 5.150 | 1,696.900 | 858.250 | 465.000 |

ZA | 831.290 | 217.520 | 0.260 | 29.070 | 0.970 | 4.640 | 1,602.500 | 832.150 | 416.500 |

KA | 693.860 | 220.470 | 0.320 | 35.300 | 1.790 | 7.680 | 1,442.000 | 653.500 | 341.800 |

LA | 573.950 | 169.340 | 0.300 | 25.530 | 1.560 | 6.560 | 1,235.100 | 536.050 | 344.800 |

MO | 691.290 | 192.570 | 0.280 | 29.030 | 0.830 | 4.310 | 1,339.900 | 676.400 | 363.000 |

PO | 504.200 | 145.310 | 0.290 | 18.450 | 1.280 | 7.170 | 1,135.400 | 486.100 | 218.500 |

MI | 444.010 | 131.700 | 0.300 | 49.780 | 0.490 | 2.020 | 660.600 | 418.500 | 305.800 |

^{a}AV, Agia Varvara; AK, Agios Kyrillos; AS, Asimi; VA, Vagionia; VO, Vorizia; GE, Gergeri; ZA, Zaros; KA, Kapetaniana; LA, Lagolio; MO, Moroni; PO, Pompia; MI, Mires.

As seen in Table 3, a total of 66 subsets were statistically compared to determine differences in the annual rainfall variations. Statistical analysis of the precipitation data showed that there was insufficient evidence for 14 (21% of total) subsets (Agia Varvara and Gergeri: *p**=* 0.1421; Agios Kyrillos and Asimi: *p**=* 0.0663; Agios Kyrillos and Vagionia: *p**=* 0.8331; Agios Kyrillos and Lagolio: *p**=* 0.3091; Agios Kyrillos and Pompia: *p**=* 0.3981; Agios Kyrillos and Mires: *p**=* 0.1633; Asimi and Vagionia: *p**=* 0.0506; Asimi and Lagolio: *p**=* 0.4021; Vagionia and Lagolio: *p**=* 0.2090; Vagionia and Pompia: *p**=* 0.5435; Vagionia and Mires: *p**=* 0.1928; Gergeri and Zaros: *p**=* 0.3004; Kapetaniana and Moroni: *p**=* 0.8125; and Pompia and Mires: *p**=* 0.2509) to suggest a significant difference in the annual rainfall values (*p* > 0.05). On the other hand, results from parametric/non-parametric analyses revealed that there were significant differences between 52 (79% of total) subsets (*p* < 0.05), and the null hypothesis was rejected in favor of the alternative hypothesis (*α* = 0.05) for these subsets.

Pairs of groups (basins)^{a} | Statistical outputs^{b,c} for annual rainfall subsets | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

AV | AK | AS | VA | VO | GE | ZA | KA | LA | MO | PO | MI | |

AV | - | p < 0.001 | p < 0.001 | p < 0.001 | p = 0.001 | p = 0.142 | p = 0.014 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 |

AK | - | p = 0.066 | p = 0.833 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p = 0.309 | p < 0.001 | p = 0.398 | p = 0.163 | |

AS | - | p = 0.051 | p < 0.001 | p < 0.001 | p < 0.001 | p = 0.010 | p = 0.402 | p = 0.011 | p = 0.003 | p = 0.018 | ||

VA | - | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p = 0.209 | p < 0.001 | p = 0.544 | p = 0.193 | |||

VO | - | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | ||||

GE | - | p = 0.300 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | |||||

ZA | - | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | ||||||

KA | - | p = 0.002 | p = 0.813 | p < 0.001 | p = 0.001 | |||||||

LA | - | p = 0.002 | p = 0.048 | p = 0.046 | ||||||||

MO | - | p < 0.001 | p = 0.001 | |||||||||

PO | - | p = 0.251 | ||||||||||

MI | - |

Pairs of groups (basins)^{a} | Statistical outputs^{b,c} for annual rainfall subsets | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

AV | AK | AS | VA | VO | GE | ZA | KA | LA | MO | PO | MI | |

AV | - | p < 0.001 | p < 0.001 | p < 0.001 | p = 0.001 | p = 0.142 | p = 0.014 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 |

AK | - | p = 0.066 | p = 0.833 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p = 0.309 | p < 0.001 | p = 0.398 | p = 0.163 | |

AS | - | p = 0.051 | p < 0.001 | p < 0.001 | p < 0.001 | p = 0.010 | p = 0.402 | p = 0.011 | p = 0.003 | p = 0.018 | ||

VA | - | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p = 0.209 | p < 0.001 | p = 0.544 | p = 0.193 | |||

VO | - | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | ||||

GE | - | p = 0.300 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | |||||

ZA | - | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | p < 0.001 | ||||||

KA | - | p = 0.002 | p = 0.813 | p < 0.001 | p = 0.001 | |||||||

LA | - | p = 0.002 | p = 0.048 | p = 0.046 | ||||||||

MO | - | p < 0.001 | p = 0.001 | |||||||||

PO | - | p = 0.251 | ||||||||||

MI | - |

^{a}AV = Agia Varvara, AK = Agios Kyrillos, AS = Asimi, VA = Vagionia, VO = Vorizia, GE = Gergeri, ZA = Zaros, KA = Kapetaniana, LA = Lagolio, MO = Moroni, PO = Pompia, and MI = Mires.

^{b}Prior to applying parametric test, the Shapiro–Wilk *W* and the Levene's tests were consecutively conducted, as preconditions, to evaluate whether the experimental subsets had either a normal or a non-normal distribution, and to evaluate whether the variances (or standard deviations) of the paired groups were either homogeneous or unequal.

^{c}*p* values > 0.05 were considered to be insignificant (indicated as bold).

Finally, all local differences between the independent groups were analyzed graphically using box-and-whisker-plots (Figure 5). As seen in Figure 5, according to the minimum, lower quartile (Q1), median (Q2), upper quartile (Q3) and the maximum categories, the box-and-whisker plots suggest that the distributions of independent samples are quantitatively different for most of the paired groups, except the above-given 14 subsets. To conclude, the statistical analysis of the annual rainfall data revealed within the 95% confidence that about 79% of the total mutual comparisons (*n* = 66) were found to be statistically noticeable (*p* < 0.05). This corroborated the non-uniform hydrological structure of the Messara basin, which is mainly ascribed to its topographic heterogeneity. A detailed hydrogeological description of the Messara Valley can be found in previous studies (Vardavas *et al*., 1996; Peterek & Schwarze, 2004; Kritsotakis & Tsanis, 2009; Kritsotakis, 2010; Varouchakis, 2016). Therefore, the mean annual rainfall variation would be obtained from independent time series data providing a statistically significant and reliable rate of dry years, depending on a case-specific threshold for the region, that would correspond to the parameter *θ*_{1} of the action *A(1)*.

## A case study application of the proposed decision-making process

Considering the methodological steps previously described for the proposed decision-making process, a realistic application was performed in the study area of the Messara watershed. The water demand in the area is *V _{Demand}* = 5 Mm

^{3}annually. The hydrological characteristics of the basin in terms of the frequency of dry years (Varouchakis

*et al*., 2016) denote a 70% probability (

*Zp*) for the demand to be covered from the reservoir. From Equation (22) and Figure 3, the annual usable volume

*V*

_{Ω}is then calculated equal to 2 Mm

^{3}. Thus, the construction cost

*C*of the reservoir from Equation (23) is calculated equal to 6 million euro, while the annual operational

*AC*was calculated equal to 1% of the construction cost (Stephenson, 2012). The probability

*θ*for the reservoir not to provide sufficient water during dry years is close to 30%, calculated from the number of dry years in the area during the last 41 years (Varouchakis

_{1}*et al*., 2016, 2018), while the cost of transferring water from another source to cover the demands, e.g., local waste water treatment plant, was calculated equal to 88,000 euro (Varouchakis & Palogos, 2017). The treated waste water is distributed free of charge, but the transferring cost with pipes needs to be covered.

On the other hand, in the area, around 1,000 wells are in operation (Kritsotakis & Tsanis, 2009). Recently, the local authorities announced a potential penalty of 300 euro for overpumping to every well operator (Local Department of Water Resources Management, 2018). According to available information, when overpumping of the aquifer occurs in the area, more than 70% of the wells (i.e., 700) violate the agreed rules. This work proposes a scaled variation of the penalty to provide a period of compliance to the users. A monitoring period of 24 months is suggested to examine the stakeholders' behavior on a monthly basis. Therefore, in Equation (7) *K1* = (100 × 700) 70,000 euro for *n* ≤ 12, *K2* = (200 × 700) 140,000 euro for 12 < *n* ≤ 18, *K3* = (300 × 700) 210,000 euro for 18 < *n* ≤ 24. The probability *θ _{0}* of overpumping follows the binomial distribution (monthly monitoring of overpumping, i.e., success or failure) with a beta prior distribution as conjugate according to the literature (Berger, 1985), identified also from the optimal fit of groundwater level variations of the last 35 years in the area to the beta distribution. In addition, the groundwater cost

*GC*to cover the demand is 0.08 euro/m

^{3}(i.e., 0.08 × 6,000,000 m

^{3}) for the area of the case study, while the lost value of groundwater as a sustainable source

*LGV*can be calculated considering its potential value as fresh potable water (Paleologos, 2008). The potable water in the area costs 0.38 euro/m

^{3}. Thus, the

*LGV*cost is equal to (0.38–0.08) × 6,000,000 m

^{3}.

Considering the available information and applying the proposed methodology, it is obtained that for up to 18 overpumping violations, action *A(0)* applying scaled penalties is more affordable compared to *A(1)* reservoir construction. For more violations, the cost of reservoir construction is lower than the overpumping penalty cost. According to the decision-making flowchart (Figure 2), an assessment follows for the impact of the groundwater level decline in the aquifer considering the withdrawn amount to cover the demand. The study area *A* is 50,000,000 m^{2} and the storativity coefficient *s* equal to 0.05 (Varouchakis, 2016). Therefore, the expected aquifer level decline was calculated, equal to 1.2 m/yr, less than 5 m/yr that may affect the set aquifer level threshold which is 25 meters above sea level (Varouchakis *et al*., 2016). Therefore, considering the sustainable water resources policy that the local authorities desire to follow and the history of aquifer overpumping in the area, investment in reservoir construction is suggested.

However, this work involves some limitations and uncertainties in terms of hydrological and economical parameters' calculation. Terms such as the pumping cost, the lost value of groundwater, the reservoir construction cost, and the annual reservoir operational cost are case specific. The latter three though are the most uncertain economic parameters as their calculation depends on assumptions that are affected from the quality of groundwater, technical, and operational costs that vary with case and time of construction and operation. On the other hand, a critical and uncertain hydrological parameter is the probability that the reservoir's usable volume is not covered because of, e.g., the presence of a dry year or reduced river flows. The presence of environmental flows is also a significant parameter in the proposed method as it is applied for the calculation of the reservoir's usable capacity. Also, overpumping penalties' cost is an uncertain parameter. Usually, they are not a function of national rules. In specific cases where the local water resources management authorities set penalties for individual overpumping, they depend on the frequency of overpumping violations and the significance of the groundwater availability problem (see Introduction section). The beta prior distribution was employed because it is the conjugate of the binomial distribution (Berger, 1985) that was applied as the sampling (monitoring) model to examine the presence or not of overpumping. However, the monitoring model can vary and, instead, to examine a success or failure situation, it can examine when a success situation will occur or for how long. For that assumption, then, a different likelihood and prior applies such as the Poisson distribution, while the conjugate prior is now the gamma distribution. Other parameters such as the normalized *Zp* variable considered in the calculation of *V*_{Ω} follow explicitly the normal distribution for various probabilities that a reservoir cover a specified water demand. Overall, this work provides a flexible and integrated decision framework that can consider a variety of parameters that affect the combined use of surface waters and groundwater under uncertainty to cover potential water needs in terms of sustainability.

## Conclusions

The output of this research work is a global decision-making framework for sustainable and economic viable groundwater management. The methodology takes into consideration the need to cover the demands by only pumping water from the aquifer with the risk of overexploiting the resources or to cover the water demands by constructing a surface reservoir. A detailed approach to designing the usable reservoir volume and to determine the associated cost is presented. Each decision is expressed by an expected loss function that is described in monetary terms. The expected losses and, consequently, the expected risks of any decision can be estimated providing the user with the potential to assess different scenarios depending on the application. The proposed framework combines historical hydrological data with statistical approaches to quantify any useful available information to offer a holistic methodology for water resources management. Therefore, decision-makers can apply or modify appropriately the proposed framework to perform a risk analysis of the potential considered actions depending on the case study. Terms such as the lost value of groundwater, the pumping cost, and the annual reservoir operational cost were appropriately considered. The authors applied the proposed framework for the described case study to provide a guide to the authorities or the stakeholders to employ it for the consideration of different management scenarios.

## References

*(amended in 2016, in Greek)*

*(in Greek)*