Risk analysis framework for the optimum remediation of a contaminated aquifer under uncertainty: application in Lake Karla aquifer, Thessaly, Greece

A risk analysis framework is proposed for the optimum remediation of a contaminated aquifer under hydrogeological uncertainty. The limited information and the spatial variation of hydraulic conductivity in a real-world large-scale aquifer create uncertain conditions for decision-making when remediation schemes ought to be accompanied by the minimum possibility of failure. The primary concern is focused on safeguarding public health when groundwater is used for urban drinking purposes from a contaminated aquifer. The proposed framework is based on the conjunctive use of stochastic simulation–optimization modelling followed up by a risk analysis application on remediation trade-offs. The framework includes three main steps/procedures: (i) the model formulation of multiple realizations of groundwater flow and contaminant transport, (ii) the optimal positioning and operation of the clean-up wells determined by the method of stochastic optimization, and (iii) the risk analysis of the optimum remediation strategies through a proposed decision model, so as the one with the minimum cost and risk of failure is chosen as the most appropriate. The proposed framework is tested for two scenarios of nitrogen fertilizer application in the cultivated areas. The strategic target is the groundwater nitrate concentration minimization in an area where exceedances of nitrate concentrations have been observed and water supply wells have been operating for the last twenty years satisfying domestic needs. The results demonstrate that, when decision-making is under hydrogeological uncertainty, the combined use of stochastic optimization and risk-based decision analysis can commend the remediation strategy with the minimum cost and the highest possibility of success.


Introduction
Groundwater resources of rural basins in the Mediterranean region are under severe pressure, with nitrate pollution being one of their most serious threats (Cinnirella et al. 2005;Rebolledo et al. 2016;Ducci et al. 2019;Cil et al. 2020;Lyra et al. 2021;Orellana-Macías et al. 2021). High nitrate concentrations are found in regions with significant agricultural production, where nitrogen fertilizers are often excessively used to maximize yield, and in depository areas of animal fecal material (MEE-GSNEW 2020). The major reason of groundwater being contaminated by nitrates seems to be the long-term and excess fertilization practices to achieve crop-yield magnification (Hass et al. 2017;Marsala et al. 2021;Ncibi et al. 2021;Lyra et al. 2022). Nitrate contamination in aquifers is primarily caused by a progressive deposition of nitrates, which can reach levels that are prohibitive for urban water supply (World Health Organization International Agency for Research on Cancer 2020). Nitrogen can pile up in groundwater, especially in confined aquifer systems and aquifers with extremely low recharge rates. Most often, it is expensive for local authorities to invest in additional nitrate clean-up treatment plans of contaminated water. The health risk is especially significant in rural locations where contaminated groundwater is pumped from groundwater wells for drinking water supply since approximately 5% of human ingested nitrates are endogenously transformed to N-nitroso compounds, that have been linked to a variety of diseases such as methemoglobinemia in new-borns (blue baby syndrome), types of cancer and thyroid problems (World Health Organization 2017; Ward et al. 2018; World Health Organization International Agency for Research on Cancer 2020). The World Health Organization (WHO) has set the short-term limit to 50 mg/L for drinking water supply based on studies to prevent infant methemoglobinemia, and this limit has also been adopted by the European Union (World Health Organization 2017; World Health Organization International Agency for Research on Cancer 2020). The Nitrates Directive 91/676/EC is one of the European Union's (EU) initial legislative actions aimed at decreasing pollution and enhancing water quality along with the Water Framework Directive (2000) and the subsidiary Groundwater Directive (2006). However, achieving ''good'' quantity and quality status in EU waters is greatly constrained by diffuse sources of pollution (European Commission 2013). Several EU countries have been convicted of failing to fulfill their obligations to the European Commission's Nitrates Directive (91/676/EEC). In each case, the European Court of Justice has ordered member states to strengthen their regulations to comply with the Nitrates Directive. In recent years, France, Greece, Poland, and Luxembourg have been taken to court over nitrate pollution and Estonia has been warned (European Commission 2015). Nitrates Directive seeks to safeguard water quality by avoiding nitrate contamination of groundwater and surface water from agricultural sources by encouraging the implementation of Best Management Practices (BMPs) in agriculture. The nitrates concentration level of 50 mg/L has been declared as the maximum limit by the Greek and international legislation. Even at lower concentrations (more than 25 mg/L), concern is aroused for the long-term use of water for drinking purposes (World Health Organization International Agency for Research on Cancer 2020; National Centre for the Environment and Sustainable Development 2018).
In cases when contaminated groundwater is used for urban drinking purposes, public health issues are raised (Bilgehan and Berktay 2006;Gutiérrez et al. 2018;Postigo et al. 2021;Sidiropoulos et al. 2021) making protection and remediation obligatory for public health safeguarding. Contaminants in water not only limit its utility value but also place an added economic burden on society through both primary treatment costs and secondary impacts on the economy (Alexander and Ndambuki 2021).
Management of groundwater pollution problems is often aimed at limiting the migration of plume off-sites, isolating, and containing a source area from further leaking, or treating polluted aquifer to bring it up to acceptable standards (Matott et al. 2006;Mategaonkar and Eldho 2012;Samadi-Darafshani et al. 2021). Pump and Treat techniques have been widely used worldwide for the remediation of polluted aquifers (Kazemzadeh-Parsi et al. 2015;Sreekanth et al. 2016;Yudao et al. 2017;Antelmi et al. 2020;Sidiropoulos et al. 2020). Coupled simulation-optimization (S-O) methods have been successfully applied to identify optimal groundwater management solutions for various pollution problems (Sreekanth and Datta 2010;Mategaonkar and Eldho 2012;Sreekanth et al. 2016;Norouzi Khatiri et al. 2020). Many types of uncertainties occur in the simulation and management of real-world, large-scale aquifer systems. The most basic source of uncertainty in groundwater flow and contaminant transport simulation originates from hydraulic conductivity heterogeneity (Freeze 1975;Wagner and Gorelick 1989;Dagan 1990;Gelhar 1993;Tiedeman and Gorelick 1993). From an engineering standpoint, protecting and remediating contaminated aquifers under uncertainty is the challenge of choosing the optimum solution from several alternatives that meet all environmental, economic, technical, political, and legal commitments. This indicates that the solution must satisfy the constraints and the objective of the proposed management problem at optimum. A contemporary and widely used methodology for such problems is the combination of stochastic groundwater flow and transport simulation models with optimization techniques targeting the selection of the optimal engineering scenario from a series of alternatives (Gorelick et al. 1984;Ahfeld et al. 1988;Mylopoulos et al. 1991;Park et al. 2014;Rona et al. 2018). But stochastic optimization uses only one technological strategy at a time (Freeze and Gorelick 1999). Decision analysis is another methodology that has been successfully applied to groundwater remediation problems under uncertainty (Massman and Freeze 1987a, b;Freeze et al. 1990aFreeze et al. , b, 1992Freeze and Gorelick 1999;Yenigül et al. 2006;Teng et al. 2019) and is more suitable for riskbased remediation design than stochastic optimization, as it allows all benefits, costs, and impacts of uncertainty to be included in a single objective function.
The study of Freeze and Gorelick (1999) was the first that compared the stochastic optimization and decision analysis frameworks. The authors indicated that their fundamental difference is the availability of the technological strategies that each framework considers. The problem was the remediation of contaminated aquifers under hydraulic conductivity uncertainty with the use of pump and treat systems. The target focused on seeking the least cost and low-risk remedial strategy. The authors indicated that the fundamental difference is the availability of the technological strategies that each framework considers. The decision analysis uses a variety of technological strategies and one of them is selected as the best. Stochastic optimization uses only one technological strategy at a time. The study indicated the problem cases for which the two approaches are suitable. But the remarkable conclusion was that they finally suggested the combined use of stochastic optimization and decision analysis for remediation problems. Mylopoulos et al. (1999) evaluated and compared these two approaches to the same groundwater quality management problem of a contaminated aquifer in Greece under hydraulic conductivity uncertainty. Stochastic optimization was evaluated as weaker than decision analysis due to two reasons: (1) the risk of failure due to uncertainties cannot be considered as probabilistic costs, but only as constraints of the management problem and (2) the optimization leads to economically non-profitable solutions, given that the beneficial effects of the dispersion factor to reduce the concentration of contaminants below the upper limit cannot be considered. Decision analysis' unique weakness is the selection of the best solution among a small, predetermined set of alternatives. For the above-mentioned reasons, Mylopoulos et al. (1999) proposed a joint stochastic optimization and decision analysis framework in a way that one methodology covers the weakness of the other. Although the flow rate of the clean-up wells was determined by the optimization problem, their locations were predefined by the researchers. While in all the mentioned studies (except that of Mylopoulos et al. 1999) stochastic optimization and decision analysis are treated as two separate techniques for risk-based groundwater quality management problems, this paper proposes their coupling in a successive row of action. The aim of this study is the remediation of a nitrate-contaminated aquifer under hydraulic conductivity uncertainty through an optimum and with the least risk of failure strategy. For that reason, a risk-based framework is proposed and applied in the Lake Karla aquifer in Eastern Thessaly, Greece. The selected methodology for nitrates removal is the stabilization of the plume with the use of pumping wells pairs aimed at controlling the hydraulic gradient along the contaminated plume boundary. The novelty of this work is the coupling of stochastic optimization and decision analysis under a two-stage procedure. In the first stage, the stochastic optimization determines the optimal positioning and operation with the minimum cost of the clean-up wells for every equally probable aquifer's realization according to the insertion of an objective function and of the constraints. In that way, stochastic optimization provides risk analysis a predetermined set of alternatives. In the second stage, risk analysis introduces environmental, economic, technical, political, and legal commitments to every optimum alternative realization, according to a proposed decision model. These commitments are valued in economic terms. The engineering remediation strategy characterized by the minimum cost and the least probability of failure is proposed for the successful removal of nitrates.

Models and methods
The proposed risk-based framework is a modelling system formed by the linking of geostatistical, simulation, and management models based on the combination of contemporary and widely used methodologies such as integrated modelling, conjunctive use of stochastic simulationoptimization modelling, and risk analysis. The steps followed in the procedure are: • Integrated modelling: the application of simulation models at the watershed level and their linking through data exchange evaluate the study area's water potential. • Stochastic simulation: multiple realizations of groundwater flow and contaminant transport are created from equally probable maps of spatially variable hydraulic conductivity generated by geostatistical tools. • Stochastic optimization: the optimal positioning and operation of the clean-up wells are determined for each aquifer's realization according to an objective function and a set of constraints.
• Risk analysis: the optimum alternative strategies generated from the previous procedure are tested by risk analysis through a proposed decision model so as the one with the less risk of failure and the minimum cost is chosen the most appropriate.

General modelling framework
The flow diagram of the proposed risk-based framework is presented in Fig. 1. The models' coupling is achieved by data linking. The components are the UTHBAL surface hydrology model (Loukas et al. 2007), the UTHRL reservoir operation model (Loukas et al. 2007), the LAK3 lakeaquifer interaction routine (Merritt and Konikow 2000), the GEPIC agronomic model (Liu et al. 2007), the Geostatistical Library (Deutsch and Journel 1998), the MODFLOW groundwater hydrology code (Harbaugh et al. 2000), the MT3DMS contaminant transport code (Zheng and Wang 1999) and the groundwater management model called GWM (Ahlfeld et al. 2005). The procedures of integrated modelling and stochastic simulation are implemented by the linking of simulation and geostatistical models. Inputs of the UTHBAL model are the mean monthly areal precipitation, temperature, and evapotranspiration. The model calculates the monthly surface runoff and the natural groundwater recharge. The average irrigation returned flow added to the recharge rate. The UTHRL model estimates the stored volume, the overflow, and the area of the lake/reservoir using as inputs the inflows from nearby water resources, the basin's surface runoff of the UTHBAL model, and the precipitation to the reservoir. The LAK3 model simulates the leakage of the lake/reservoir to groundwater using as inputs the Fig. 1 Flow diagram of the proposed risk analysis framework indicating the models, the data exchange, and the methodological procedures estimates of the stored volume in the lake/reservoir. The Geostatistical Library produces equally probable realizations of the spatially distributed hydraulic conductivity and is presented in the next section. The MODFLOW code receives the groundwater recharge, the irrigation returned flow, and the seepage from the lake/reservoir and calculates the groundwater flow for each equally probable realization of the hydraulic conductivity in a stochastic mode. The GEPIC model uses as inputs the groundwater recharge and the irrigation return flow, and according to the nitrogen fertilizer application rate of each crop estimates the nitrates leached into the groundwater. The MT3DMS code uses the nitrate leaching rates from GEPIC and the multiple groundwater flow realizations from MODFLOW and determines additional realizations of the extent and the spatial distribution of the nitrate plumes in a stochastic mode. Since the study focuses on the nitrates removal from groundwater under hydrogeologic uncertainty, the core of the modelling system is the MODFLOW and MT3DMS codes.
The stochastic optimization is achieved by coupling the two groundwater codes (MODFLOW, MT3DMS) with the groundwater management model (GWM). This model has been developed to be compatible with MODFLOW and uses a response-matrix approach to solve several types of linear, nonlinear, and mixed-binary linear groundwater management problems. An optimization exercise is formulated by a set of decision variables, an objective function, and a set of constraints targeting the optimal positioning and operation with the minimum cost of the clean-up wells. This exercise is executed by GWM for each multiple equally probable aquifer realization generated by stochastic simulation. So, GWM receives the multiple realizations of groundwater flow to produce an equal number of optimum solutions. These solutions are imported into the two groundwater models to generate the optimum alternatives for aquifer quality and quantity status. A reliability analysis is achieved to evaluate the clean-up strategies.
The risk analysis procedure is the final step of the proposed risk-based framework. It uses a predefined set of remediation strategies from the stochastic optimization procedure. Environmental, economic, technical, political, and legal commitments, valued in economic terms, are introduced to every optimum alternative according to a proposed decision model. The engineering remediation strategy with the minimum cost and the least possibility of failure is proposed for the successful removal of nitrates. A more detailed description of the aforementioned models can be found in Sidiropoulos et al. (2013), Mylopoulos and Sidiropoulos (2016), and Sidiropoulos et al. (2019).

Stochastic simulation
The geological heterogeneity and the limited scatter data of hydraulic conductivity generate uncertainty in the simulation of groundwater flow and the nitrates' concentration and transport. In this case, the uncertainty analysis is treated in a geostatistical fashion with the use of spatially distributed, autocorrelated, and conditioned conductivities (Freeze and Gorelick 1999). The geostatistical approach has been widely and successfully used for hydraulic conductivity uncertainty assessment. In the proposed framework the Geostatistical Library is chosen and applied through a two-step process: (i) the values of hydraulic conductivity at unsampled locations are estimated with Simple Kriging, and (ii) equally probable realizations of the spatially distributed hydraulic conductivity are produced by the conditional Sequential Gaussian SIMulation (SGSIM) algorithm (Journel and Alabert 1989;Deutsch and Journel 1998).
She Conditional Probability Distribution Function (CPDF) of a random variable is determined by the scatter data, their histogram, and the model variogram. Since hydraulic conductivity is the random variable, the secondorder stationarity hypothesis is assumed, and its spatial distribution is approached by a log-normal function characterized by the mean value and the covariance (Dagan 1982;Hoeksema and Kitanidis 1985). For this reason, the values of hydraulic conductivity are transformed into a lognormal distribution. The application of the conditional SGSIM algorithm is achieved with the use of a spaced grid that covers the whole study area, and a random pattern that runs through it. Simple Kriging uses the log-normal transformed values of hydraulic conductivity to assess the expected value of the probability distribution function for the grid nodes' lacking sample. The shape of the CPDF is known based on the second-order stationarity assumption, allowing the estimation of any possible value probability. Thus, the conditioning is performed by using both scattered data and every simulated value at any unsampled location. This introduces spatial continuity to the simulation process, as needed by the spatial covariance function. The equally probable hydraulic conductivity realizations have the same statistical characteristics.
Graphical and statistical features of the scatter data are performed. A series of experimental semi-variograms is calculated in different directions so that the hydraulic conductivity's anisotropy direction if exists, to be identified. GSLIB generates a variogram map for each experimental semi-variogram. The general anisotropy trend is revealed by the direction in which the variogram map elongates. Finally, several model semi-variograms are tested and the one that best fits the experimental semi-variogram and gives the minimal variance of the estimation error is selected.
MODFLOW and MT3DMS models use the equally probable realizations of the spatially distributed hydraulic conductivity generated from the Geostatistical Library and produce additional realizations of groundwater flow and the extent and the spatial distribution of the nitrate plumes.

Stochastic optimization
A stochastic optimization procedure is chosen to be executed before risk analysis due to its great advantage, which is the determination of a globally optimal solution, and due to other reasons, that will be presented below in the manuscript. The stochastic optimization is applied to the multiple realizations of groundwater flow and nitrate concentration and dispersion as generated by groundwater simulation codes (MODFLOW, MT3DMS). This is achieved by the combination of groundwater models and groundwater management model (GWM) applying the response matrix method. The selected methodology for nitrates removal regards the stabilization of the plume with the use of pumping wells pairs that aim at the control of the hydraulic gradient along the boundary of the contaminated plume. Monte Carlo approach is implemented to solve the optimization exercise for each realization as proposed by Wagner and Gorelick (1989). In that way, an equal number of optimum solutions result. The target of the optimization exercise is the minimization of pumping wells' operational cost while ensuring the hydraulic gradient inversion toward the interior of the plume and the prevention of high hydraulic head drawdown. A general form of the optimization equations can be presented as: Realization 2 : Realization n : is the total pumping rate vector of the clean-up wells, A, is the operational cost of clean-up wells per unit pumping, . .:m is the vector of a lower limit of hydraulic heads difference between the M checkpoint pairs outside and inside the plume ensuring a specific limit value of the local hydraulic gradient, which always has a direction towards the inside of the plume. This lower limit is defined as the product of the minimum required decline at the specific control location and the distance between the two points. hðK i Þ; i ¼ 1. . .n, is the vector of the simulated hydraulic heads, h i ; i ¼ 1. . .n, is the vector of the hydraulic heads' constraint.

Risk analysis
The main advantage of risk analysis is the introduction of all benefits, costs, and effects of uncertainty to the objective function. This is achieved via the coupling of the stochastic application of groundwater models with a decision model. To overcome the decision analysis' unique weakness, the decision model is applied to the optimum alternatives generated by stochastic optimization.
Risk analysis quantifies the probability of failure for each of the optimal decontamination scenarios of the aquifer. At this point, the proposal made to define the decision model, is to separate the failure risk into two components. The first component concerns the environmental risk that arises exclusively from the presence of the contaminant plume in the aquifer. The second component concerns the most realistic risk which is the inability to supply the community (village or town) with drinking water for a short-term period. Although this component is not related to the standard concept of risk, the fact that it is a risk arising from the failure of projects and has a probability of occurrence, gives this term a concept similar to the classic risk, and can be studied as such. The total risk is calculated as the sum of the two components.
The first component (R e ) addresses the environmental hazard, since it estimates the general degradation of the environment, i.e., the aquifer's pollution. Failure is defined as the presence of contaminants in the aquifer. Therefore, the corresponding probability of failure, P fe , is calculated by the probability related to the presence of pollutants' mass at the water supply wells. Since the degradation of the environment lasts as long as the plume remains in the aquifer, it makes sense that the cost of failure (C fe ) is a function of time. Cost of failure (C fe ) is a function of many factors, such as existing legislation, the state's concern for environmental protection, political will, etc. The partition of the risk into two components requires their consecutive distributivity in the utility function c(C f ). Thus, the utility function of the environmental risk is characterized as c(C fe ) and corresponds to the tendency of the decision-maker to avoid this type of risk. The first component of failure risk is formulated as: Regarding the second component of failure risk (R w ), failure has been previously defined as the presence of contamination in the supply wells. Then, the probability of failure (P fw ) is expressed as the percentage of failed aquifer remediation solutions by the total of realizations. As already mentioned, the assumption of immediate cessation of water intake excludes the relation of the risk to human health, or even life. Thus, the cost of failure (C fw ) is defined as the cost of water volume (V p ) needed for the community on a time basis. Cost of water volume (V p ) varies over time and depends on the number of contaminated water supply wells. The utility function is defined as c(C fw ). The second component of failure risk is formulated as: and the total risk is expressed by the equation: The comparison of the alternatives that will facilitate the final decision is performed within the framework of the cost-benefit-risk theory. As there is no benefit in this case, the cost function only includes the cost of each alternative. For the random alternative j, the objective function U takes the general form: where C j is the total cost of alternative j, and R j is the total risk. The obvious goal of the decision analysis method is the selection of the alternative that minimizes the objective function of Eq. 13. The total cost for such a problem is divided as follows: where C jcon is the construction cost and C jop is the operational cost of the clean-up wells used in each alternative j.

Study area
The Lake Karla watershed, with an area of 1171 km 2 , is in Eastern Thessaly of central Greece. The mountains of Pelion (S-E), Mavrovouni (E), Ossa (N-E), and Xalkodonio (S-W) form a closed and elongated basin. The Pinios River forms the northwest boundary, where the Thessaly plain exists (Fig. 2). The climate is characterized as Mediterranean with the mean annual of precipitation, temperature, and potential evapotranspiration to value at 560 mm, 14°C, and 775 mm, respectively (Sidiropoulos et al. 2015). The area is divided in two zones regarding the soil cover: the first zone is located at the Mavrovouni foothills and the second one covers the entire plain (Zalidis et al. 2004). Swo types of soils prevail in the first zone depending on the parent material of the semi-mountainous parts that surround them: the first type (60% of the zone area) prevails at the south-east, next to the area of former Lake Karla, containing Entisols and Inceptisols and the second one (40% of the zone area), located at the northeast, contains soils originating from the decomposition of limestone rocks. Two soil types are included in the second zone as well: the first is located at the west edge of the plain (31% of the zone area) characterizing by heavy soils of Mollisols and Inceptisols, which are quite fertile but small patches of alkaline soils are observed and the second one occupies the central part of the plain from north to south (69% of the zone area) containing heavy soils of Vertisols, Mollisols and Inceptisols which are alkaline and saline.
The hydrogeological regime of the study area consists of aquiferous sand intercalations split up by layers of Neogene and clay to silt-clay deposits. The Neogene deposits consist of marbles and conglomerates. Permeable structures, which appear mainly on the plain, cover 54.9% of the total area, impermeable geological structures cover 30.6%, and karstic aquifers cover 14.5%. Interest was focused on the aquifer of contemporary sediments out of all the karstic and sedimentary aquifers in the area since its resources are harvested from the boreholes. The selection was supported by the widespread belief that there is no hydraulic conduct between karstic and sedimentary aquifers, which is based on earlier hydrogeological investigations. Due to their piezometric continuity, the four zones of the hydrogeological basin of eastern Thessaly are treated as a single aquifer (Constantinidis 1978). Three different aquifer types are found in the vertical plane: the sandy-clay lake deposits, the coarse sediments, and the marbles. The three aquifers mentioned above are stacked one on top of the other (successive horizons). The top two layers were treated as one because of their high hydraulic conduct. The research was focused on these two layers since their hydraulic conduct with the aquifer of the marbles is negligible. The aquifer lies in the lower part of the basin covering an area of 500 km 2 . It is a phreatic grainy aquifer, with a geological structure made up mostly of recent grains of varied sizes and thickness up to 200 m (Mylopoulos and Sidiropoulos 2016). The basement rocks are composed of impermeable marbles and schists (Constantinidis 1978). Due to the existence of impermeable bedrocks on Mavrovouni Mountain, the eastern (N-S) aquifer edge is a noflow boundary, but the southern, western, and northern (S-N) margins have a hydraulic connection with adjacent aquifers (Fig. 2). The maximum absolute water level is found at the north and is about 30 m, while the minimum one is at the south which topically reaches the -30 m. The direction of the groundwater flow is from north-west to south-east.
Lake Karla catchment is a rural basin characterized by water-demanding crop cultivations for the last fifty years. According to the Corine Land Cover Database, over 80% of the aquifer's area is characterized as farmland. Cotton, wheat, maize, alfalfa, trees, and vegetables are the main crops. After the natural lake was drained in the 1960s, the agricultural intensification began. The cultivation of the region was one of the reasons for the lake's draining. Due to the lack of a surface water body and an irrigation network, about 97% of irrigation demands are met by groundwater through private extraction wells . Their yield varies from 20 to 100 m 3 /h. Groundwater resources are renewed by natural recharge, the irrigation return flow, and by reservoir seepage. Overexploitation of groundwater resources has resulted in a 100-m drop in the water table during the last three decades (Sidiropoulos et al. 2013). In 2018, a reservoir was built up at the basin's lowest part (with an average area of 38 km 2 ) to prevent further degradation. The construction project started in 2000, while the water supply to the reservoir began in 2012. The reservoir was built up to meet the irrigation needs of surrounding agricultural areas (about 46 hm 3 /y), although it does not yet fully operate. As a complementary project is the establishment and operation of forty additional supply boreholes which will provide the city of Volos with 15.7 hm 3 per year for covering of domestic urban water needs. The boreholes are located in the southern part of the aquifer, but they are not operational yet (Fig. 2). Sidiropoulos et al. (2015) and Mylopoulos and Sidiropoulos (2016) provide further details on the restoration project.
According to Directive 91/676/EEC and the Joint Ministerial Decision 161,890/1335/1997, the study area is one of the seven Nitrate Vulnerable Zones (NVZs) of Greece from agricultural run-offs. However, nitrogen overfertilization practices are often attempted by the farmers for crop yield maximization. The amounts of fertilizer applied in crops are greater than those recommended by the European Union's Nitrates Directive . In recent years, Nitrogen-based fertilizers are used widely in agriculture because they are the most absorbed and productive fertilizers, and they have high solubility and biodegradation. On the other hand, nitrate anions are leached easily due to their high mobility and solubility (Topcu and Kirda 2013) and along with the depletion of the groundwater table and the insufficient natural replenishment, a steady rise in nitrate concentrations has been triggered . A network of nine observation wells for groundwater quality monitoring has been established in the area (Fig. 2). The groundwater quality is characterized by nitrates presence throughout the entire study area, with the maximum concentration values (40 mg/L) observed in the central and southern part of the aquifer. The aquifer's quantitative and qualitative deterioration has caused local authorities to shut down most of the drinking water wells. Eleven of them operate for the last two decades: ten on the south side of the aquifer that supply with drinking water the Municipality of Volos, and one on the eastern part that supplies the Municipality of Rigas Fereos.

Results and discussion
The presentation and discussion of the results follow the structure of the procedure presented in the previous section of this study. But an important note should be made at this point for the time periods of applicating the procedure's steps. The simulation of all the time-variable models is conducted in a monthly time step. Both the integrated modelling and the stochastic simulation are applied for the historical period of June 1995 to January 2017. The period after January 2017 the problem of decontamination of the aquifer begins with the use of stochastic optimization and risk analysis methodologies. The meteorological, hydrological, and water use data of the management period were assumed to be the historical ones. The optimization exercise is applied to the qualitative and quantitative status of the aquifer as simulated by the MODFLOW and MT3DMS codes for January 2017. The management period lasts until the decreasing of the mean nitrate concentration values in the water supply zone below 25 mg/L (Fig. 2). The duration of the management period depends on the management scenarios described below.
Stochastic optimization and risk analysis are tested for two scenarios of nitrogen fertilizer applications in the cultivated areas: 1. Scenario A-Current nitrogen fertilizer application: it represents the real fertilizer applications. The farmers try to maximize crop production by applying greater amounts of nitrogen fertilization than those imposed by the E.U. Nitrates Directive (Table 1). 2. Scenario B-Compliance with the E.U. Nitrates Directive: it is a hypothetical scenario in which the farmers apply the amounts of nitrogen fertilization that are imposed by the E.U. Nitrates Directive requirements (Table 1).
Detailed information for the development of the models at the Lake Karla aquifer has been documented in previous studies (Sidiropoulos et al. 2019. Statistical parameters were used to evaluate the calibration and validation of the simulation models. The multistart Generalized Reduced Gradient algorithm and the split sample test were used to calibrate the six parameters of UTHBAL model using the observed monthly values of basin runoff to Pagasitikos Gulf for the period October 1960-September 2002. For model calibration and validation, the Nash-Sutcliffe Model Efficiency was equal to 0.66. The LAK3 and UTHRL models were calibrated against the reservoir's observed water levels for the operational period 2012-2017. The LAK3 model was calibrated for the ratio of the vertical hydraulic conductivity of the reservoir' bottom to the depth, which determines the leakage. The Coefficient of Determination (R 2 ) was equal to 0.78. The UTHRL model was calibrated for the water transfer from Pinios river, and the Coefficient of Determination (R 2 ) was equal to 0.84. It was not possible a validation to be achieved for these two models because of the short period of reservoir operation. The calibration of MODFLOW model was implemented for the spatially distributed hydraulic conductivity adjustment against twenty-seven

Integrated modelling
This chapter presents the application approach to the Lake Karla aquifer and the results of the simulation system's models that do not apply in a stochastic mode, the UTH-BAL, the UTHRL, the LAK3, and the GEPIC.
The UTHBAL model is approached in a semi-distributed mode since the aquifer is entirely located in the lower zone of the basin (Sidiropoulos et al. 2010). Monthly time series of recharge and surface runoff are generated as presented in Fig. 3a. The mean annual groundwater recharge volume is estimated equal to 41.7 hm 3 .
The UTHRL and LAK3 models participate in the simulation system in the period when the reservoir operates even partially. This period is from 2012 to 2017. The inflows are the precipitation to the reservoir, the natural surface runoff through 1 T ditch as simulated by the UTHBAL model, the water transfer from the Pinios river as simulated using the pumping station's operating schedule, and the water transfer system's conveyance. The outflows are the seepage to groundwater as calculated by the LAK3 model and the evaporation of surface water. No withdrawals for crop irrigation were achieved for this 5-year period since the reservoir flooding was in progress. Unfortunately, due to significant damage to the pumping stations and the uncompleted works, the water supply of the reservoir was not fully successful (Sidiropoulos et al. 2017). The combined use of the UTHRL and the LAK3 models estimates the water budget of Lake Karla (Fig. 3b). The mean annual volume of the reservoir's seepage to the aquifer is estimated to 10.64 hm 3 . GEPIC (GIS Environmental Policy Integrated Climate) is an agroecosystem simulation tool coupled with the Environmental Policy Integrated Climate (EPIC) biophysical simulation tool. GEPIC abstracts necessary data from GIS raster maps and edits them to the EPIC data format. The data relocated to the EPIC core, where the EPIC code runs. The EPIC simulation results are redirected to the GEPIC interface environment to produce output maps. The successful application of the GEPIC model utilizes information on Lake Karla aquifer location, geomorphology data such as DEM and slope, climate and land use data, plant phenology parameters, soil physical and chemical characteristics, and finally agricultural management practices such as the annual irrigation depth and the annual fertilizer application per crop. Land use data are classified into 3 categories from 0 to 2. 0 indicates the absence of a specific crop, while 1 indicates crops that meet their water needs under the rain, and 2 indicates crops that meet their water needs with irrigation conditions. GEPIC is applied in a monthly time step and spatial resolution of 0.01°9 0.01°d egrees (about 1 km 9 1 km nearby the equator) for raster datasets. During the simulation period (1995-2017) alterations in the annual crops are observed, so the raster map of GIS indicating the location of each crop differs for each year. Table 2 presents the mean percentage area of land uses and crop classification of the aquifer study area. Leaching depends on the spatially varying mass nitrogen loading according to the crop pattern, the cultivated area, the natural recharge to groundwater, and the nitrogen losses to the unsaturated zone. The cultivations' crop pattern is obtained from the Greek Integrated Management System, while the nitrogen loading is determined using literature (Roser and Ritchie 2019). The GEPIC model calculates the nitrate leaching of the spatially varying mass nitrogen loading to groundwater (Lyra et al. 2021). Figure 4 presents the nitrate leaching as a percent of the nitrogen loading mass results field by field, for two indicative years (2005,2015). Areas with no value (white color) refer to land uses where no agricultural activity exists such as water bodies, urban areas, and pastures. The values range from 0 to 50%. According to Fig. 4 high mass values of nitrate leaching, which reach 50% of the nitrogen loading mass, prevail in the aquifer area. Percentage values of nitrate leaching higher than 40% are observed at almost 88% of aquifer's cultivated area for 2009 and at 84% for 2015. This status does not only illustrate these two indicative years, but it is a phenomenon that occurs for the last thirty years.

Stochastic simulation
In the area of the Lake Karla aquifer, only fifteen scatter data measurements of hydraulic conductivity are available (Fig. 2). Since the aquifer is real and large-scaled, Fig. 3 Graphs of a monthly time series of recharge and surface runoff generated by UTHBAL model, b annual water budget of Lake Karla and its stored volume heterogeneity prevails. Thus, uncertainty arises in the groundwater flow and nitrate transport simulation. The range of the hydraulic conductivity samples is 0.936 m/d. The mean value is 0.378 m/d and the standard deviation is 0.312 m/d. Figure 5a presents the raw data histogram (in m/d). For second-order stationarity, the values were lognormally transformed with a variance of 1.00 and a mean value of 0.00. The anisotropy direction of hydraulic conductivity was investigated by calculating many experimental semi-variograms in different directions. A variogram map has been created with GSLIB for each semi-variogram. In general, the main anisotropy trend is revealed by the direction in which the variogram map elongates. However, there was no such direction on any map, thus, indicating that hydraulic conductivity has no anisotropy direction. But this is because the number of samples is too limited to determine hydraulic conductivity anisotropy (Hoseini et al. 2018). Many types of models semi-variograms were tested and the spherical one of Eq. 15 was selected because it produced the best fitting to the experimental semi-variogram and the estimation error with the minimum variance for the study area (Mylopoulos and Sidiropoulos 2016).
One hundred hydraulic conductivity simulation maps have been created with GSLIB. This number of realizations is evaluated as satisfactory for the statistical stability between log-transform hydraulic conductivity and simulated heads, due to: (i) the use of Kriging as the interpolation method, (ii) the conditioning of simulated data, and (ii) the calibration of the simulation models, . The frequency percentage histogram of the stochastically simulated hydraulic conductivity values (Fig. 5b) has the same statistical pattern as the raw data histogram (Fig. 5a) confirming the stochastic procedure's reliability and stability, as well as the correct selection of the spherical model semi-variogram as the most appropriate.
The groundwater flow was simulated with the MOD-FLOW model by creating a single layer with 12,500-cells grid (200 m 9 200 m) refined in the location of the wells. The hydraulic contact with the adjacent aquifer at west and south-west was addressed as a General-Head boundary with a value of conductance equal to 0.08 m 2 /d. At the eastern and south-eastern boundary, impermeable formations are present, and a no-flow boundary is assigned. The initial conditions were set up in January 1987 through the spatial interpolation of the observed hydraulic head values . Due to the absence of any field measurements, the specific storativity is assumed uniform for the aquifer's area and equal to 0.3 which equates to the average value between fine sand and clay (Johnson 1967). The hydraulic conductivity varies from 0.017 to 0.953 m/day and is treated as a stochastic parameter. The mean annual recharge to the aquifer is equal to 63.35 hm 3 and is the sum of the natural recharge from the surface, the irrigation return flow, and the reservoir's seepage. Mylopoulos and Sidiropoulos (2016) analytically presented the assessment of extracted groundwater from wells and its purpose of use. Specifically, for the given simulation period, the mean annual groundwater volume for domestic use was estimated at 3.1 hm 3 and for irrigation use 110.1 hm 3 .
The MODFLOW model received the 100 stochastic maps of spatially distributed hydraulic conductivity from GSLIB and produced additional equally probable realizations of groundwater flow. From the results of the groundwater model application, three important remarks are to be mentioned and discussed: (1) The non-renewable groundwater resources are under exploitation and their mean annual volume is estimated equal to 49.85 hm 3 . (2) The water table drawdown reaches up the 100 m during the simulation period (Fig. 6a). (3) The large number and the high pumping rates of the extraction wells cause the distortion of the natural direction of groundwater flow. Groundwater flow has been redirected from north-west to south-east to north-central and south-southeast.
The MT3DMS model was used to stochastically simulate the nitrate transport and dispersion, based on the 100 equally probable groundwater flow realizations and all the aquifer's characteristics from the MODFLOW model. Since low values of nitrate concentration characterize the adjacent aquifers, the boundaries were considered as zero nitrate transport. The hydrodynamic advection and dispersion parameters could not be studied as spatially variable due to the lack of scatter data. For that reason, uniform values were assumed based on the literature (Gelhar et al. 1992;Shamrukh et al. 2001;Konikow 2011). The longitudinal and transverse dispersivity and porosity coefficients were assessed equal to 20 m, 2 m, and 0.3, respectively, with the MT3DMS calibration process by using the zonal automated parameter estimation process. The initial conditions were set up for January 1995 with the spatial interpolation of the observed nitrate concentration values (Sidiropoulos et al. 2019). Nitrate leaching has been calculated by the GEPIC model. The MT3DMS model generated 100 equally probable maps of spatially distributed nitrate concentrations for the simulation period . Figure 6b presents the mean nitrate concentration values of the 100 stochastic maps for January 2017.
Nitrate concentration values higher than the 25 mg/L limit (Directive 91/676/EC) are observed in the central and southern part of the aquifer where intense water table drawdown is observed (Antonakos and Lambrakis 2000;Sidiropoulos et al. 2021). Special emphasis must be given to the fact that one of the above-mentioned locations is the area where the operating supply wells are located (continuous black line). This result raises safeguarding public health issues. The MT3DMS model indicates that the nitrate pollutants occur in places outside the supply wells area and due to the groundwater flow (Fig. 6b) are directed towards it.

Stochastic optimization
The remediation scheme focuses only on the area of the supply wells (continuous black line of Fig. 6b) where sets of clean-up wells are introduced in the optimization exercise at the boundary locations the plume crosses (Fig. 6b).
Controlling the hydraulic gradient along the boundary will prevent the nitrates' transport to the supply wells. The optimal positioning of the wells was derived from the optimization procedure, since for the unnecessary ones a zero pumping rate was set. The operational cost of the clean-up wells per unit pumping (A) was configured at 0.10 €/m 3 (Sidiropoulos 2014), and the allowable hydraulic head depletion at 30 m, according to the Environmental Terms for the construction and operation of the former lake Karla re-flooding project (Ministry of Environment, Planning and Public Works, 2000). The optimization problem was solved for each one of the 100 realizations, with the conjunctive use of GWM models, and optimum solutions were derived. These optimum solutions were imported in groundwater models to generate the optimum realizations of the aquifer's quantity and quality status and to form the predetermined set of alternatives needed for the application of the risk analysis framework. The extensive processing of the optimization results is presented for the two scenarios in the next paragraphs.

Scenario A
Twenty-nine pairs of clean-up wells were introduced in the optimization exercise, but only thirteen were chosen as operational and necessary (Fig. 6b). The management period lasted ten years (2017)(2018)(2019)(2020)(2021)(2022)(2023)(2024)(2025)(2026)(2027) and this resulted in the decontamination procedure to be applied for a whole The pumping rate of the clean-up wells varied for each realization, and the mean values of the total capacity of the thirty-three pairs range from 10.271 to 167.836 m 3 /d (Fig. 7a). The clean-up wells characterized by low pumping rates are located at the plumes' boundary, while the ones with the high rates are at the center of the plumes. Figure 7b presents the total pumping rate of the clean-up wells system for every realization. The range of the objective function varies from 369.27 * 10 3 € (Realization 15) to 374.96 * 10 3 € (Realization 72) as it is presented in Fig. 7b. This is a significant amount of money that must be spent on nitrates removal due to the long duration of the decontamination period. This result leads to a very important remark: if the farmers apply greater amounts of nitrogen fertilization than those imposed by the E.U. Nitrates Directive any nitrates removal attempt is so costly and time-consuming that it may not even be eligible.
To complete the processing of the optimization results, it was deemed necessary to check the reliability of the optimum solutions. The methodology is as follows: each optimal solution resulting from the corresponding realization was checked if it met the requirements of the optimization exercise for the remaining 99 realizations, with a reliability percentage for every optimum solution, the degree of reliability. Figure 7c presents the degree of reliability for each optimum solution. Figure 7d shows the degree of reliability in relation to the total pumping rate of the clean-up wells system. It is noted that there are solutions with a low degree of reliability equal to 34 successes out of 100. There are also solutions with a very high degree of reliability. The maximum degree of reliability observed is 96 successes, out of 100. It is no coincidence that the optimum solution of realization 15, with the lowest value of the total pumping rate, has the lowest degree of reliability, while the optimum solution of realization 72, with the highest value of the total pumping rate, has the highest degree of reliability. This is attributed to the higher the pumping rate of the wells, the more probable the successful decontamination of the aquifer. But the optimum solution of realization 72 is the most expensive alternative.

Scenario B
Twenty-nine pairs of clean-up wells were introduced in the optimization exercise, but only thirteen were chosen as operational and necessary (Fig. 6b). The management period lasted four years (2017)(2018)(2019)(2020)(2021), which is a significantly shorter duration than that of Scenario A. This is because it is assumed that the farmers applied the amounts of nitrogen fertilization that are imposed by the E.U. Nitrates Directive during the decontamination period. The pumping rate of the clean-up wells varied for each Fig. 7 Scenario A graphs of a mean values of the total capacity of the clean-up well pairs, b total pumping rate of the clean-up wells system and objective function values for every realization, c degree of reliability for each optimum solution, and d degree of reliability for total pumping rate of the clean-up wells system realization. The mean values of the total capacity of the thirty-three pairs range from 9.66 to 166.62 m 3 /d (Fig. 8a). The clean-up wells characterized by low pumping rates are located at the plumes' boundary, while the ones with the high rates are at the center of the plumes. Figure 8b presents the total pumping rate of the clean-up wells system for every realization. The range of the objective function varies from 145.85 * 10 3 € (Realization 50) to 148.36 * 10 3 € (Realization 93) as is presented in Fig. 8b. The cost of the decontamination method is less than half that of scenario A. This result shows the importance of the need for full and faithful implementation of the E.U. Nitrates Directive by farmers, otherwise, the decontamination of the aquifer will have huge financial costs. Figure 8c presents the degree of reliability for each optimum solution. Figure 8d shows the degree of reliability in relation to the total pumping rate of the clean-up wells system. The optimum solution of realization 50 has the lowest degree of reliability equal to 31%, while the optimum solution of realization 93, has the highest degree of reliability equal to 97%.
So, the manager would face the following dilemma for both scenarios: To choose the solution with the lowest cost and the lowest reliability or the one with the highest cost and the highest reliability? Optimization cannot provide a solution, and this is an important disadvantage. That is why decision analysis comes to give a solution to this dilemma of choosing the strategy with a low cost and the minimum risk of failure.

Risk analysis
The decision model components and their terms depend on the pollution problem and the environmental, economic, technical, political, and legal commitments ought to be defined. The probability of failure (P fe ) of Eq. 10 refers to the percentage of wells that abstract groundwater with nitrate concentration greater than 25 mg/L. The cost of failure (C fe ) was estimated to be 100 €/day which corresponds to a conservative view of the fines imposed in such cases, as reported by the River Basin Management Plans (Special Secretariat for Water 2017). For the estimation of the utility function c(C fe ), the assumption of a neutral attitude of the manager towards the risk was made where c(C fe ) = 1, a fact that covers the multitude of cases in Greece (Almpouras 2015). The estimation of the second component R w (t) terms is as follows. Since the corresponding probability of failure (P fw ) is expressed as the percentage of failed aquifer remediation solutions from the one hundred realizations, it is given by the equation: According to the pricing policy of the Municipal Water Supply and Sewerage Company of Volos, a conservative Fig. 8 Scenario B graphs of a mean values of the total capacity of the clean-up well pairs, b total pumping rate of the clean-up wells system and objective function values for every realization, c degree of reliability for each optimum solution, and d degree of reliability for total pumping rate of the clean-up wells system purchase price of water, is considered at 0.41 €/m 3 , and the cost of failure is given by the equation: where V p (t) is the volume of the contaminated water. The utility function c(C fw ) is equal to one according to the assumption made for the c(C fe ). For the estimation of the objective function U (Eq. 13) the assessment of the construction cost (C con ) and the operational cost (C op ) of the clean-up wells is needed (Eq. 14). The construction cost is set at zero because existing irrigation wells are used as clean-up wells, as there are plenty in the area that are programmed to stop operating (Ministry of Environment, Planning and Public Works 2000). The operational cost is equal to the objective function of the optimization exercise given by Eqs. 1, 4, and 7 and presented in Fig. 7b (Scenario A) and Fig. 8b (Scenario B). From this approach, it is clear the reason why stochastic optimization precedes risk analysis in the proposed framework. Stochastic optimization provides risk analysis with values of the objective function's (U) terms (P fw , C op ), so as the second is applied. This is the way how these two procedures link and the meaning behind the prementioned phrase that stochastic optimization provides to risk analysis a predetermined set of alternatives.

Scenario A
The percentage of wells that abstract groundwater with nitrate concentration greater than 25 mg/L varies among the realizations from eight out of ten (P fe = 0.8) to ten out of ten (P fe = 1). The fact that a decade is needed for the decontamination to be accomplished affects the values of C fe (t) and C fw (t). The cost of failure C fe (t) is the same for each realization and equal to 365,000 €. The probability of failure (P fw ) varies from 0.04 to 0.66. The cost of failure C fw (t) varies from 37,884 € to 44,772 €. So, the objective function (U) can be calculated for every optimum solution (Fig. 9a).

Scenario B
The probability of failure (P fe ) varies among the realizations from 0.5 to 0.8. The fact that four years are needed for the decontamination to be accomplished affects the value of C fe (t) and C fw (t). The cost of failure C fe (t) is the same for each realization and equal to 146,000 €. The probability of failure (P fw ) varies from 0.03 to 0.69. The cost of failure C fw (t) varies from 37,884 € to 44,772 € and the objective function (U) is calculated for every optimum solution (Fig. 9b).
Risk analysis creates larger costs than stochastic optimization as expressed by Eq. 13. This is due to the costs of failure (C fe (t), C fw (t)) that are inserted in the total risk (R j ), and express the environmental, economic, technical, political, and legal commitments. So, the objective function of the risk analysis (U) is larger than that of the stochastic optimization (F) by the term of the total risk. In stochastic optimization, the realizations with the greatest degree of reliability for Scenario A are the 87th (96%) and the 72nd (96%) but are the most expensive since F 87 = 374.94 * 10 3 € and F 72 = 374.96 * 10 3 €. In risk analysis, these two realizations have the lowest cost since U 87 = 668.45 * 10 3 € and U 72 = 668.47 * 10 3 €. The same fact exists for realizations 29 (96% degree of reliability, F 29 = 148.15 * 10 3 €, U 29 = 221.53 * 10 3 €), 91 (96% degree of reliability, F 91 = 148.28 * 10 3 €, U 91 = 221.66 * 10 3 €), and 93 (97% degree of reliability, F 93 = 148.36 * 10 3 €, U 93 = 221.64 * 10 3 €), for the scenario B. So, stochastic optimization treats the most reliable solutions as the most expensive, while risk analysis treats them as the most inexpensive. This is because the most reliable solutions have the lowest risk, and the risk is evaluated as a cost through the total risk term (R j ). The higher the degree of reliability, the lower the values of the probability of failure (P fw , P fe ). Despite the different treatment that the two procedures give to the cost, the risk analysis finally validates the usefulness of the optimization procedure as it promotes the optimum solution with the highest degree of reliability which carries the least risk. So, at this point, the risk analysis answers the dilemma of which optimum solution to be selected. Specifically, for scenario A, the optimum solution 87 is the remediation strategy with the minimum cost and the least possibility of failure. As for scenario B, the manager can choose the optimum solution 93 which has the highest degree of reliability and therefore the lowest risk of failure and is only 115 € more expensive than the 29th. Finally, risk analysis highlights to a greater degree than the stochastic optimization the total cost difference between the 2 scenarios since the longer the remediation period, the higher the costs of failure (C fe (t), C fw (t)). A remediation strategy without the simultaneous nitrogen fertilizer application reduction in cultivated areas fails because of the high cost.
The results of the two proposed solutions application for both scenarios are presented in Fig. 10. Figure 10a, b present the maps of nitrates concentration of solution 87 for scenario A on January 2027 and of solution 93 for scenario B on January 2021, respectively. The decontamination of the aquifer in the water supply wells' area is successful for both scenarios since the nitrate concentration values are lower than 25 mg/L. But, for scenario A, since large application of nitrogen amounts take place, this limit is exceeded in areas of the aquifer such as in the center. On the contrary, no exceedance of the limit of 25 mg/L is observed for the entire study area for scenario B.

Conclusion
The mismanaged groundwater exploitation and the longterm and intense fertilization practices in regions with significant agricultural production have led to quantitative and qualitative degradation of aquifers. When groundwater is used for urban drinking purposes, aquifers' remediation is mandatory for protecting public health from a variety of diseases. But the simulation and management of real, largescale aquifers encounter uncertainty mainly generated by the heterogeneity of hydraulic conductivity. Under this status, the decision-making must be focused on choosing the optimum remediation solution from several alternatives that meet all environmental, economic, technical, political, and legal commitments.
In this study, a risk analysis framework is proposed for the optimum nitrate removal of degraded aquifers under hydraulic conductivity uncertainty. The framework is based on integrated modelling, stochastic simulation, stochastic optimization, and risk analysis approaches, which are implemented through the coupling of geostatistical, simulation, and management models. The selected methodology for nitrates removal is the stabilization of the plume with the use of pumping well pairs aimed at controlling the hydraulic gradient along the contaminated plume boundary. The final target of the proposed framework is the selection of an aquifer decontamination strategy with the least risk of failure and the minimum cost from several alternatives. The coupling of stochastic optimization and risk analysis has been proven successful for the nitrate removal of an aquifer. The application of these two approaches is achieved in a way that one methodology strengthens the weakness of the other. Stochastic optimization must proceed risk analysis as it provides a predetermined set of alternatives, and the values of the objective function's (U) terms (P fw , C op ), so as the second to be applied. Stochastic optimization treats the most reliable solutions as the most expensive, while risk analysis treats them as the most inexpensive. Risk analysis based on the optimum solutions of stochastic optimization gives a solution to the dilemma of choosing the strategy with a low cost and the minimum risk of failure. Nevertheless, any remediation strategy without the simultaneous nitrogen fertilizer application reduction in cultivated areas could be unprofitable. This conclusion highlights the importance of the need for full and faithful implementation of the E.U. Nitrates Directive by farmers and local stakeholders.
The limitations of the proposed methodology have to do with the application of the modelling system since all its components ought to be set up and calibrated for every study area. This requires time and effort and the availability of data. Furthermore, the use of the modelling system components depends on the type of water resources of an area. For example, the limitations of the modelling system are: (a) MODFLOW model cannot simulate karstic processes, (b) GEPIC model is used for rural basins, (c) LAK3 and UTHRL models are used for reservoirs. However, since the linking of the models is based on outputs/inputs interconnections of the models, there is the option components (i.e. models) to be removed or added with the necessary coupling. Thus, the work presented in this paper could be more general and not site specific.
Nevertheless, the successful application of the proposed risk analysis framework on the Lake Karla aquifer proves that it could become an effective tool for determining that nitrate removal strategy with the less risk of failure and the minimum cost and it could be implemented for rural basin degraded aquifers.
Acknowledgements Dr. Pantelis Sidiropoulos was a post-doctoral scholar of the University of Thessaly. The scientific publication was held within the framework of the invitation ''Granting of scholarship for Post-Doctoral Research'' of the University of Thessaly, which was being implemented by the University of Thessaly and was funded by the Stavros Niarchos Foundation.
Funding Open access funding provided by HEAL-Link Greece. The research was funded by the Stavros Niarchos Foundation.
Data availability The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Declarations
Conflict of interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.